MATH60082
Examples Sheet 1
Analytic Pricing
Contents
1.1 Coding up the Cumulative Normal Distribution . . . . . . . . . . . . . . . . . . . . 1
1.1.(i) Using a polynomial approximation . . . . . . . . . . . . . . . . . . . . . . . 1
1.1.(ii) Using integration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.1.(iii)Using a higher order polynomial . . . . . . . . . . . . . . . . . . . . . . . . 8
1.1.(iv)Using a inbuilt functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.1.(v) Comparing efficiency . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
1.2 Black Scholes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
1.3 Coursework Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
1.1 Coding up the Cumulative Normal Distribution
Write a function to calculate N(x) the cumulative normal distribution with mean
zero and variance one given by the formula
N(x) =
1√
2pi
∫ x
−∞
e−t
2/2dt
1.1.(i) Using a polynomial approximation
Using a polynomial approximation to N(x) (less than single precision!):
N(x) = 1−N ′(x)(a1k + a2k2 + a3k3), x ≥ 0,
N(x) = 1−N(−x), x < 0.
Here k = 11+γx , γ = 0.33267, a1 = 0.43618, a2 = −0.12017, and a3 = 0.93730.
First, start with a new project and an empty cpp file.
#inc lude
#inc lude
#inc lude
us ing namespace std ;
1
i n t main ( )
{
// do nothing
return 0 ;
}
download
Compile and check it runs. Looking at the algorithm in the question we see that we need several
variables, k, γ, a1, a2 and a3. Declare all those variables along with a default value of x = 1.
#inc lude
#inc lude
#inc lude
us ing namespace std ;
i n t main ( )
{
double x=1;
double gamma=0.33267 , a1=0.43618 , a2=−0.12017 , a3=0.93730;
double k=1./(1.+gamma∗x ) ;
// do nothing
return 0 ;
}
download
Compile and check it runs. Now enter the calculation for N(x) with x ≥ 0 and check it works.
#inc lude
#inc lude
#inc lude
us ing namespace std ;
i n t main ( )
{
double x=1;
double gamma=0.33267 , a1=0.43618 , a2=−0.12017 , a3=0.93730;
double k=1./(1.+gamma∗x ) ;
cout << 1.−1./ sq r t ( 8 .∗ atan ( 1 . ) ) ∗exp(−x∗x /2 . ) ∗( a1∗k+a2∗k∗k+a3∗k∗k∗k ) << endl ;
// do nothing
return 0 ;
}
download
Compile and check the output, it should be 0.841352. Now you need to make it work when x < 0,
this means using an textttif statement and changing the value of γ as well as the formula. Check
your code works in both cases.
#inc lude
#inc lude
#inc lude
us ing namespace std ;
2
i n t main ( )
{
double x=−1;
double gamma=0.33267 , a1=0.43618 , a2=−0.12017 , a3=0.93730;
i f (x>=0)
{
double k=1./(1.+gamma∗x ) ;
cout << 1.−1./ sq r t ( 8 .∗ atan ( 1 . ) ) ∗exp(−x∗x /2 . ) ∗( a1∗k+a2∗k∗k+a3∗k∗k∗k ) << endl ;
}
e l s e
{
double k=1./(1.−gamma∗x ) ;
cout << 1 ./ sq r t ( 8 .∗ atan ( 1 . ) ) ∗exp(−x∗x /2 . ) ∗( a1∗k+a2∗k∗k+a3∗k∗k∗k ) << endl ;
}
// do nothing
return 0 ;
}
download
Finally, move your code into a function and test the result. You will need to use the return
keyword to return the value of the function instead of printing to screen.
#inc lude
#inc lude
#inc lude
us ing namespace std ;
double no rma lD i s t r ibut i on po ly ( double x )
{
s t a t i c double one ove r sq r t 2 =1./ sq r t ( 8 .∗ atan ( 1 . ) ) ;
double gamma=0.33267 , a1=0.43618 , a2=−0.12017 , a3=0.93730;
i f (x>=0)
{
double k=1./(1.+gamma∗x ) ;
r e turn 1.− one ove r sq r t 2 ∗exp(−x∗x /2 . ) ∗( a1∗k+a2∗k∗k+a3∗k∗k∗k ) ;
}
e l s e
{
double k=1./(1.−gamma∗x ) ;
r e turn one ove r sq r t 2 ∗exp(−x∗x /2 . ) ∗( a1∗k+a2∗k∗k+a3∗k∗k∗k ) ;
}
}
i n t main ( )
{
f o r ( i n t x=−3;x<=3;x++)
cout << ”N( ”<// do nothing
return 0 ;
}
download
The output from this code is:
N(-3) = 0.00135489
3
N(-2) = 0.0227587
N(-1) = 0.158648
N(0) = 0.500002
N(1) = 0.841352
N(2) = 0.977241
N(3) = 0.998645
1.1.(ii) Using integration
Use Simpson’s rule to evaluate the integral defining N(x). This states that when there are values
of y = f(x) for values of x at equal intervals h apart, an approximate numerical integration is
given by ∫ b
x=a
f(x)dx ≈ h
3
[y0 + 4y1 + 2y2 + 4y3 + 2y4 + . . .+ 4yn−1 + yn]
where h = x1 − x0 = (b − a)/n. Note this formula involves an even number of intervals, and
the (truncation) error diminishes as h4. You will need around 2000 points or more for machine
accuracy, making this a very expensive method.
Again, first we think about the stages of our program. We must choose a method to implement
the integration, test and verify it. Once it is then applied to the problem we need to decide how
to deal with the concept of infinity in this setting. First let use code up the integration method.
To be able to check the method, you should choose a simple function with a known integral so
that we can compare accuracy. Start with a new project and an empty cpp file.
#inc lude
#inc lude
us ing namespace std ;
i n t main ( )
{
// do nothing
return 0 ;
}
download
Compile and check it runs. We are going to implement Simpson’s method for integration firstly
on the sin function so that we can check the results. The method is given by
∫ b
a
f(x)dx ≈ h
3
[
f(a) + 2
n/2−1∑
j=1
f(a+ 2jh)
+4
n/2∑
j=1
f(a+ (2j − 1)h) + f(b)
] (1)
Start by declaring the input parameters required by the method, including the number of
steps and the integration range, and follow with the variables used inside such as the dummy
4
integration variable s and the step size h. We should supply default values and initialise variables
if required
// number o f s t ep s
i n t N=10;
// range o f i n t e g r a t i o n
double a=0,b=1. ;
// l o c a l v a r i a b l e s
double s , h , sum=0. ;
// i n i a l i s e the v a r i a b l e s
h=(b−a ) /N;
download
Next we need to add the first few terms and the last one to sum, before adding a loop to go over
the rest.
// add in the f i r s t few terms
sum = sum + s in ( a ) + 4 .∗ s i n ( a+h) ;
// and the l a s t one
sum = sum + s in (b) ;
download
Note that we should really check that N is even as this method relies on that fact. Since we have
already taken care of the terms for 0, 1 and N, the loop must run over the terms 2 up to and
including N-1. Check your loop runs over the correct values before continuing
// loop over terms 2 up to N−1
f o r ( i n t i =1; i{
s = a + 2∗ i ∗h ;
// check even terms
cout << 2∗ i << ” ” << s << endl ;
s = s + h ;
// check odd terms
cout << 2∗ i+1 << ” ” << s << endl ;
}
download
Now add them in and check against the analytical solution.
i n t main ( )
{
// number o f s t ep s
i n t N=10;
// range o f i n t e g r a t i o n
double a=0,b=1. ;
// l o c a l v a r i a b l e s
double s , h , sum=0. ;
// i n i a l i s e the v a r i a b l e s
h=(b−a ) /N;
// add in the f i r s t few terms
5
sum = sum + s in ( a ) + 4 .∗ s i n ( a+h) ;
// and the l a s t one
sum = sum + s in (b) ;
// loop over terms 2 up to N−1
f o r ( i n t i =1; i{
s = a + 2∗ i ∗h ;
sum = sum + 2.∗ s i n ( s ) ;
s = s + h ;
sum = sum + 4.∗ s i n ( s ) ;
}
// complete the i n t e g r a l
sum = h∗sum/ 3 . ;
// output r e s u l t s
cout << sum << ” ” << 1.− cos (b) << endl ;
r e turn 0 ;
}
download
You should check accuracy by varying the value of N.
Next we are going to change this algorithm to solve the normal distribution. Do this in 2
steps, first go through and change all calls to the sin function by exp(-s*s/2.), then multiply the
final result by 1/sqrt(2*pi). The result (for a=0 and b=1) should be 0.341345. Now we have to
deal with infinity, so how can we integrate over infinity? For the normal distribution this is quite
easy as it is an even function and we know that the integral between 0 and infinity is one half.
So simply adding one half to the result of the integration between 0 and b gives the value of N(b)
for all b. At this stage your code should look like
// add in the f i r s t few terms
sum = sum + exp(−a∗a /2 . ) + 4 .∗ exp(−(a+h) ∗( a+h) /2 . ) ;
// and the l a s t one
sum = sum + exp(−b∗b /2 . ) ;
// loop over terms 2 up to N−1
f o r ( i n t i =1; i{
s = a + 2∗ i ∗h ;
sum = sum + 2.∗ exp(−s ∗ s / 2 . ) ;
s = s + h ;
sum = sum + 4.∗ exp(−s ∗ s / 2 . ) ;
}
// complete the i n t e g r a l
sum = 0.5 + h∗sum/3 ./ sq r t ( 8 .∗ atan ( 1 . ) ) ;
// output r e s u l t s
cout << sum << endl ;
download
Now move this algorithm into a function with the definition
double normalDi s t r ibut ion ( double x )
download
6
The last thing to take account of is what happens when x is either very large and positive or very
large and negative. The best thing to do here is use analytical results and return 0 if x is large
and negative and 1 if x if large and positive. Don’t forget to set the value of b to be x and choose
a value of N sufficiently large to maintain accuracy.
#inc lude
#inc lude
us ing namespace std ;
double normalDi s t r ibut ion ( double x )
{
i f (x<−10.) re turn 0 . ;
i f (x>10.) re turn 1 . ;
// number o f s t ep s
i n t N=2000;
// range o f i n t e g r a t i o n
double a=0,b=x ;
// l o c a l v a r i a b l e s
double s , h , sum=0. ;
// i n i a l i s e the v a r i a b l e s
h=(b−a ) /N;
// add in the f i r s t few terms
sum = sum + exp(−a∗a /2 . ) + 4 .∗ exp(−(a+h) ∗( a+h) /2 . ) ;
// and the l a s t one
sum = sum + exp(−b∗b /2 . ) ;
// loop over terms 2 up to N−1
f o r ( i n t i =1; i{
s = a + 2∗ i ∗h ;
sum = sum + 2.∗ exp(−s ∗ s / 2 . ) ;
s = s + h ;
sum = sum + 4.∗ exp(−s ∗ s / 2 . ) ;
}
// complete the i n t e g r a l
sum = 0.5 + h∗sum/3 ./ sq r t ( 8 .∗ atan ( 1 . ) ) ;
// re turn r e s u l t
re turn sum ;
}
i n t main ( )
{
cout . p r e c i s i o n (16) ;
cout << ”N(0)=” << normalDi s t r ibut ion ( 0 . ) << endl ;
cout << ”N(1)=” << normalDi s t r ibut ion ( 1 . ) << endl ;
cout << ”N(2)=” << normalDi s t r ibut ion ( 2 . ) << endl ;
r e turn 0 ;
}
download
Here we needed to choose N=2000 to get close to double precision, meaning that a single
function call is incredibly expensive. There are other, more efficient approximations out there.
If you are using a compiler that is math library C99 compliant you will have access to the erfc
7
function, which will give you double precision accuracy using the formula
N(x) =
1
2
erfc
(
− x√
2pi
)
1.1.(iii) Using a higher order polynomial
Using a fast and accurate polynomial approximation to N(x) (a double precision approximation,
see West, 2005, for more details):
N(x) = 1−
√
2piN ′(x)
P (x)
Q(x)
, 0 ≤ x ≤ 10√
2
,
N(x) = 1−N ′(x) 1
F4(x)
,
10√
2
< x ≤ 37,
N(x) = 1, x > 37,
N(x) = 1−N(−x), x < 0.
Here
P (x) =
6∑
i=0
aix
i, Q(x) =
7∑
i=0
bix
i
and the coefficients are:
a0 220.206867912376 b0 440.413735824752
a1 221.213596169931 b1 793.826512519948
a2 112.079291497871 b2 637.333633378831
a3 33.912866078383 b3 296.564248779674
a4 6.37396220353165 b4 86.7807322029461
a5 0.700383064443688 b5 16.064177579207
a6 0.0352624965998911 b6 1.75566716318264
b7 0.0883883476483184
For the function F4(x) use the recurrence relation
F0(x) = x+
13
20
Fi(x) = x+
5− i
Fi−1(x)
for i = 1, 2, 3, 4
Unfortunately if you are using Visual Studio it is not yet implemented (might yet be in VS2013)
so we need another method. One such method is described in West (2005) and they give a VB
code implementation of the method. We are going to follow that method closely. Start with your
empty program, I have included some benchmark values in comments to check against
8
#inc lude
#inc lude
us ing namespace std ;
i n t main ( )
{
// Benchmark va lue s o f N(x )
// N(0) = 0 .5
// N(1) = 0.841344746068543
// N(2) = 0.977249868051821
return 0 ;
}
download
Compile and check the code runs with no output.
The algorithm is given by the following expressions
N(x) = 1−
√
2piN ′(x)
P (x)
Q(x)
, 0 ≤ x ≤ 10√
2
,
N(x) = 1−N ′(x) 1
F (x)
,
10√
2
< x ≤ 37,
N(x) = 1, x > 37,
N(x) = 1−N(−x), x < 0.
where P , Q and F are polynomials, more details can be found in the examples sheet (click here).
First declare and assign some of the variables needed for the calculation
// c a l c u l a t e \ s q r t {2\ pi } upfront once
double RT2PI = sq r t ( 4 . 0∗ acos ( 0 . 0 ) ) ;
// c a l c u l a t e 10/\ s q r t {2} upfront once
double SPLIT = 10 ./ sq r t (2 ) ;
// ente r c o e f f i c i e n t s f o r the polynomia l s P and Q
double a [ ] =
{220 .206867912376 ,221 .213596169931 ,112 .079291497871 ,33 .912866078383 ,6 .37396220353165 ,0 .700383064443688 ,3 .52624965998911
e−02};
double b [ ] =
{440 .413735824752 ,793 .826512519948 ,637 .333633378831 ,296 .564248779674 ,86 .7807322029461 ,16 .064177579207 ,1 .75566716318264 ,8 .83883476483184
e−02};
// c a l c u l a t e the value o f N at x=1
double x=1. ;
download
where we have chosen x = 1 to test. Since we know what x is we don’t need to consider the
3 cases. Now calculate the polynomial expressions using a nested multiplication (important for
accuracy) and we get
double NDash = exp(−x∗x /2 . 0 ) /RT2PI ;
double Px = ( ( ( ( ( a [ 6 ] ∗ x + a [ 5 ] ) ∗x + a [ 4 ] ) ∗x + a [ 3 ] ) ∗x + a [ 2 ] ) ∗x + a [ 1 ] ) ∗x + a [ 0 ] ;
9
double Qx = ( ( ( ( ( ( b [ 7 ] ∗ x + b [ 6 ] ) ∗x + b [ 5 ] ) ∗x + b [ 4 ] ) ∗x + b [ 3 ] ) ∗x + b [ 2 ] ) ∗x + b
[ 1 ] ) ∗x + b [ 0 ] ;
cout << Px << ” ” << Qx << ” ” << NDash << endl ;
// output i s 594 .522 2272.83 0.241971
download
the output you should expect is shown. We can now estimate the value of N(x)
cout . p r e c i s i o n (16) ;
cout << 1 − exp(−x∗x /2 . 0 ) ∗Px/Qx << endl ;
// output i s 0.841344746068543
download
You now need to include the different cases, such as negative x, etc, and then move everything
into a function. After some testing you should arrive at an algorithm that looks something like:
#inc lude
#inc lude
us ing namespace std ;
/∗
∗ A ∗ ∗ f a s t ∗ and accurate polynomial approximation to $N(x ) $ ( double p r e c i s i o n
) :
∗ $$ N(x )= 1−\ s q r t {2\ pi }N ’( x ) \ f r a c {P(x ) }{Q(x) } , \quad 0 \ l e q x\ l e q \ f r a c {10}{\
s q r t {2}} , $$
∗ $$ N(x )= 1−N ’( x ) \ f r a c {1}{F 4 (x ) } , \quad \ f r a c {10}{\ s q r t {2}} < x\ l e q 37 , $$
∗ $$ N(x )= 1 , \quad x > 37 , $$
∗ $$N(x )=1−N(−x ) ,\quad x<0.$$
∗ Here
∗ $$
∗ P(x ) = \sum { i =0}ˆ6 a i xˆ i , \quad \quad Q(x ) = \sum { i =0}ˆ7 b i xˆ i
∗ $$
∗ and the c o e f f i c i e n t s are :
∗ \begin { tabu la r }{ l c | | l c }
∗ \ h l i n e
∗ $a 0$ & 220.206867912376& $b 0$&440.413735824752 \\
∗ $a 1$ & 221.213596169931& $b 1$&793.826512519948 \\
∗ $a 2$ & 112.079291497871& $b 2$& 637.333633378831\\
∗ $a 3$ & 33.912866078383& $b 3$&296.564248779674 \\
∗ $a 4$ & 6.37396220353165& $b 4$&86.7807322029461 \\
∗ $a 5$ & 0.700383064443688& $b 5$& 16.064177579207\\
∗ $a 6$ & 0.0352624965998911& $b 6$&1.75566716318264 \\
∗ & & $b 7$& 0.0883883476483184\\
∗ \end{ tabu la r }
∗ For the func t i on $F 4 (x ) $ we use the r e cu r r ence r e l a t i o n
∗ $$
∗ F 0 (x ) = x + \ f r a c {13}{20}
∗ $$
∗ $$
∗ F { i }( x ) = x + \ f r a c {5− i }{F { i −1}(x ) } \quad \quad \ t ex t { f o r } i =1 ,2 ,3 ,4
∗ $$
∗
∗ On entry x i s a r e a l value , and we return double p r e c i s i o n es t imate o f the
cummulative normal d i s t r i b u t i o n
∗/
10
double normalDi s t r ibut ion ( double x )
{
// c a l c u l a t e \ s q r t {2\ pi } upfront once
s t a t i c const double RT2PI = sq r t ( 4 . 0∗ acos ( 0 . 0 ) ) ;
// c a l c u l a t e 10/\ s q r t {2} upfront once
s t a t i c const double SPLIT = 10 ./ sq r t (2 ) ;
s t a t i c const double a [ ] =
{220 .206867912376 ,221 .213596169931 ,112 .079291497871 ,33 .912866078383 ,6 .37396220353165 ,0 .700383064443688 ,3 .52624965998911
e−02};
s t a t i c const double b [ ] =
{440 .413735824752 ,793 .826512519948 ,637 .333633378831 ,296 .564248779674 ,86 .7807322029461 ,16 .064177579207 ,1 .75566716318264 ,8 .83883476483184
e−02};
const double z = fabs (x ) ;
// Now N(x ) = 1 − N(−x ) = 1−\ s q r t {2\ pi }N ’( x ) \ f r a c {P(x ) }{Q(x) }
// so N(−x ) = \ s q r t {2\ pi }N ’( x ) \ f r a c {P(x ) }{Q(x) }
// now l e t \ s q r t {2\ pi }N ’( z ) \ f r a c {P(x ) }{Q( z ) } = Nz
// There fore we have
// Nxm = N(x ) = \ s q r t {2\ pi }N ’( z ) \ f r a c {P(x ) }{Q( z ) } = Nz i f x<0
// Nxp = N(x ) = 1 − \ s q r t {2\ pi }N ’( z ) \ f r a c {P(x ) }{Q( z ) } = 1−Nz i f x>=0
double Nz = 0 . 0 ;
// i f z ou t s id e these l im i t s then value e f f e c t i v e l y 0 or 1 f o r machine p r e c i s i o n
i f ( z<=37.0)
{
// NDash = N ’ ( z ) ∗ s q r t {2\ pi }
const double NDash = exp(−z∗z /2 . 0 ) /RT2PI ;
i f ( z{
// here Pz = P( z ) i s a polynomial
const double Pz = ( ( ( ( ( a [ 6 ] ∗ z + a [ 5 ] ) ∗z + a [ 4 ] ) ∗z + a [ 3 ] ) ∗z + a [ 2 ] ) ∗z + a [ 1 ] )
∗z + a [ 0 ] ;
// and Qz = Q( z ) i s a polynomial
const double Qz = ( ( ( ( ( ( b [ 7 ] ∗ z + b [ 6 ] ) ∗z + b [ 5 ] ) ∗z + b [ 4 ] ) ∗z + b [ 3 ] ) ∗z + b
[ 2 ] ) ∗z + b [ 1 ] ) ∗z + b [ 0 ] ;
// use polynomia l s to c a l c u l a t e N( z ) = \ s q r t {2\ pi }N ’( x ) \ f r a c {P(x ) }{Q(x) }
Nz = RT2PI∗NDash∗Pz/Qz ;
}
e l s e
{
// implement r e cu r r ence r e l a t i o n on F 4 ( z )
const double F4z = z + 1 .0/ ( z + 2 . 0/ ( z + 3 . 0/ ( z + 4 . 0/ ( z + 13 . 0/20 . 0 ) ) ) ) ;
// use polynomia l s to c a l c u l a t e N( z ) , note here that Nz = N’ / F
Nz = NDash/F4z ;
}
}
//
return x>=0.0 ? 1−Nz : Nz ;
}
i n t main ( )
{
cout . p r e c i s i o n (16) ;
cout << ”N(0)=” << normalDi s t r ibut ion ( 0 . ) << endl ;
cout << ”N(1)=” << normalDi s t r ibut ion ( 1 . ) << endl ;
11
cout << ”N(2)=” << normalDi s t r ibut ion ( 2 . ) << endl ;
r e turn 0 ;
}
download
1.1.(iv) Using a inbuilt functions
Use the built in cmath function erfc to create an approximation to the normal distribution. Use
the relationship
N(x) =
1
2
erfc
(
− x√
2
)
.
First, start with a new project and an empty cpp file.
#inc lude
#inc lude
#inc lude
us ing namespace std ;
i n t main ( )
{
// do nothing
return 0 ;
}
download
Compile and check it runs. Now we can use the erfc function to calculate N(x). Output your
calculation with x = 1.
#inc lude
#inc lude
#inc lude
us ing namespace std ;
i n t main ( )
{
double x=1. ;
cout << 0 .5∗ e r f c (−x/ sq r t ( 2 . ) ) << endl ;
r e turn 0 ;
}
download
The output should be 0.841345. Put this into a function for ease of use.
#inc lude
#inc lude
#inc lude
us ing namespace std ;
12
double no rma lD i s t r i bu t i on bu i l t i n ( double x )
{
re turn 0 .5∗ e r f c (−x/ sq r t ( 2 . ) ) ;
}
i n t main ( )
{
f o r ( i n t x=−3;x<=3;x++)
cout << ”N( ”<}
download
1.1.(v) Comparing efficiency
Create a table of values to compare each of the approximations (i)-(iv). Which of the calculations
is most accurate? Create a loop calling each function 1 million times, which of the function is the
most efficient?
Here we use an example code as the start point. This code time how long N = 107 sin func-
tion evaluations take (plus the one or two additions/checks to run the loop). This can give you
a measure of how long functions take relative to each other. For functions that take longer to
evaluate you will need to lower the value of N .
#inc lude
#inc lude
#inc lude
#inc lude
us ing namespace std ;
i n t main ( )
{
i n t N=10000000;
// get s t a r t time
auto s t a r t = std : : chrono : : s t e ady c l o ck : : now( ) ;
// code in here i s timed
double sum=0. ;
f o r ( i n t i =0; i{
sum = sum + s in ( i ) ;
}
cout << sum << endl ;
// get f i n i s h time
auto f i n i s h = std : : chrono : : s t e ady c l o ck : : now( ) ;
// convert i n to r e a l time in seconds
auto e lapsed = std : : chrono : : dura t i on ca s t >( f i n i s h
− s t a r t ) ;
// output va lue s
cout << ” Total time e lapsed f o r ” ;
cout << N << ” c a l c u l a t i o n s i s ” << e lapsed . count ( ) << endl ;
}
13
download
Try adding in the normalDistribution codes from above in place of the the sin function and
see how long they each take – which one is fastest? Make sure you change the
build environment to Release when checking code ef-
ficiency. The code here implements this for the polynomial approximation, note that we
evaluate the function over the interval [−5, 5] to make it a fair test.
#inc lude
#inc lude
#inc lude
#inc lude
us ing namespace std ;
double no rma lD i s t r ibut i on po ly ( double x )
{
double on e ov e r s q r t 2p i = 1 . / sq r t ( 8 .∗ atan ( 1 . ) ) ;
double gamma = 0 .33267 ;
double a1 = 0 .43618 ;
double a2 = −0.12017;
double a3 = 0 .93730 ;
i f ( x >= 0)
{
double k = 1 / ( 1 . + gamma∗x ) ;
r e turn 1 . − on e ov e r s q r t 2p i ∗exp(−x∗x ∗0 . 5 ) ∗( a1∗k + a2∗k∗k + a3∗k∗k∗k ) ;
}
e l s e
{
double k = 1 / ( 1 . − gamma∗x ) ;
r e turn one ov e r s q r t 2p i ∗exp(−x∗x ∗0 . 5 ) ∗( a1∗k + a2∗k∗k + a3∗k∗k∗k ) ;
}
}
i n t main ( )
{
i n t N=1000000;
// get s t a r t time
auto s t a r t = std : : chrono : : s t e ady c l o ck : : now( ) ;
// code in here i s timed
double sum=0. ;
f o r ( i n t i =0; i{
f o r ( i n t x=−5;x<=5;x++)
sum = sum + norma lD i s t r ibut i on po ly (x ) ;
}
cout << sum << endl ;
// get f i n i s h time
auto f i n i s h = std : : chrono : : s t e ady c l o ck : : now( ) ;
// convert i n to r e a l time in seconds
auto e lapsed = std : : chrono : : dura t i on ca s t >( f i n i s h
− s t a r t ) ;
// output va lue s
cout << ” Total time e lapsed f o r ” ;
cout << 11 .∗N << ” c a l c u l a t i o n s i s ” << e lapsed . count ( ) << endl ;
14
}download
1.2 Black Scholes
Write a C++ code to calculate the option values for a call C(t = 0, S), and a put P (t = 0, S)
using the Black-Scholes formula:
C(t = 0, S) = SN(d1)−Xe−r(T−t)N(d2), (2)
P (t = 0, S) = Xe−r(T−t)N(−d2)− SN(−d1) (3)
with
d1 =
log(S/X) + (r + σ2/2)(T − t)
σ
√
T − t (4)
d2 =
log(S/X) + (r − σ2/2)(T − t)
σ
√
T − t (5)
You can use the cumulative distribution function N(x) from earlier on. So rather than starting
with an empty project copy and paste the code from the previous project into your new project,
you can find a version of the code on my website (click here). Now we need to create a function
for the call option, so we need to think about what arguments need to be supplied and what are
the local variables to the function. Obviously we need to supply the asset price, current time
both of which may vary, and the parameters for the function are the strike price, interest rate,
volatility and maturity, and we want to return a real number as the answer. Inside the function,
we will need to calculate the value of d1 and d2. This leads to a function definition of the form
double ca l lOpt i onPr i c e ( double S , double t , double X, double r , double sigma , double T)
{
double d1 ;
double d2 ;
re turn 0 ;
}
download
You must place this function header in between the main function and the normalDistribution
function, since we will be calling ”normalDistribution” inside ”callOptionPrice” and then ”cal-
lOptionPrice” inside ”main”. Next we can simply fill in the mathematical formulas as follows
double ca l lOpt i onPr i c e ( double S , double t , double X, double r , double sigma , double T)
{
double d1=( log (S/X) + ( r+sigma∗ sigma /2 . ) ∗(T−t ) ) /( sigma∗ s q r t (T−t ) ) ;
double d2=( log (S/X) + ( r−sigma∗ sigma /2 . ) ∗(T−t ) ) /( sigma∗ s q r t (T−t ) ) ;
r e turn normalDi s t r ibut ion ( d1 ) ∗S − normalDi s t r ibut ion ( d2 ) ∗X∗exp(−r ∗(T−t ) ) ;
15
}i n t main ( )
{
cout << ”Cal l Option Pr i ce = ” << ca l lOpt i onPr i c e ( 1 , 0 , 1 , 0 . 0 5 , 0 . 2 , 1 ) << endl ;
r e turn 0 ;
}
download
Now we need to test the function under different settings. There are obviously going to be
problem in each of the following cases:
• S = 0
• σ = 0
• t = T
since all will result in undefined mathematical values given the way that d1 and d2 are calculated.
However, returning a value of plus or minus infinity for d1 and d2 does not result in an undefined
value for the cumulative normal distribution since the function returns finite values from infinite
limits. Using our knowledge of what happens to this function we can go through each case and
decide what are the appropriate values to return.
Case S=0
In this situation we must make use of the boundary condition of the problem, which is that
the option value is worthless if the asset value is zero. The only slight caveat here is that you
should check whether S is smaller than a very small number relative to the strike price rather than
comparing it with zero. In code you add the following line to your function before the calculations.
i f (S<1.e−14∗X) return 0 . ;
download
By executing a return the function will stop and return the value without executing any more
code.
Case sigma=0
This case is very similar to when t=T. Depending on the sign of the numerator in the calculation
for d1 and d2 the function will either return zero or the asset minus the discounted strike price.
In code this looks like
i f ( sigma∗ sigma∗ f abs (T−t )<1.e−14)
{
i f (Se l s e re turn S−X∗exp(−r ∗(T−t ) ) ;
}
download
16
Case t=T
Finally if t and T are almost equal then we are at maturity and we return the payoff. The code
might look like
i f ( f abs (T−t )<1.e−14)
{
i f (Se l s e re turn S−X;
}
download
On adding each of these parts to the code you should test and validate each part using a range
of parameters. Note here that another case that could cause problems is if t > T . There is no
sensible value that can be returned in this case so if you add in a check for it you would be looking
to exit the program if this happened or at least print a warning to the screen. Adding this in is
left as an exercise. Your final code should look like
#inc lude
#inc lude
us ing namespace std ;
double normalDi s t r ibut ion ( double x )
{
s t a t i c const double RT2PI = sq r t ( 4 . 0∗ acos ( 0 . 0 ) ) ;
s t a t i c const double SPLIT = 10 ./ sq r t (2 ) ;
s t a t i c const double a [ ] =
{220 .206867912376 ,221 .213596169931 ,112 .079291497871 ,33 .912866078383 ,6 .37396220353165 ,0 .700383064443688 ,3 .52624965998911
e−02};
s t a t i c const double b [ ] =
{440 .413735824752 ,793 .826512519948 ,637 .333633378831 ,296 .564248779674 ,86 .7807322029461 ,16 .064177579207 ,1 .75566716318264 ,8 .83883476483184
e−02};
const double z = fabs (x ) ;
double Nz = 0 . 0 ;
// i f z ou t s id e these l im i t s then value e f f e c t i v e l y 0 or 1 f o r machine p r e c i s i o n
i f ( z<=37.0)
{
// NDash = N ’ ( z ) ∗ s q r t {2\ pi }
const double NDash = exp(−z∗z /2 . 0 ) /RT2PI ;
i f ( z{
const double Pz = ( ( ( ( ( a [ 6 ] ∗ z + a [ 5 ] ) ∗z + a [ 4 ] ) ∗z + a [ 3 ] ) ∗z + a [ 2 ] ) ∗z + a [ 1 ] )
∗z + a [ 0 ] ;
const double Qz = ( ( ( ( ( ( b [ 7 ] ∗ z + b [ 6 ] ) ∗z + b [ 5 ] ) ∗z + b [ 4 ] ) ∗z + b [ 3 ] ) ∗z + b
[ 2 ] ) ∗z + b [ 1 ] ) ∗z + b [ 0 ] ;
Nz = RT2PI∗NDash∗Pz/Qz ;
}
e l s e
{
const double F4z = z + 1 .0/ ( z + 2 . 0/ ( z + 3 . 0/ ( z + 4 . 0/ ( z + 13 . 0/20 . 0 ) ) ) ) ;
Nz = NDash/F4z ;
17
}
}
re turn x>=0.0 ? 1−Nz : Nz ;
}
// return the value o f a c a l l opt ion us ing the black s cho l e s formula
double ca l lOpt i onPr i c e ( double S , double t , double X, double r , double sigma , double T)
{
i f ( f abs (T−t )<1.e−14) // check i f we are at maturity
{
i f (Se l s e re turn S−X;
}
i f ( (T−t )<=−1.e−14) re turn 0 . ; // opt ion exp i red
i f (X<1.e−14∗S) re turn S−X∗exp(−r ∗(T−t ) ) ; // check i f s t r i k e << a s s e t then e x e r c i s e
with c e r t a i n t y
i f (S<1.e−14∗X) return 0 . ; // check i f a s s e t << s t r i k e then wor th l e s s
i f ( sigma∗ sigma ∗(T−t )<1.e−14) // check i f var i ance very smal l then no d i f f u s i o n
{
i f (Se l s e re turn S−X∗exp(−r ∗(T−t ) ) ;
}
// c a l c u l a t e opt ion p r i c e
double d1=( log (S/X) + ( r+sigma∗ sigma /2 . ) ∗(T−t ) ) /( sigma∗ s q r t (T−t ) ) ;
double d2=( log (S/X) + ( r−sigma∗ sigma /2 . ) ∗(T−t ) ) /( sigma∗ s q r t (T−t ) ) ;
r e turn normalDi s t r ibut ion ( d1 ) ∗S − normalDi s t r ibut ion ( d2 ) ∗X∗exp(−r ∗(T−t ) ) ;
}
i n t main ( )
{
cout << ”Cal l Option Pr i ce = ” << ca l lOpt i onPr i c e ( 1 , 0 , 1 , 0 . 0 5 , 0 . 2 , 1 ) << endl ;
r e turn 0 ;
}
download
Similarly for a put option we would have
#inc lude
#inc lude
us ing namespace std ;
double normalDi s t r ibut ion ( double x )
{
s t a t i c const double RT2PI = sq r t ( 4 . 0∗ acos ( 0 . 0 ) ) ;
s t a t i c const double SPLIT = 10 ./ sq r t (2 ) ;
s t a t i c const double a [ ] =
{220 .206867912376 ,221 .213596169931 ,112 .079291497871 ,33 .912866078383 ,6 .37396220353165 ,0 .700383064443688 ,3 .52624965998911
e−02};
s t a t i c const double b [ ] =
{440 .413735824752 ,793 .826512519948 ,637 .333633378831 ,296 .564248779674 ,86 .7807322029461 ,16 .064177579207 ,1 .75566716318264 ,8 .83883476483184
e−02};
const double z = fabs (x ) ;
double Nz = 0 . 0 ;
// i f z ou t s id e these l im i t s then value e f f e c t i v e l y 0 or 1 f o r machine p r e c i s i o n
18
i f ( z<=37.0)
{
// NDash = N ’ ( z ) ∗ s q r t {2\ pi }
const double NDash = exp(−z∗z /2 . 0 ) /RT2PI ;
i f ( z{
const double Pz = ( ( ( ( ( a [ 6 ] ∗ z + a [ 5 ] ) ∗z + a [ 4 ] ) ∗z + a [ 3 ] ) ∗z + a [ 2 ] ) ∗z + a [ 1 ] )
∗z + a [ 0 ] ;
const double Qz = ( ( ( ( ( ( b [ 7 ] ∗ z + b [ 6 ] ) ∗z + b [ 5 ] ) ∗z + b [ 4 ] ) ∗z + b [ 3 ] ) ∗z + b
[ 2 ] ) ∗z + b [ 1 ] ) ∗z + b [ 0 ] ;
Nz = RT2PI∗NDash∗Pz/Qz ;
}
e l s e
{
const double F4z = z + 1 .0/ ( z + 2 . 0/ ( z + 3 . 0/ ( z + 4 . 0/ ( z + 13 . 0/20 . 0 ) ) ) ) ;
Nz = NDash/F4z ;
}
}
re turn x>=0.0 ? 1−Nz : Nz ;
}
// return the value o f a put opt ion us ing the black s cho l e s formula
double putOptionPrice ( double S , double t , double X, double r , double sigma , double T)
{
i f ( f abs (T−t )<1.e−14) // check i f we are at maturity
{
i f (Se l s e re turn 0 ;
}
i f ( (T−t )<=−1.e−14) re turn 0 . ; // opt ion exp i red
i f (X<1.e−14∗S) re turn 0 . ; // check i f s t r i k e << a s s e t then e x e r c i s e with c e r t a i n t y
i f (S<1.e−14∗X) return X∗exp(−r ∗(T−t ) ) − S ; // check i f a s s e t << s t r i k e then
wor th l e s s
i f ( sigma∗ sigma ∗(T−t )<1.e−14) // check i f var i ance very smal l then no d i f f u s i o n
{
i f (Se l s e re turn 0 . ;
}
// c a l c u l a t e opt ion p r i c e
double d1=( log (S/X) + ( r+sigma∗ sigma /2 . ) ∗(T−t ) ) /( sigma∗ s q r t (T−t ) ) ;
double d2=( log (S/X) + ( r−sigma∗ sigma /2 . ) ∗(T−t ) ) /( sigma∗ s q r t (T−t ) ) ;
r e turn normalDi s t r ibut ion (−d2 ) ∗X∗exp(−r ∗(T−t ) ) − normalDi s t r ibut ion (−d1 ) ∗S ;
}
i n t main ( )
{
cout << ”Put Option Pr i ce = ” << putOptionPrice ( 1 , 0 , 1 , 0 . 0 5 , 0 . 2 , 1 ) << endl ;
r e turn 0 ;
}
download
Other Options
For a binary call we get
19
#inc lude
#inc lude
us ing namespace std ;
double normalDi s t r ibut ion ( double x )
{
s t a t i c const double RT2PI = sq r t ( 4 . 0∗ acos ( 0 . 0 ) ) ;
s t a t i c const double SPLIT = 10 ./ sq r t (2 ) ;
s t a t i c const double a [ ] =
{220 .206867912376 ,221 .213596169931 ,112 .079291497871 ,33 .912866078383 ,6 .37396220353165 ,0 .700383064443688 ,3 .52624965998911
e−02};
s t a t i c const double b [ ] =
{440 .413735824752 ,793 .826512519948 ,637 .333633378831 ,296 .564248779674 ,86 .7807322029461 ,16 .064177579207 ,1 .75566716318264 ,8 .83883476483184
e−02};
const double z = fabs (x ) ;
double Nz = 0 . 0 ;
// i f z ou t s id e these l im i t s then value e f f e c t i v e l y 0 or 1 f o r machine p r e c i s i o n
i f ( z<=37.0)
{
// NDash = N ’ ( z ) ∗ s q r t {2\ pi }
const double NDash = exp(−z∗z /2 . 0 ) /RT2PI ;
i f ( z{
const double Pz = ( ( ( ( ( a [ 6 ] ∗ z + a [ 5 ] ) ∗z + a [ 4 ] ) ∗z + a [ 3 ] ) ∗z + a [ 2 ] ) ∗z + a [ 1 ] )
∗z + a [ 0 ] ;
const double Qz = ( ( ( ( ( ( b [ 7 ] ∗ z + b [ 6 ] ) ∗z + b [ 5 ] ) ∗z + b [ 4 ] ) ∗z + b [ 3 ] ) ∗z + b
[ 2 ] ) ∗z + b [ 1 ] ) ∗z + b [ 0 ] ;
Nz = RT2PI∗NDash∗Pz/Qz ;
}
e l s e
{
const double F4z = z + 1 .0/ ( z + 2 . 0/ ( z + 3 . 0/ ( z + 4 . 0/ ( z + 13 . 0/20 . 0 ) ) ) ) ;
Nz = NDash/F4z ;
}
}
re turn x>=0.0 ? 1−Nz : Nz ;
}
// return the value o f a put opt ion us ing the black s cho l e s formula
double b inaryCal lOpt ionPr ice ( double S , double t , double X, double r , double sigma ,
double T)
{
i f ( f abs (T−t )<1.e−14) // check i f we are at maturity
{
i f (S>X) return 1 ;
e l s e re turn 0 ;
}
i f ( (T−t )<=−1.e−14) re turn 0 . ; // opt ion exp i red
i f (X<1.e−14∗S) re turn exp(−r ∗(T−t ) ) ; // check i f s t r i k e << a s s e t then e x e r c i s e
with c e r t a i n t y
i f (S<1.e−14∗X) return 0 . ; // check i f a s s e t << s t r i k e then wor th l e s s
i f ( sigma∗ sigma ∗(T−t )<1.e−14) // check i f var i ance very smal l then no d i f f u s i o n
{
i f (S>X∗exp(−r ∗(T−t ) ) ) re turn exp(−r ∗(T−t ) ) ;
20
e l s e re turn 0 . ;
}
// c a l c u l a t e opt ion p r i c e
double d2=( log (S/X) + ( r−sigma∗ sigma /2 . ) ∗(T−t ) ) /( sigma∗ s q r t (T−t ) ) ;
r e turn normalDi s t r ibut ion ( d2 ) ∗exp(−r ∗(T−t ) ) ;
}
i n t main ( )
{
cout << ”Binary Ca l l Option Pr i ce = ” << binaryCal lOpt ionPr ice ( 1 , 0 , 1 , 0 . 0 5 , 0 . 2 , 1 )
<< endl ;
r e turn 0 ;
}
download
and for a binary put
#inc lude
#inc lude
us ing namespace std ;
double normalDi s t r ibut ion ( double x )
{
s t a t i c const double RT2PI = sq r t ( 4 . 0∗ acos ( 0 . 0 ) ) ;
s t a t i c const double SPLIT = 10 ./ sq r t (2 ) ;
s t a t i c const double a [ ] =
{220 .206867912376 ,221 .213596169931 ,112 .079291497871 ,33 .912866078383 ,6 .37396220353165 ,0 .700383064443688 ,3 .52624965998911
e−02};
s t a t i c const double b [ ] =
{440 .413735824752 ,793 .826512519948 ,637 .333633378831 ,296 .564248779674 ,86 .7807322029461 ,16 .064177579207 ,1 .75566716318264 ,8 .83883476483184
e−02};
const double z = fabs (x ) ;
double Nz = 0 . 0 ;
// i f z ou t s id e these l im i t s then value e f f e c t i v e l y 0 or 1 f o r machine p r e c i s i o n
i f ( z<=37.0)
{
// NDash = N ’ ( z ) ∗ s q r t {2\ pi }
const double NDash = exp(−z∗z /2 . 0 ) /RT2PI ;
i f ( z{
const double Pz = ( ( ( ( ( a [ 6 ] ∗ z + a [ 5 ] ) ∗z + a [ 4 ] ) ∗z + a [ 3 ] ) ∗z + a [ 2 ] ) ∗z + a [ 1 ] )
∗z + a [ 0 ] ;
const double Qz = ( ( ( ( ( ( b [ 7 ] ∗ z + b [ 6 ] ) ∗z + b [ 5 ] ) ∗z + b [ 4 ] ) ∗z + b [ 3 ] ) ∗z + b
[ 2 ] ) ∗z + b [ 1 ] ) ∗z + b [ 0 ] ;
Nz = RT2PI∗NDash∗Pz/Qz ;
}
e l s e
{
const double F4z = z + 1 .0/ ( z + 2 . 0/ ( z + 3 . 0/ ( z + 4 . 0/ ( z + 13 . 0/20 . 0 ) ) ) ) ;
Nz = NDash/F4z ;
}
}
re turn x>=0.0 ? 1−Nz : Nz ;
}
21
// return the value o f a put opt ion us ing the black s cho l e s formula
double binaryPutOptionPrice ( double S , double t , double X, double r , double sigma , double
T)
{
i f ( f abs (T−t )<1.e−14) // check i f we are at maturity
{
i f (Se l s e re turn 0 ;
}
i f ( (T−t )<=−1.e−14) re turn 0 . ; // opt ion exp i red
i f (X<1.e−14∗S) re turn 0 . ; // check i f s t r i k e << a s s e t then e x e r c i s e with c e r t a i n t y
i f (S<1.e−14∗X) return exp(−r ∗(T−t ) ) ; // check i f a s s e t << s t r i k e then wor th l e s s
i f ( sigma∗ sigma ∗(T−t )<1.e−14) // check i f var i ance very smal l then no d i f f u s i o n
{
i f (Se l s e re turn 0 . ;
}
// c a l c u l a t e opt ion p r i c e
double d2=( log (S/X) + ( r−sigma∗ sigma /2 . ) ∗(T−t ) ) /( sigma∗ s q r t (T−t ) ) ;
r e turn normalDi s t r ibut ion (−d2 ) ∗exp(−r ∗(T−t ) ) ;
}
i n t main ( )
{
cout << ”Binary Put Option Pr i ce = ” << binaryPutOptionPrice ( 1 , 0 , 1 , 0 . 0 5 , 0 . 2 , 1 ) <<
endl ;
r e turn 0 ;
}
download
1.3 Coursework Example
A trader has asked you to price the value of the call option C(S, t) at time t = 0 according to the
standard Black-Scholes formula, where T = 1, X = 1, r = 0.05 and σ = 0.2. Write a program to
calculate C and output the results to screen. You must generate four columns of data:
• the value of S;
• the value of d1;
• the value of d1;
• the value of C(S, t = 0).
Output each of the values when the stock price is
S ∈ {0.8, 0.9, 1, 1.1, 1.2}.
You should use a for loop to generate the data.
22
We could do this with the call option code above, but to make outputting d1 and d2 easier lets
just calculate them again. First, open up a new project and copy in the code from earlier on with
normal distribution, get it by (clicking here). Then simply start by editing the main function so
that has variables for all the required inputs (S, t, T , X, r and σ), and then calculate d1 and N(d1).
i n t main ( )
{
double S=1. , t =0. ,X=1. , r =0.05 , sigma=0.2 ,T=1. ;
double d 1=( log (S/X) + ( r+sigma∗ sigma /2 . ) ∗(T−t ) ) / sigma/ sq r t (T−t ) ;
cout << ”d 1 = ” << d 1 << endl ;
cout << ”N( d 1 ) = ” << normalDi s t r ibut ion ( d 1 ) << endl ;
r e turn 0 ;
}
download
Once you have checked this result, do the same for d2 and N(d2) then calculate the call option
price
i n t main ( )
{
double S=1. , t =0. ,X=1. , r =0.05 , sigma=0.2 ,T=1. ;
double d 1=( log (S/X) + ( r+sigma∗ sigma /2 . ) ∗(T−t ) ) / sigma/ sq r t (T−t ) ;
double d 2=( log (S/X) + ( r−sigma∗ sigma /2 . ) ∗(T−t ) ) / sigma/ sq r t (T−t ) ;
cout << ”d 1 = ” << d 1 << endl ;
cout << ”N( d 1 ) = ” << normalDi s t r ibut ion ( d 1 ) << endl ;
cout << ”d 2 = ” << d 2 << endl ;
cout << ”N( d 2 ) = ” << normalDi s t r ibut ion ( d 2 ) << endl ;
cout << ”C = ” << S∗ normalDi s t r ibut ion ( d 1 ) − X∗exp(−r ∗(T−t ) ) ∗ normalDi s t r ibut ion (
d 2 ) << endl ;
r e turn 0 ;
}
download
If you are happy everything is working then run a loop to calculate values for different values
of S, something like this should work
i n t main ( )
{
double S=1. , t =0. ,X=1. , r =0.05 , sigma=0.2 ,T=1. ;
f o r ( i n t i =8; i <=12; i++)
{
S = i ∗ 0 . 1 ;
cout << S ;
double d 1=( log (S/X) + ( r+sigma∗ sigma /2 . ) ∗(T−t ) ) / sigma/ sq r t (T−t ) ;
double d 2=( log (S/X) + ( r−sigma∗ sigma /2 . ) ∗(T−t ) ) / sigma/ sq r t (T−t ) ;
cout << ” ” << d 1 ;
cout << ” ” << d 2 ;
23
cout << ” ” << S∗ normalDi s t r ibut ion ( d 1 ) − X∗exp(−r ∗(T−t ) ) ∗
normalDi s t r ibut ion ( d 2 ) << endl ;
}
re turn 0 ;
}
download
Now you have your solution, export this manually (from terminal output) or using a filestream
in c++ so that you can create a table. Your final result should look something like this:-
References
West, Graeme. Better approximations to cumulative normal functions. Wilmott Magazine, 9:
70–76, 2005. URL https://lyle.smu.edu/~aleskovs/emis/sqc2/accuratecumnorm.pdf.
24




























#inc lude
#inc lude
us ing namespace std ;
1
i n t main ( )
{
// do nothing
return 0 ;
}
download
Compile and check it runs. Looking at the algorithm in the question we see that we need several
variables, k, γ, a1, a2 and a3. Declare all those variables along with a default value of x = 1.
#inc lude
#inc lude
#inc lude
us ing namespace std ;
i n t main ( )
{
double x=1;
double gamma=0.33267 , a1=0.43618 , a2=−0.12017 , a3=0.93730;
double k=1./(1.+gamma∗x ) ;
// do nothing
return 0 ;
}
download
Compile and check it runs. Now enter the calculation for N(x) with x ≥ 0 and check it works.
#inc lude
#inc lude
#inc lude
us ing namespace std ;
i n t main ( )
{
double x=1;
double gamma=0.33267 , a1=0.43618 , a2=−0.12017 , a3=0.93730;
double k=1./(1.+gamma∗x ) ;
cout << 1.−1./ sq r t ( 8 .∗ atan ( 1 . ) ) ∗exp(−x∗x /2 . ) ∗( a1∗k+a2∗k∗k+a3∗k∗k∗k ) << endl ;
// do nothing
return 0 ;
}
download
Compile and check the output, it should be 0.841352. Now you need to make it work when x < 0,
this means using an textttif statement and changing the value of γ as well as the formula. Check
your code works in both cases.
#inc lude
#inc lude
#inc lude
us ing namespace std ;
2
i n t main ( )
{
double x=−1;
double gamma=0.33267 , a1=0.43618 , a2=−0.12017 , a3=0.93730;
i f (x>=0)
{
double k=1./(1.+gamma∗x ) ;
cout << 1.−1./ sq r t ( 8 .∗ atan ( 1 . ) ) ∗exp(−x∗x /2 . ) ∗( a1∗k+a2∗k∗k+a3∗k∗k∗k ) << endl ;
}
e l s e
{
double k=1./(1.−gamma∗x ) ;
cout << 1 ./ sq r t ( 8 .∗ atan ( 1 . ) ) ∗exp(−x∗x /2 . ) ∗( a1∗k+a2∗k∗k+a3∗k∗k∗k ) << endl ;
}
// do nothing
return 0 ;
}
download
Finally, move your code into a function and test the result. You will need to use the return
keyword to return the value of the function instead of printing to screen.
#inc lude
#inc lude
#inc lude
us ing namespace std ;
double no rma lD i s t r ibut i on po ly ( double x )
{
s t a t i c double one ove r sq r t 2 =1./ sq r t ( 8 .∗ atan ( 1 . ) ) ;
double gamma=0.33267 , a1=0.43618 , a2=−0.12017 , a3=0.93730;
i f (x>=0)
{
double k=1./(1.+gamma∗x ) ;
r e turn 1.− one ove r sq r t 2 ∗exp(−x∗x /2 . ) ∗( a1∗k+a2∗k∗k+a3∗k∗k∗k ) ;
}
e l s e
{
double k=1./(1.−gamma∗x ) ;
r e turn one ove r sq r t 2 ∗exp(−x∗x /2 . ) ∗( a1∗k+a2∗k∗k+a3∗k∗k∗k ) ;
}
}
i n t main ( )
{
f o r ( i n t x=−3;x<=3;x++)
cout << ”N( ”<// do nothing
return 0 ;
}
download
The output from this code is:
N(-3) = 0.00135489
3
N(-2) = 0.0227587
N(-1) = 0.158648
N(0) = 0.500002
N(1) = 0.841352
N(2) = 0.977241
N(3) = 0.998645
1.1.(ii) Using integration
Use Simpson’s rule to evaluate the integral defining N(x). This states that when there are values
of y = f(x) for values of x at equal intervals h apart, an approximate numerical integration is
given by ∫ b
x=a
f(x)dx ≈ h
3
[y0 + 4y1 + 2y2 + 4y3 + 2y4 + . . .+ 4yn−1 + yn]
where h = x1 − x0 = (b − a)/n. Note this formula involves an even number of intervals, and
the (truncation) error diminishes as h4. You will need around 2000 points or more for machine
accuracy, making this a very expensive method.
Again, first we think about the stages of our program. We must choose a method to implement
the integration, test and verify it. Once it is then applied to the problem we need to decide how
to deal with the concept of infinity in this setting. First let use code up the integration method.
To be able to check the method, you should choose a simple function with a known integral so
that we can compare accuracy. Start with a new project and an empty cpp file.
#inc lude
#inc lude
us ing namespace std ;
i n t main ( )
{
// do nothing
return 0 ;
}
download
Compile and check it runs. We are going to implement Simpson’s method for integration firstly
on the sin function so that we can check the results. The method is given by
∫ b
a
f(x)dx ≈ h
3
[
f(a) + 2
n/2−1∑
j=1
f(a+ 2jh)
+4
n/2∑
j=1
f(a+ (2j − 1)h) + f(b)
] (1)
Start by declaring the input parameters required by the method, including the number of
steps and the integration range, and follow with the variables used inside such as the dummy
4
integration variable s and the step size h. We should supply default values and initialise variables
if required
// number o f s t ep s
i n t N=10;
// range o f i n t e g r a t i o n
double a=0,b=1. ;
// l o c a l v a r i a b l e s
double s , h , sum=0. ;
// i n i a l i s e the v a r i a b l e s
h=(b−a ) /N;
download
Next we need to add the first few terms and the last one to sum, before adding a loop to go over
the rest.
// add in the f i r s t few terms
sum = sum + s in ( a ) + 4 .∗ s i n ( a+h) ;
// and the l a s t one
sum = sum + s in (b) ;
download
Note that we should really check that N is even as this method relies on that fact. Since we have
already taken care of the terms for 0, 1 and N, the loop must run over the terms 2 up to and
including N-1. Check your loop runs over the correct values before continuing
// loop over terms 2 up to N−1
f o r ( i n t i =1; i{
s = a + 2∗ i ∗h ;
// check even terms
cout << 2∗ i << ” ” << s << endl ;
s = s + h ;
// check odd terms
cout << 2∗ i+1 << ” ” << s << endl ;
}
download
Now add them in and check against the analytical solution.
i n t main ( )
{
// number o f s t ep s
i n t N=10;
// range o f i n t e g r a t i o n
double a=0,b=1. ;
// l o c a l v a r i a b l e s
double s , h , sum=0. ;
// i n i a l i s e the v a r i a b l e s
h=(b−a ) /N;
// add in the f i r s t few terms
5
sum = sum + s in ( a ) + 4 .∗ s i n ( a+h) ;
// and the l a s t one
sum = sum + s in (b) ;
// loop over terms 2 up to N−1
f o r ( i n t i =1; i{
s = a + 2∗ i ∗h ;
sum = sum + 2.∗ s i n ( s ) ;
s = s + h ;
sum = sum + 4.∗ s i n ( s ) ;
}
// complete the i n t e g r a l
sum = h∗sum/ 3 . ;
// output r e s u l t s
cout << sum << ” ” << 1.− cos (b) << endl ;
r e turn 0 ;
}
download
You should check accuracy by varying the value of N.
Next we are going to change this algorithm to solve the normal distribution. Do this in 2
steps, first go through and change all calls to the sin function by exp(-s*s/2.), then multiply the
final result by 1/sqrt(2*pi). The result (for a=0 and b=1) should be 0.341345. Now we have to
deal with infinity, so how can we integrate over infinity? For the normal distribution this is quite
easy as it is an even function and we know that the integral between 0 and infinity is one half.
So simply adding one half to the result of the integration between 0 and b gives the value of N(b)
for all b. At this stage your code should look like
// add in the f i r s t few terms
sum = sum + exp(−a∗a /2 . ) + 4 .∗ exp(−(a+h) ∗( a+h) /2 . ) ;
// and the l a s t one
sum = sum + exp(−b∗b /2 . ) ;
// loop over terms 2 up to N−1
f o r ( i n t i =1; i{
s = a + 2∗ i ∗h ;
sum = sum + 2.∗ exp(−s ∗ s / 2 . ) ;
s = s + h ;
sum = sum + 4.∗ exp(−s ∗ s / 2 . ) ;
}
// complete the i n t e g r a l
sum = 0.5 + h∗sum/3 ./ sq r t ( 8 .∗ atan ( 1 . ) ) ;
// output r e s u l t s
cout << sum << endl ;
download
Now move this algorithm into a function with the definition
double normalDi s t r ibut ion ( double x )
download
6
The last thing to take account of is what happens when x is either very large and positive or very
large and negative. The best thing to do here is use analytical results and return 0 if x is large
and negative and 1 if x if large and positive. Don’t forget to set the value of b to be x and choose
a value of N sufficiently large to maintain accuracy.
#inc lude
#inc lude
us ing namespace std ;
double normalDi s t r ibut ion ( double x )
{
i f (x<−10.) re turn 0 . ;
i f (x>10.) re turn 1 . ;
// number o f s t ep s
i n t N=2000;
// range o f i n t e g r a t i o n
double a=0,b=x ;
// l o c a l v a r i a b l e s
double s , h , sum=0. ;
// i n i a l i s e the v a r i a b l e s
h=(b−a ) /N;
// add in the f i r s t few terms
sum = sum + exp(−a∗a /2 . ) + 4 .∗ exp(−(a+h) ∗( a+h) /2 . ) ;
// and the l a s t one
sum = sum + exp(−b∗b /2 . ) ;
// loop over terms 2 up to N−1
f o r ( i n t i =1; i{
s = a + 2∗ i ∗h ;
sum = sum + 2.∗ exp(−s ∗ s / 2 . ) ;
s = s + h ;
sum = sum + 4.∗ exp(−s ∗ s / 2 . ) ;
}
// complete the i n t e g r a l
sum = 0.5 + h∗sum/3 ./ sq r t ( 8 .∗ atan ( 1 . ) ) ;
// re turn r e s u l t
re turn sum ;
}
i n t main ( )
{
cout . p r e c i s i o n (16) ;
cout << ”N(0)=” << normalDi s t r ibut ion ( 0 . ) << endl ;
cout << ”N(1)=” << normalDi s t r ibut ion ( 1 . ) << endl ;
cout << ”N(2)=” << normalDi s t r ibut ion ( 2 . ) << endl ;
r e turn 0 ;
}
download
Here we needed to choose N=2000 to get close to double precision, meaning that a single
function call is incredibly expensive. There are other, more efficient approximations out there.
If you are using a compiler that is math library C99 compliant you will have access to the erfc
7
function, which will give you double precision accuracy using the formula
N(x) =
1
2
erfc
(
− x√
2pi
)
1.1.(iii) Using a higher order polynomial
Using a fast and accurate polynomial approximation to N(x) (a double precision approximation,
see West, 2005, for more details):
N(x) = 1−
√
2piN ′(x)
P (x)
Q(x)
, 0 ≤ x ≤ 10√
2
,
N(x) = 1−N ′(x) 1
F4(x)
,
10√
2
< x ≤ 37,
N(x) = 1, x > 37,
N(x) = 1−N(−x), x < 0.
Here
P (x) =
6∑
i=0
aix
i, Q(x) =
7∑
i=0
bix
i
and the coefficients are:
a0 220.206867912376 b0 440.413735824752
a1 221.213596169931 b1 793.826512519948
a2 112.079291497871 b2 637.333633378831
a3 33.912866078383 b3 296.564248779674
a4 6.37396220353165 b4 86.7807322029461
a5 0.700383064443688 b5 16.064177579207
a6 0.0352624965998911 b6 1.75566716318264
b7 0.0883883476483184
For the function F4(x) use the recurrence relation
F0(x) = x+
13
20
Fi(x) = x+
5− i
Fi−1(x)
for i = 1, 2, 3, 4
Unfortunately if you are using Visual Studio it is not yet implemented (might yet be in VS2013)
so we need another method. One such method is described in West (2005) and they give a VB
code implementation of the method. We are going to follow that method closely. Start with your
empty program, I have included some benchmark values in comments to check against
8
#inc lude
#inc lude
us ing namespace std ;
i n t main ( )
{
// Benchmark va lue s o f N(x )
// N(0) = 0 .5
// N(1) = 0.841344746068543
// N(2) = 0.977249868051821
return 0 ;
}
download
Compile and check the code runs with no output.
The algorithm is given by the following expressions
N(x) = 1−
√
2piN ′(x)
P (x)
Q(x)
, 0 ≤ x ≤ 10√
2
,
N(x) = 1−N ′(x) 1
F (x)
,
10√
2
< x ≤ 37,
N(x) = 1, x > 37,
N(x) = 1−N(−x), x < 0.
where P , Q and F are polynomials, more details can be found in the examples sheet (click here).
First declare and assign some of the variables needed for the calculation
// c a l c u l a t e \ s q r t {2\ pi } upfront once
double RT2PI = sq r t ( 4 . 0∗ acos ( 0 . 0 ) ) ;
// c a l c u l a t e 10/\ s q r t {2} upfront once
double SPLIT = 10 ./ sq r t (2 ) ;
// ente r c o e f f i c i e n t s f o r the polynomia l s P and Q
double a [ ] =
{220 .206867912376 ,221 .213596169931 ,112 .079291497871 ,33 .912866078383 ,6 .37396220353165 ,0 .700383064443688 ,3 .52624965998911
e−02};
double b [ ] =
{440 .413735824752 ,793 .826512519948 ,637 .333633378831 ,296 .564248779674 ,86 .7807322029461 ,16 .064177579207 ,1 .75566716318264 ,8 .83883476483184
e−02};
// c a l c u l a t e the value o f N at x=1
double x=1. ;
download
where we have chosen x = 1 to test. Since we know what x is we don’t need to consider the
3 cases. Now calculate the polynomial expressions using a nested multiplication (important for
accuracy) and we get
double NDash = exp(−x∗x /2 . 0 ) /RT2PI ;
double Px = ( ( ( ( ( a [ 6 ] ∗ x + a [ 5 ] ) ∗x + a [ 4 ] ) ∗x + a [ 3 ] ) ∗x + a [ 2 ] ) ∗x + a [ 1 ] ) ∗x + a [ 0 ] ;
9
double Qx = ( ( ( ( ( ( b [ 7 ] ∗ x + b [ 6 ] ) ∗x + b [ 5 ] ) ∗x + b [ 4 ] ) ∗x + b [ 3 ] ) ∗x + b [ 2 ] ) ∗x + b
[ 1 ] ) ∗x + b [ 0 ] ;
cout << Px << ” ” << Qx << ” ” << NDash << endl ;
// output i s 594 .522 2272.83 0.241971
download
the output you should expect is shown. We can now estimate the value of N(x)
cout . p r e c i s i o n (16) ;
cout << 1 − exp(−x∗x /2 . 0 ) ∗Px/Qx << endl ;
// output i s 0.841344746068543
download
You now need to include the different cases, such as negative x, etc, and then move everything
into a function. After some testing you should arrive at an algorithm that looks something like:
#inc lude
#inc lude
us ing namespace std ;
/∗
∗ A ∗ ∗ f a s t ∗ and accurate polynomial approximation to $N(x ) $ ( double p r e c i s i o n
) :
∗ $$ N(x )= 1−\ s q r t {2\ pi }N ’( x ) \ f r a c {P(x ) }{Q(x) } , \quad 0 \ l e q x\ l e q \ f r a c {10}{\
s q r t {2}} , $$
∗ $$ N(x )= 1−N ’( x ) \ f r a c {1}{F 4 (x ) } , \quad \ f r a c {10}{\ s q r t {2}} < x\ l e q 37 , $$
∗ $$ N(x )= 1 , \quad x > 37 , $$
∗ $$N(x )=1−N(−x ) ,\quad x<0.$$
∗ Here
∗ $$
∗ P(x ) = \sum { i =0}ˆ6 a i xˆ i , \quad \quad Q(x ) = \sum { i =0}ˆ7 b i xˆ i
∗ $$
∗ and the c o e f f i c i e n t s are :
∗ \begin { tabu la r }{ l c | | l c }
∗ \ h l i n e
∗ $a 0$ & 220.206867912376& $b 0$&440.413735824752 \\
∗ $a 1$ & 221.213596169931& $b 1$&793.826512519948 \\
∗ $a 2$ & 112.079291497871& $b 2$& 637.333633378831\\
∗ $a 3$ & 33.912866078383& $b 3$&296.564248779674 \\
∗ $a 4$ & 6.37396220353165& $b 4$&86.7807322029461 \\
∗ $a 5$ & 0.700383064443688& $b 5$& 16.064177579207\\
∗ $a 6$ & 0.0352624965998911& $b 6$&1.75566716318264 \\
∗ & & $b 7$& 0.0883883476483184\\
∗ \end{ tabu la r }
∗ For the func t i on $F 4 (x ) $ we use the r e cu r r ence r e l a t i o n
∗ $$
∗ F 0 (x ) = x + \ f r a c {13}{20}
∗ $$
∗ $$
∗ F { i }( x ) = x + \ f r a c {5− i }{F { i −1}(x ) } \quad \quad \ t ex t { f o r } i =1 ,2 ,3 ,4
∗ $$
∗
∗ On entry x i s a r e a l value , and we return double p r e c i s i o n es t imate o f the
cummulative normal d i s t r i b u t i o n
∗/
10
double normalDi s t r ibut ion ( double x )
{
// c a l c u l a t e \ s q r t {2\ pi } upfront once
s t a t i c const double RT2PI = sq r t ( 4 . 0∗ acos ( 0 . 0 ) ) ;
// c a l c u l a t e 10/\ s q r t {2} upfront once
s t a t i c const double SPLIT = 10 ./ sq r t (2 ) ;
s t a t i c const double a [ ] =
{220 .206867912376 ,221 .213596169931 ,112 .079291497871 ,33 .912866078383 ,6 .37396220353165 ,0 .700383064443688 ,3 .52624965998911
e−02};
s t a t i c const double b [ ] =
{440 .413735824752 ,793 .826512519948 ,637 .333633378831 ,296 .564248779674 ,86 .7807322029461 ,16 .064177579207 ,1 .75566716318264 ,8 .83883476483184
e−02};
const double z = fabs (x ) ;
// Now N(x ) = 1 − N(−x ) = 1−\ s q r t {2\ pi }N ’( x ) \ f r a c {P(x ) }{Q(x) }
// so N(−x ) = \ s q r t {2\ pi }N ’( x ) \ f r a c {P(x ) }{Q(x) }
// now l e t \ s q r t {2\ pi }N ’( z ) \ f r a c {P(x ) }{Q( z ) } = Nz
// There fore we have
// Nxm = N(x ) = \ s q r t {2\ pi }N ’( z ) \ f r a c {P(x ) }{Q( z ) } = Nz i f x<0
// Nxp = N(x ) = 1 − \ s q r t {2\ pi }N ’( z ) \ f r a c {P(x ) }{Q( z ) } = 1−Nz i f x>=0
double Nz = 0 . 0 ;
// i f z ou t s id e these l im i t s then value e f f e c t i v e l y 0 or 1 f o r machine p r e c i s i o n
i f ( z<=37.0)
{
// NDash = N ’ ( z ) ∗ s q r t {2\ pi }
const double NDash = exp(−z∗z /2 . 0 ) /RT2PI ;
i f ( z{
// here Pz = P( z ) i s a polynomial
const double Pz = ( ( ( ( ( a [ 6 ] ∗ z + a [ 5 ] ) ∗z + a [ 4 ] ) ∗z + a [ 3 ] ) ∗z + a [ 2 ] ) ∗z + a [ 1 ] )
∗z + a [ 0 ] ;
// and Qz = Q( z ) i s a polynomial
const double Qz = ( ( ( ( ( ( b [ 7 ] ∗ z + b [ 6 ] ) ∗z + b [ 5 ] ) ∗z + b [ 4 ] ) ∗z + b [ 3 ] ) ∗z + b
[ 2 ] ) ∗z + b [ 1 ] ) ∗z + b [ 0 ] ;
// use polynomia l s to c a l c u l a t e N( z ) = \ s q r t {2\ pi }N ’( x ) \ f r a c {P(x ) }{Q(x) }
Nz = RT2PI∗NDash∗Pz/Qz ;
}
e l s e
{
// implement r e cu r r ence r e l a t i o n on F 4 ( z )
const double F4z = z + 1 .0/ ( z + 2 . 0/ ( z + 3 . 0/ ( z + 4 . 0/ ( z + 13 . 0/20 . 0 ) ) ) ) ;
// use polynomia l s to c a l c u l a t e N( z ) , note here that Nz = N’ / F
Nz = NDash/F4z ;
}
}
//
return x>=0.0 ? 1−Nz : Nz ;
}
i n t main ( )
{
cout . p r e c i s i o n (16) ;
cout << ”N(0)=” << normalDi s t r ibut ion ( 0 . ) << endl ;
cout << ”N(1)=” << normalDi s t r ibut ion ( 1 . ) << endl ;
11
cout << ”N(2)=” << normalDi s t r ibut ion ( 2 . ) << endl ;
r e turn 0 ;
}
download
1.1.(iv) Using a inbuilt functions
Use the built in cmath function erfc to create an approximation to the normal distribution. Use
the relationship
N(x) =
1
2
erfc
(
− x√
2
)
.
First, start with a new project and an empty cpp file.
#inc lude
#inc lude
#inc lude
us ing namespace std ;
i n t main ( )
{
// do nothing
return 0 ;
}
download
Compile and check it runs. Now we can use the erfc function to calculate N(x). Output your
calculation with x = 1.
#inc lude
#inc lude
#inc lude
us ing namespace std ;
i n t main ( )
{
double x=1. ;
cout << 0 .5∗ e r f c (−x/ sq r t ( 2 . ) ) << endl ;
r e turn 0 ;
}
download
The output should be 0.841345. Put this into a function for ease of use.
#inc lude
#inc lude
#inc lude
us ing namespace std ;
12
double no rma lD i s t r i bu t i on bu i l t i n ( double x )
{
re turn 0 .5∗ e r f c (−x/ sq r t ( 2 . ) ) ;
}
i n t main ( )
{
f o r ( i n t x=−3;x<=3;x++)
cout << ”N( ”<}
download
1.1.(v) Comparing efficiency
Create a table of values to compare each of the approximations (i)-(iv). Which of the calculations
is most accurate? Create a loop calling each function 1 million times, which of the function is the
most efficient?
Here we use an example code as the start point. This code time how long N = 107 sin func-
tion evaluations take (plus the one or two additions/checks to run the loop). This can give you
a measure of how long functions take relative to each other. For functions that take longer to
evaluate you will need to lower the value of N .
#inc lude
#inc lude
#inc lude
#inc lude
us ing namespace std ;
i n t main ( )
{
i n t N=10000000;
// get s t a r t time
auto s t a r t = std : : chrono : : s t e ady c l o ck : : now( ) ;
// code in here i s timed
double sum=0. ;
f o r ( i n t i =0; i{
sum = sum + s in ( i ) ;
}
cout << sum << endl ;
// get f i n i s h time
auto f i n i s h = std : : chrono : : s t e ady c l o ck : : now( ) ;
// convert i n to r e a l time in seconds
auto e lapsed = std : : chrono : : dura t i on ca s t >( f i n i s h
− s t a r t ) ;
// output va lue s
cout << ” Total time e lapsed f o r ” ;
cout << N << ” c a l c u l a t i o n s i s ” << e lapsed . count ( ) << endl ;
}
13
download
Try adding in the normalDistribution codes from above in place of the the sin function and
see how long they each take – which one is fastest? Make sure you change the
build environment to Release when checking code ef-
ficiency. The code here implements this for the polynomial approximation, note that we
evaluate the function over the interval [−5, 5] to make it a fair test.
#inc lude
#inc lude
#inc lude
#inc lude
us ing namespace std ;
double no rma lD i s t r ibut i on po ly ( double x )
{
double on e ov e r s q r t 2p i = 1 . / sq r t ( 8 .∗ atan ( 1 . ) ) ;
double gamma = 0 .33267 ;
double a1 = 0 .43618 ;
double a2 = −0.12017;
double a3 = 0 .93730 ;
i f ( x >= 0)
{
double k = 1 / ( 1 . + gamma∗x ) ;
r e turn 1 . − on e ov e r s q r t 2p i ∗exp(−x∗x ∗0 . 5 ) ∗( a1∗k + a2∗k∗k + a3∗k∗k∗k ) ;
}
e l s e
{
double k = 1 / ( 1 . − gamma∗x ) ;
r e turn one ov e r s q r t 2p i ∗exp(−x∗x ∗0 . 5 ) ∗( a1∗k + a2∗k∗k + a3∗k∗k∗k ) ;
}
}
i n t main ( )
{
i n t N=1000000;
// get s t a r t time
auto s t a r t = std : : chrono : : s t e ady c l o ck : : now( ) ;
// code in here i s timed
double sum=0. ;
f o r ( i n t i =0; i{
f o r ( i n t x=−5;x<=5;x++)
sum = sum + norma lD i s t r ibut i on po ly (x ) ;
}
cout << sum << endl ;
// get f i n i s h time
auto f i n i s h = std : : chrono : : s t e ady c l o ck : : now( ) ;
// convert i n to r e a l time in seconds
auto e lapsed = std : : chrono : : dura t i on ca s t >( f i n i s h
− s t a r t ) ;
// output va lue s
cout << ” Total time e lapsed f o r ” ;
cout << 11 .∗N << ” c a l c u l a t i o n s i s ” << e lapsed . count ( ) << endl ;
14
}download
1.2 Black Scholes
Write a C++ code to calculate the option values for a call C(t = 0, S), and a put P (t = 0, S)
using the Black-Scholes formula:
C(t = 0, S) = SN(d1)−Xe−r(T−t)N(d2), (2)
P (t = 0, S) = Xe−r(T−t)N(−d2)− SN(−d1) (3)
with
d1 =
log(S/X) + (r + σ2/2)(T − t)
σ
√
T − t (4)
d2 =
log(S/X) + (r − σ2/2)(T − t)
σ
√
T − t (5)
You can use the cumulative distribution function N(x) from earlier on. So rather than starting
with an empty project copy and paste the code from the previous project into your new project,
you can find a version of the code on my website (click here). Now we need to create a function
for the call option, so we need to think about what arguments need to be supplied and what are
the local variables to the function. Obviously we need to supply the asset price, current time
both of which may vary, and the parameters for the function are the strike price, interest rate,
volatility and maturity, and we want to return a real number as the answer. Inside the function,
we will need to calculate the value of d1 and d2. This leads to a function definition of the form
double ca l lOpt i onPr i c e ( double S , double t , double X, double r , double sigma , double T)
{
double d1 ;
double d2 ;
re turn 0 ;
}
download
You must place this function header in between the main function and the normalDistribution
function, since we will be calling ”normalDistribution” inside ”callOptionPrice” and then ”cal-
lOptionPrice” inside ”main”. Next we can simply fill in the mathematical formulas as follows
double ca l lOpt i onPr i c e ( double S , double t , double X, double r , double sigma , double T)
{
double d1=( log (S/X) + ( r+sigma∗ sigma /2 . ) ∗(T−t ) ) /( sigma∗ s q r t (T−t ) ) ;
double d2=( log (S/X) + ( r−sigma∗ sigma /2 . ) ∗(T−t ) ) /( sigma∗ s q r t (T−t ) ) ;
r e turn normalDi s t r ibut ion ( d1 ) ∗S − normalDi s t r ibut ion ( d2 ) ∗X∗exp(−r ∗(T−t ) ) ;
15
}i n t main ( )
{
cout << ”Cal l Option Pr i ce = ” << ca l lOpt i onPr i c e ( 1 , 0 , 1 , 0 . 0 5 , 0 . 2 , 1 ) << endl ;
r e turn 0 ;
}
download
Now we need to test the function under different settings. There are obviously going to be
problem in each of the following cases:
• S = 0
• σ = 0
• t = T
since all will result in undefined mathematical values given the way that d1 and d2 are calculated.
However, returning a value of plus or minus infinity for d1 and d2 does not result in an undefined
value for the cumulative normal distribution since the function returns finite values from infinite
limits. Using our knowledge of what happens to this function we can go through each case and
decide what are the appropriate values to return.
Case S=0
In this situation we must make use of the boundary condition of the problem, which is that
the option value is worthless if the asset value is zero. The only slight caveat here is that you
should check whether S is smaller than a very small number relative to the strike price rather than
comparing it with zero. In code you add the following line to your function before the calculations.
i f (S<1.e−14∗X) return 0 . ;
download
By executing a return the function will stop and return the value without executing any more
code.
Case sigma=0
This case is very similar to when t=T. Depending on the sign of the numerator in the calculation
for d1 and d2 the function will either return zero or the asset minus the discounted strike price.
In code this looks like
i f ( sigma∗ sigma∗ f abs (T−t )<1.e−14)
{
i f (Se l s e re turn S−X∗exp(−r ∗(T−t ) ) ;
}
download
16
Case t=T
Finally if t and T are almost equal then we are at maturity and we return the payoff. The code
might look like
i f ( f abs (T−t )<1.e−14)
{
i f (Se l s e re turn S−X;
}
download
On adding each of these parts to the code you should test and validate each part using a range
of parameters. Note here that another case that could cause problems is if t > T . There is no
sensible value that can be returned in this case so if you add in a check for it you would be looking
to exit the program if this happened or at least print a warning to the screen. Adding this in is
left as an exercise. Your final code should look like
#inc lude
#inc lude
us ing namespace std ;
double normalDi s t r ibut ion ( double x )
{
s t a t i c const double RT2PI = sq r t ( 4 . 0∗ acos ( 0 . 0 ) ) ;
s t a t i c const double SPLIT = 10 ./ sq r t (2 ) ;
s t a t i c const double a [ ] =
{220 .206867912376 ,221 .213596169931 ,112 .079291497871 ,33 .912866078383 ,6 .37396220353165 ,0 .700383064443688 ,3 .52624965998911
e−02};
s t a t i c const double b [ ] =
{440 .413735824752 ,793 .826512519948 ,637 .333633378831 ,296 .564248779674 ,86 .7807322029461 ,16 .064177579207 ,1 .75566716318264 ,8 .83883476483184
e−02};
const double z = fabs (x ) ;
double Nz = 0 . 0 ;
// i f z ou t s id e these l im i t s then value e f f e c t i v e l y 0 or 1 f o r machine p r e c i s i o n
i f ( z<=37.0)
{
// NDash = N ’ ( z ) ∗ s q r t {2\ pi }
const double NDash = exp(−z∗z /2 . 0 ) /RT2PI ;
i f ( z{
const double Pz = ( ( ( ( ( a [ 6 ] ∗ z + a [ 5 ] ) ∗z + a [ 4 ] ) ∗z + a [ 3 ] ) ∗z + a [ 2 ] ) ∗z + a [ 1 ] )
∗z + a [ 0 ] ;
const double Qz = ( ( ( ( ( ( b [ 7 ] ∗ z + b [ 6 ] ) ∗z + b [ 5 ] ) ∗z + b [ 4 ] ) ∗z + b [ 3 ] ) ∗z + b
[ 2 ] ) ∗z + b [ 1 ] ) ∗z + b [ 0 ] ;
Nz = RT2PI∗NDash∗Pz/Qz ;
}
e l s e
{
const double F4z = z + 1 .0/ ( z + 2 . 0/ ( z + 3 . 0/ ( z + 4 . 0/ ( z + 13 . 0/20 . 0 ) ) ) ) ;
Nz = NDash/F4z ;
17
}
}
re turn x>=0.0 ? 1−Nz : Nz ;
}
// return the value o f a c a l l opt ion us ing the black s cho l e s formula
double ca l lOpt i onPr i c e ( double S , double t , double X, double r , double sigma , double T)
{
i f ( f abs (T−t )<1.e−14) // check i f we are at maturity
{
i f (Se l s e re turn S−X;
}
i f ( (T−t )<=−1.e−14) re turn 0 . ; // opt ion exp i red
i f (X<1.e−14∗S) re turn S−X∗exp(−r ∗(T−t ) ) ; // check i f s t r i k e << a s s e t then e x e r c i s e
with c e r t a i n t y
i f (S<1.e−14∗X) return 0 . ; // check i f a s s e t << s t r i k e then wor th l e s s
i f ( sigma∗ sigma ∗(T−t )<1.e−14) // check i f var i ance very smal l then no d i f f u s i o n
{
i f (Se l s e re turn S−X∗exp(−r ∗(T−t ) ) ;
}
// c a l c u l a t e opt ion p r i c e
double d1=( log (S/X) + ( r+sigma∗ sigma /2 . ) ∗(T−t ) ) /( sigma∗ s q r t (T−t ) ) ;
double d2=( log (S/X) + ( r−sigma∗ sigma /2 . ) ∗(T−t ) ) /( sigma∗ s q r t (T−t ) ) ;
r e turn normalDi s t r ibut ion ( d1 ) ∗S − normalDi s t r ibut ion ( d2 ) ∗X∗exp(−r ∗(T−t ) ) ;
}
i n t main ( )
{
cout << ”Cal l Option Pr i ce = ” << ca l lOpt i onPr i c e ( 1 , 0 , 1 , 0 . 0 5 , 0 . 2 , 1 ) << endl ;
r e turn 0 ;
}
download
Similarly for a put option we would have
#inc lude
#inc lude
us ing namespace std ;
double normalDi s t r ibut ion ( double x )
{
s t a t i c const double RT2PI = sq r t ( 4 . 0∗ acos ( 0 . 0 ) ) ;
s t a t i c const double SPLIT = 10 ./ sq r t (2 ) ;
s t a t i c const double a [ ] =
{220 .206867912376 ,221 .213596169931 ,112 .079291497871 ,33 .912866078383 ,6 .37396220353165 ,0 .700383064443688 ,3 .52624965998911
e−02};
s t a t i c const double b [ ] =
{440 .413735824752 ,793 .826512519948 ,637 .333633378831 ,296 .564248779674 ,86 .7807322029461 ,16 .064177579207 ,1 .75566716318264 ,8 .83883476483184
e−02};
const double z = fabs (x ) ;
double Nz = 0 . 0 ;
// i f z ou t s id e these l im i t s then value e f f e c t i v e l y 0 or 1 f o r machine p r e c i s i o n
18
i f ( z<=37.0)
{
// NDash = N ’ ( z ) ∗ s q r t {2\ pi }
const double NDash = exp(−z∗z /2 . 0 ) /RT2PI ;
i f ( z{
const double Pz = ( ( ( ( ( a [ 6 ] ∗ z + a [ 5 ] ) ∗z + a [ 4 ] ) ∗z + a [ 3 ] ) ∗z + a [ 2 ] ) ∗z + a [ 1 ] )
∗z + a [ 0 ] ;
const double Qz = ( ( ( ( ( ( b [ 7 ] ∗ z + b [ 6 ] ) ∗z + b [ 5 ] ) ∗z + b [ 4 ] ) ∗z + b [ 3 ] ) ∗z + b
[ 2 ] ) ∗z + b [ 1 ] ) ∗z + b [ 0 ] ;
Nz = RT2PI∗NDash∗Pz/Qz ;
}
e l s e
{
const double F4z = z + 1 .0/ ( z + 2 . 0/ ( z + 3 . 0/ ( z + 4 . 0/ ( z + 13 . 0/20 . 0 ) ) ) ) ;
Nz = NDash/F4z ;
}
}
re turn x>=0.0 ? 1−Nz : Nz ;
}
// return the value o f a put opt ion us ing the black s cho l e s formula
double putOptionPrice ( double S , double t , double X, double r , double sigma , double T)
{
i f ( f abs (T−t )<1.e−14) // check i f we are at maturity
{
i f (Se l s e re turn 0 ;
}
i f ( (T−t )<=−1.e−14) re turn 0 . ; // opt ion exp i red
i f (X<1.e−14∗S) re turn 0 . ; // check i f s t r i k e << a s s e t then e x e r c i s e with c e r t a i n t y
i f (S<1.e−14∗X) return X∗exp(−r ∗(T−t ) ) − S ; // check i f a s s e t << s t r i k e then
wor th l e s s
i f ( sigma∗ sigma ∗(T−t )<1.e−14) // check i f var i ance very smal l then no d i f f u s i o n
{
i f (Se l s e re turn 0 . ;
}
// c a l c u l a t e opt ion p r i c e
double d1=( log (S/X) + ( r+sigma∗ sigma /2 . ) ∗(T−t ) ) /( sigma∗ s q r t (T−t ) ) ;
double d2=( log (S/X) + ( r−sigma∗ sigma /2 . ) ∗(T−t ) ) /( sigma∗ s q r t (T−t ) ) ;
r e turn normalDi s t r ibut ion (−d2 ) ∗X∗exp(−r ∗(T−t ) ) − normalDi s t r ibut ion (−d1 ) ∗S ;
}
i n t main ( )
{
cout << ”Put Option Pr i ce = ” << putOptionPrice ( 1 , 0 , 1 , 0 . 0 5 , 0 . 2 , 1 ) << endl ;
r e turn 0 ;
}
download
Other Options
For a binary call we get
19
#inc lude
#inc lude
us ing namespace std ;
double normalDi s t r ibut ion ( double x )
{
s t a t i c const double RT2PI = sq r t ( 4 . 0∗ acos ( 0 . 0 ) ) ;
s t a t i c const double SPLIT = 10 ./ sq r t (2 ) ;
s t a t i c const double a [ ] =
{220 .206867912376 ,221 .213596169931 ,112 .079291497871 ,33 .912866078383 ,6 .37396220353165 ,0 .700383064443688 ,3 .52624965998911
e−02};
s t a t i c const double b [ ] =
{440 .413735824752 ,793 .826512519948 ,637 .333633378831 ,296 .564248779674 ,86 .7807322029461 ,16 .064177579207 ,1 .75566716318264 ,8 .83883476483184
e−02};
const double z = fabs (x ) ;
double Nz = 0 . 0 ;
// i f z ou t s id e these l im i t s then value e f f e c t i v e l y 0 or 1 f o r machine p r e c i s i o n
i f ( z<=37.0)
{
// NDash = N ’ ( z ) ∗ s q r t {2\ pi }
const double NDash = exp(−z∗z /2 . 0 ) /RT2PI ;
i f ( z{
const double Pz = ( ( ( ( ( a [ 6 ] ∗ z + a [ 5 ] ) ∗z + a [ 4 ] ) ∗z + a [ 3 ] ) ∗z + a [ 2 ] ) ∗z + a [ 1 ] )
∗z + a [ 0 ] ;
const double Qz = ( ( ( ( ( ( b [ 7 ] ∗ z + b [ 6 ] ) ∗z + b [ 5 ] ) ∗z + b [ 4 ] ) ∗z + b [ 3 ] ) ∗z + b
[ 2 ] ) ∗z + b [ 1 ] ) ∗z + b [ 0 ] ;
Nz = RT2PI∗NDash∗Pz/Qz ;
}
e l s e
{
const double F4z = z + 1 .0/ ( z + 2 . 0/ ( z + 3 . 0/ ( z + 4 . 0/ ( z + 13 . 0/20 . 0 ) ) ) ) ;
Nz = NDash/F4z ;
}
}
re turn x>=0.0 ? 1−Nz : Nz ;
}
// return the value o f a put opt ion us ing the black s cho l e s formula
double b inaryCal lOpt ionPr ice ( double S , double t , double X, double r , double sigma ,
double T)
{
i f ( f abs (T−t )<1.e−14) // check i f we are at maturity
{
i f (S>X) return 1 ;
e l s e re turn 0 ;
}
i f ( (T−t )<=−1.e−14) re turn 0 . ; // opt ion exp i red
i f (X<1.e−14∗S) re turn exp(−r ∗(T−t ) ) ; // check i f s t r i k e << a s s e t then e x e r c i s e
with c e r t a i n t y
i f (S<1.e−14∗X) return 0 . ; // check i f a s s e t << s t r i k e then wor th l e s s
i f ( sigma∗ sigma ∗(T−t )<1.e−14) // check i f var i ance very smal l then no d i f f u s i o n
{
i f (S>X∗exp(−r ∗(T−t ) ) ) re turn exp(−r ∗(T−t ) ) ;
20
e l s e re turn 0 . ;
}
// c a l c u l a t e opt ion p r i c e
double d2=( log (S/X) + ( r−sigma∗ sigma /2 . ) ∗(T−t ) ) /( sigma∗ s q r t (T−t ) ) ;
r e turn normalDi s t r ibut ion ( d2 ) ∗exp(−r ∗(T−t ) ) ;
}
i n t main ( )
{
cout << ”Binary Ca l l Option Pr i ce = ” << binaryCal lOpt ionPr ice ( 1 , 0 , 1 , 0 . 0 5 , 0 . 2 , 1 )
<< endl ;
r e turn 0 ;
}
download
and for a binary put
#inc lude
#inc lude
us ing namespace std ;
double normalDi s t r ibut ion ( double x )
{
s t a t i c const double RT2PI = sq r t ( 4 . 0∗ acos ( 0 . 0 ) ) ;
s t a t i c const double SPLIT = 10 ./ sq r t (2 ) ;
s t a t i c const double a [ ] =
{220 .206867912376 ,221 .213596169931 ,112 .079291497871 ,33 .912866078383 ,6 .37396220353165 ,0 .700383064443688 ,3 .52624965998911
e−02};
s t a t i c const double b [ ] =
{440 .413735824752 ,793 .826512519948 ,637 .333633378831 ,296 .564248779674 ,86 .7807322029461 ,16 .064177579207 ,1 .75566716318264 ,8 .83883476483184
e−02};
const double z = fabs (x ) ;
double Nz = 0 . 0 ;
// i f z ou t s id e these l im i t s then value e f f e c t i v e l y 0 or 1 f o r machine p r e c i s i o n
i f ( z<=37.0)
{
// NDash = N ’ ( z ) ∗ s q r t {2\ pi }
const double NDash = exp(−z∗z /2 . 0 ) /RT2PI ;
i f ( z{
const double Pz = ( ( ( ( ( a [ 6 ] ∗ z + a [ 5 ] ) ∗z + a [ 4 ] ) ∗z + a [ 3 ] ) ∗z + a [ 2 ] ) ∗z + a [ 1 ] )
∗z + a [ 0 ] ;
const double Qz = ( ( ( ( ( ( b [ 7 ] ∗ z + b [ 6 ] ) ∗z + b [ 5 ] ) ∗z + b [ 4 ] ) ∗z + b [ 3 ] ) ∗z + b
[ 2 ] ) ∗z + b [ 1 ] ) ∗z + b [ 0 ] ;
Nz = RT2PI∗NDash∗Pz/Qz ;
}
e l s e
{
const double F4z = z + 1 .0/ ( z + 2 . 0/ ( z + 3 . 0/ ( z + 4 . 0/ ( z + 13 . 0/20 . 0 ) ) ) ) ;
Nz = NDash/F4z ;
}
}
re turn x>=0.0 ? 1−Nz : Nz ;
}
21
// return the value o f a put opt ion us ing the black s cho l e s formula
double binaryPutOptionPrice ( double S , double t , double X, double r , double sigma , double
T)
{
i f ( f abs (T−t )<1.e−14) // check i f we are at maturity
{
i f (Se l s e re turn 0 ;
}
i f ( (T−t )<=−1.e−14) re turn 0 . ; // opt ion exp i red
i f (X<1.e−14∗S) re turn 0 . ; // check i f s t r i k e << a s s e t then e x e r c i s e with c e r t a i n t y
i f (S<1.e−14∗X) return exp(−r ∗(T−t ) ) ; // check i f a s s e t << s t r i k e then wor th l e s s
i f ( sigma∗ sigma ∗(T−t )<1.e−14) // check i f var i ance very smal l then no d i f f u s i o n
{
i f (Se l s e re turn 0 . ;
}
// c a l c u l a t e opt ion p r i c e
double d2=( log (S/X) + ( r−sigma∗ sigma /2 . ) ∗(T−t ) ) /( sigma∗ s q r t (T−t ) ) ;
r e turn normalDi s t r ibut ion (−d2 ) ∗exp(−r ∗(T−t ) ) ;
}
i n t main ( )
{
cout << ”Binary Put Option Pr i ce = ” << binaryPutOptionPrice ( 1 , 0 , 1 , 0 . 0 5 , 0 . 2 , 1 ) <<
endl ;
r e turn 0 ;
}
download
1.3 Coursework Example
A trader has asked you to price the value of the call option C(S, t) at time t = 0 according to the
standard Black-Scholes formula, where T = 1, X = 1, r = 0.05 and σ = 0.2. Write a program to
calculate C and output the results to screen. You must generate four columns of data:
• the value of S;
• the value of d1;
• the value of d1;
• the value of C(S, t = 0).
Output each of the values when the stock price is
S ∈ {0.8, 0.9, 1, 1.1, 1.2}.
You should use a for loop to generate the data.
22
We could do this with the call option code above, but to make outputting d1 and d2 easier lets
just calculate them again. First, open up a new project and copy in the code from earlier on with
normal distribution, get it by (clicking here). Then simply start by editing the main function so
that has variables for all the required inputs (S, t, T , X, r and σ), and then calculate d1 and N(d1).
i n t main ( )
{
double S=1. , t =0. ,X=1. , r =0.05 , sigma=0.2 ,T=1. ;
double d 1=( log (S/X) + ( r+sigma∗ sigma /2 . ) ∗(T−t ) ) / sigma/ sq r t (T−t ) ;
cout << ”d 1 = ” << d 1 << endl ;
cout << ”N( d 1 ) = ” << normalDi s t r ibut ion ( d 1 ) << endl ;
r e turn 0 ;
}
download
Once you have checked this result, do the same for d2 and N(d2) then calculate the call option
price
i n t main ( )
{
double S=1. , t =0. ,X=1. , r =0.05 , sigma=0.2 ,T=1. ;
double d 1=( log (S/X) + ( r+sigma∗ sigma /2 . ) ∗(T−t ) ) / sigma/ sq r t (T−t ) ;
double d 2=( log (S/X) + ( r−sigma∗ sigma /2 . ) ∗(T−t ) ) / sigma/ sq r t (T−t ) ;
cout << ”d 1 = ” << d 1 << endl ;
cout << ”N( d 1 ) = ” << normalDi s t r ibut ion ( d 1 ) << endl ;
cout << ”d 2 = ” << d 2 << endl ;
cout << ”N( d 2 ) = ” << normalDi s t r ibut ion ( d 2 ) << endl ;
cout << ”C = ” << S∗ normalDi s t r ibut ion ( d 1 ) − X∗exp(−r ∗(T−t ) ) ∗ normalDi s t r ibut ion (
d 2 ) << endl ;
r e turn 0 ;
}
download
If you are happy everything is working then run a loop to calculate values for different values
of S, something like this should work
i n t main ( )
{
double S=1. , t =0. ,X=1. , r =0.05 , sigma=0.2 ,T=1. ;
f o r ( i n t i =8; i <=12; i++)
{
S = i ∗ 0 . 1 ;
cout << S ;
double d 1=( log (S/X) + ( r+sigma∗ sigma /2 . ) ∗(T−t ) ) / sigma/ sq r t (T−t ) ;
double d 2=( log (S/X) + ( r−sigma∗ sigma /2 . ) ∗(T−t ) ) / sigma/ sq r t (T−t ) ;
cout << ” ” << d 1 ;
cout << ” ” << d 2 ;
23
cout << ” ” << S∗ normalDi s t r ibut ion ( d 1 ) − X∗exp(−r ∗(T−t ) ) ∗
normalDi s t r ibut ion ( d 2 ) << endl ;
}
re turn 0 ;
}
download
Now you have your solution, export this manually (from terminal output) or using a filestream
in c++ so that you can create a table. Your final result should look something like this:-
References
West, Graeme. Better approximations to cumulative normal functions. Wilmott Magazine, 9:
70–76, 2005. URL https://lyle.smu.edu/~aleskovs/emis/sqc2/accuratecumnorm.pdf.
24

学霸联盟