## Inverse Normal CDF

Now that I’ve got some Monte Carlo code up, it’s inevitable that I will eventually need an implementation of the Inverse of the Normal Cumulative Density Function (CDF). The inverse of a CDF is called a Quantile function by the way, so I’ll often refer to this as the Normal Quantile function. The reason this is so important is because, as I discussed in my post on Random Number Generators, it can be used to convert uniformly distributed random numbers into normally distributed random numbers, which we will need for Monte Carlo simulations if we believe that underlyings are distributed according to a lognormal distribution.

As a reminder, I’ve repeated the graph from the previous post here to show the behaviour of the Normal Quantile:

As with the implementation of an algorithm for the Normal CDF (discussed here), there were several possibilities for implementation. In the end I’ve decided to go with the analytic approximation due to Beasley, Springer and Moro as presented in Joshi’s text on computational derivative pricing:

a1 = 2.50662823884
a2 = -18.61500062529
a3 = 41.39119773534
a4 = -25.44106049637

b1 = -8.47351093090
b2 = 23.08336743743
b3 = -21.06224101826
b4 = 3.13082909833

c1 = 0.3374754822726147
c2 = 0.9761690190917186
c3 = 0.1607979714918209
c4 = 0.0276438810333863
c5 = 0.0038405729373609
c6= 0.0003951896511919
c7 = 0.0000321767881768
c8 = 0.0000002888167364
c9 = 0.0000003960315187

A polynomial form is used for the central region of the quantile, where $\inline&space;0.08&space;<&space;x&space;<&space;0.92$:

$y&space;=&space;x&space;-&space;0.5&space;\&space;;&space;\quad&space;z&space;=&space;y^2$

$\Phi^{-1}(x)&space;\simeq&space;y&space;\cdot&space;{&space;(a_1&space;+&space;a_2&space;z&space;+&space;a_3&space;z^2&space;+&space;a_4&space;z^3&space;)&space;\over&space;(1&space;+&space;b_1&space;z&space;+&space;b_2&space;z^2&space;+&space;b_3&space;z^3&space;+&space;b_4&space;z^4&space;)&space;}$

And if $\inline&space;x&space;\leq&space;0.08$ or $\inline&space;x&space;\geq&space;0.92$ then:

$y&space;=&space;x&space;-&space;0.5&space;\&space;;&space;\quad&space;z&space;=&space;\Bigl\{&space;\begin{matrix}&space;x&space;&&space;y\leq&space;0\\&space;1&space;-&space;x&space;&&space;y&space;>&space;0&space;\end{matrix}&space;\&space;;&space;\quad&space;\kappa&space;=&space;\ln&space;(-\ln&space;z)$

$\Phi^{-1}(x)&space;=&space;\pm&space;(c_1&space;+&space;c_2\kappa&space;+&space;c_3\kappa^2&space;+&space;c_4\kappa^3&space;+&space;c_5\kappa^4&space;+&space;c_6\kappa^5&space;+&space;c_7\kappa^6&space;+&space;c_8\kappa^7&space;+&space;c_9\kappa^8)$

With a plus at the front if $\inline&space;y&space;\geq&space;0$ and a minus in front if $\inline&space;y&space;<&space;0$.

This is fairly swift as an algorithm, and the quantile method of generating normal variates will turn out to have some key advantages over the various other methods discussed before which will become apparent in future. In a future post I’ll show how to integrate this method of converting uniform variates to gaussian variates into the Monte Carlo code that I’ve made available here, thanks to the factory pattern that I’ve used it will turn out to be incredibly straight-forward!

Finally, if you are wondering why the numbers in these numerical algorithms seem so arbitrary, it’s because they [sort-of] are! Once I’ve implemented some optimisation routines (this might not be for a while…) I’ll have a look at how we can find our own expressions for these functions, and hopefully be able to back out some of the numbers above!!