options - How to derive the implied probability distribution from B-S volatilities?

The general problem I have is visualization of the implied distribution of returns of a currency pair.

I usually use QQplots for historical returns, so for example versus the normal distribution:

Now I would like to see the same QQplot, but for implied returns given a set of implied BS volatilities, for example here are the surfaces:

USDZAR  1month  3month  6month  12month 2year
10dPut  15.82   14.59   14.51   14.50   15.25
25dPut  16.36   15.33   15.27   15.17   15.66
ATMoney 17.78   17.01   16.94   16.85   17.36
25dCall 20.34   20.06   20.24   20.38   20.88
10dCall 22.53   22.65   23.39   24.23   24.84

EURPLN  1month  3month  6month  12month 2year
10dPut  9.10    9.06    9.10    9.43    9.53

25dPut  9.74    9.54    9.51    9.68    9.77
ATMoney 10.89   10.75   10.78   10.92   11.09
25dCall 12.83   12.92   13.22   13.55   13.68
10dCall 14.44   15.08   15.57   16.16   16.34

EURUSD  1month  3month  6month  12month 2year
10dPut  19.13   19.43   19.61   19.59   18.90
25dPut  16.82   16.71   16.67   16.49   15.90
AtMoney 14.77   14.52   14.44   14.22   13.87
25dCall 13.56   13.30   13.23   13.04   12.85

10dCall 12.85   12.85   12.90   12.89   12.78

Anybody know how I could go about doing this? Any R packages or hints where to start on this? It doesn't necessarily have to be a qqplot, it could just be a plot of the density function; that would help me too. Thanks.

Answer

You can directly imply a probability distribution from a volatility skew.

Note that, for any terminal probability distribution $p(S)$ at tenor $T$, we have the model-free formula for the call price $C(K)$ as a function of strike $K$

\begin{equation} C=e^{-rT} \int_0^\infty (S-K)^+ p(S) dS \end{equation}

Therefore we can write

\begin{equation} e^{rT} \frac{\partial C}{\partial K}=\int_K^\infty (-1) \cdot p(S) dS \end{equation} and by the fundamental theorem of calculus \begin{equation} e^{rT} \frac{\partial^2 C}{\partial K^2} = p(K) \end{equation}

Therefore, all you need, in order to find the value of $p(x)$ for any $x$, is the second derivative of call prices at strike $x$.

Usually, one uses a fitted skew (such as a polynomial fit) to the available volatility values at the given tenor. In your case, with just 5 points, I would recommend fitting a parabola in log strike space. Once you have a continuous skew $\sigma(K)$ then you just need to find

\begin{equation} {\left. \frac{\partial^2 }{\partial x^2}\right|} BS_{\text{Call}}(S_0, x, \sigma(x), r, T, q) \end{equation}

evaluated at $x=K$ which can be done either with a bunch of symbol-jiggling, or by simply finite differencing. In your case I recommend the latter.

Once you have probability distribution values, of course, the process of generating the qq plots is one you have already mastered.

Edit: Sign error correction, per @Robino