Where is this function coming from?

[Alex1993]

Hello!

I found this function online and I cannot trace the research paper from which is coming from. It calculates the Cumulative Normal Distribution. Any clue where I should look? I heard someone mentioning "Joshi", but so far I had no success to see where this is coming from and where are the constants are coming from.

Thank you!

Numeric CND(const Numeric& d)//there are many functions for generating a random number from a Gaussian Distribution. BoxMuller is one and this is a 
							 //second function that generates random numbers
{
	const double A1 = 0.31938153; const double A2 = -0.356563782;//defines constant
	const double A3 = 1.781477937; const double A4 = -1.821255978;
	const double A5 = 1.330274429;
	const double RSQRT2PI = 0.39894228040143267793994605993438;
	double tmp = 1.0 / (1.0 + 0.2316419 * fabs(d));
	double cnd = RSQRT2PI*exp(-0.5*d*d)*(tmp*(A1 + tmp*(A2 + tmp*(A3 + tmp*(A4 + tmp*A5)))));
	if (d > 0) cnd = 1.0 - cnd;
	return cnd;
}

[mbozzi]

Joshi is an author, IIRC:
A quick Google confirms that:
https://www.amazon.com/Patterns-Derivatives-Pricing-Mathematics-Finance/dp/0521721628

Seems to be relevant (you'll find a little bit of statistics in finance.)

I don't have the book, though, so I cannot definitively provide a source for the code you've provided.

Edit, re. the constants (just speculation).
RSQRT2PI is the reciprocal square root of 2π, i.e., (1 / sqrt(2π)).

The Phi function (the CND) is defined as (1 / sqrt(2π)) * the integral from -infinity to x of exp(-x^2 / 2), or obviously the integral of the standard normal distribution.

See:
https://en.wikipedia.org/wiki/Normal_distribution#Numerical_approximations_for_the_normal_CDF

The A1, A2, A3, ... are oscillating.
I will note that the form (tmp * (A1 + tmp * (A2 + tmp * ( ... )))
Is a right-fold (or reduction) expression of the function λx, y: x + tmp*y over A1..A5.

That looks like a continued fraction to me, given that tmp is fractional. Perhaps look into continued fraction approximations to the CND. Maybe that's where A1..5 come from.

See this paper:
http://qspace.qu.edu.qa/bitstream/handle/10576/6929/049506-0009-fulltext.pdf;sequence=8

Or maybe someone will share more than I can about this. My math & statistics are pretty weak.

[lastchance]

The cumulative density function of the normal distribution (CND, or Phi) is defined by @mbozzi as above.

Instead of your function you could compute it from the C++ library using the error function erf(x), which was apparently introduced in the C++11 standard (http://www.cplusplus.com/reference/cmath/erf/). The relationship is

CND = (1/2)( 1 + erf( x / sqrt(2) ) )

The approximation almost certainly comes from an original approximation for the error function given in Abramowitz and Stegun (their equation 7.1.26) - see also https://en.wikipedia.org/wiki/Error_function
(Note that to compare the constants you will have to play with lots of factors of sqrt(2.pi), sqrt(2); however, they do match.)

Abramowitz and Stegun cite the original reference as:
C.Hastings, Jr, Approximations for digital computers, Princeton University Press, 1955 - I believe this has recently been reprinted. You may find other information on possible original sources at http://dlmf.nist.gov/7.24

You might also find approximations like this (in C++) in
Numerical Recipes: The Art of Scientific Computing (3rd ed.), New York: Cambridge University Press, 2007 (but I don't have a copy to hand).

Abramowitz and Stegun can be found online at
http://www.cs.bham.ac.uk/~aps/research/projects/as/book.php
It is still a treasure trove of numerical methods for special functions.

I'm not sure that too much can be read into the form
tmp*(A1 + tmp*(A2 + tmp*(A3 + tmp*(A4 + tmp*A5))))
It is simply an efficient way of evaluating polynomials like
a₁t+a₂t²+a₃t³+a₄t⁴+a₅t⁵

EDIT. The precise numbers that you are using come from "Options, futures, and other derivatives" by John Hull - many editions - who has a section on the Cumulative Standard Normal Function. However, he cites Abramowicz and Stegun.