You cannot derive the probability distribution you require, because for VaR you need a real-world probability distribution. From the options prices, it is only possible to obtain a risk-neutral distribution.

Now, if you are willing to assume some kind of parametric relationship between the risk-neutral and real-world distributions, then you might find the options prices useful. The resulting mathematics for a stochastic volatility model is somewhat tricky, however. You can find most of it in Jim Gatheral's books. A sloppy treatment would just take the risk-neutral distribution and shift its mean.

Obtaining the approximate risk-neutral distribution is fairly simple. Let p(S) be the time-T risk-neutral probability density. Then we see that (TeX notation alert)

• C := Call(T) = B(0,T) \int_0^\infty Max(0,S-K) p(S) dS
• dC/dK = B(0,T) \int_0^\infty 1{S>=K} (-1) p(S) dS [differentiate under integral]
• dC/dK = B(0,T) \int_K^\infty (-1) p(S) dS
• d^2C/dK^2 = B(0,T) p(K) [Fundamental thm of calculus]

Alternatively, you could say that p(S) is the density, and is the derivative of the cumulative distribution function P(S), and write

• C := Call(T) = B(0,T) \int_0^\infty Max(0,S-K) p(S) dS
• dC/dK = B(0,T) \int_0^\infty 1{S>=K} (-1) p(S) dS [differentiate under integral]
• dC/dK = B(0,T) \int_K^\infty (-1) p(S) dS
• dC/dK = B(0,T) (-1) ( P(\infty) - P(K))
• d^2C/dK^2 = B(0,T) p(K)

Either way you end up finding the density p(x) = (1/B(0,T)) d^2C(x)/dx^2 where x is the strike. So an approximate density comes from using the actual option prices available to you. You can spline interpolate, or if you have a regular grid of strikes spaced by dK you can make a histogram of values

  ( C(K+dK) - 2C(K) +C(K-dK) ) / dK^2

and divide by the discount factor to find your risk-neutral distribution.

You cannot derive the probability distribution you require, because for VaR you need a real-world probability distribution. From the options prices, it is only possible to obtain a risk-neutral distribution.

Now, if you are willing to assume some kind of parametric relationship between the risk-neutral and real-world distributions, then you might find the options prices useful. The resulting mathematics for a stochastic volatility model is somewhat tricky, however. You can find most of it in Jim Gatheral's books. A sloppy treatment would just take the risk-neutral distribution and shift its mean.

Obtaining the approximate risk-neutral distribution is fairly simple. Let p(S) be the time-T risk-neutral probability density. Then we see that (TeX notation alert) \begin{equation} C := Call(T) = B(0,T) \int_0^\infty Max(0,S-K) p(S) dS \end{equation} \begin{equation} \frac{dC}{dK} = B(0,T) \int_0^\infty 1[S>=K] (-1) p(S) dS \qquad\text{[differentiate under integral] } \end{equation} \begin{equation} \frac{dC}{dK} = B(0,T) \int_K^\infty (-1) p(S) dS \end{equation} \begin{equation} \frac{d^2C}{dK^2} = B(0,T) p(K) \qquad \text{ [Fundamental thm of calculus]} \end{equation} Alternatively, you could say that p(S) is the density, and is the derivative of the cumulative distribution function P(S), and write

\begin{equation} C := Call(T) = B(0,T) \int_0^\infty Max(0,S-K) p(S) dS \end{equation} \begin{equation} \frac{dC}{dK} = B(0,T) \int_0^\infty 1[S>=K] (-1) p(S) dS \qquad\text{[differentiate under integral] } \end{equation} \begin{equation} \frac{dC}{dK} = B(0,T) \int_K^\infty (-1) p(S) dS \end{equation} \begin{equation} \frac{dC}{dK} = B(0,T) (-1) ( P(\infty) - P(K)) \end{equation} \begin{equation} \frac{d^2C}{dK^2} = B(0,T) p(K) \end{equation}

Either way you end up finding the density

\begin{equation} p(x) = \frac{1}{B(0,T)} \frac{d^2C(x)}{dx^2} \end{equation} where $x$ is the strike. So an approximate density comes from using the actual option prices available to you. You can spline interpolate, or if you have a regular grid of strikes spaced by dK you can make a histogram of values \begin{equation} \frac{ C(K+dK) - 2C(K) +C(K-dK) }{ dK^2} \end{equation} and divide by the discount factor to find your risk-neutral distribution.