LINEBURG

 << Пред. стр. страница 17(всего 20)ОГЛАВЛЕНИЕ След. стр. >>
j
r
rr0 rud = udr0
ВЁ
B
ВЁ
rr ВЁ
ВЁ
r ВЁ
rr ВЁ
r ВЁВЁ
r rd = dr0
j
rr
rr
r
rr
rdd = ddr0
r
j

The п¬Ѓgure illustrates the building of an interest rate tree of one period spot rates by assuming that for any given period t the next
period interest rate can only take on two values, rt+1 = urt or rt+1 = drt , where u and d are constants. r0 is the initial spot rate.

Code 20.1 shows how one can construct such an interest rate tree.
Code 20.2 implements the original algorithm for a call option on a (long maturity) zero coupon bond.

124
#include <vector>
#include <cmath>
using namespace std;

vector< vector<double> >
build interest rate tree rendleman bartter(const double& r0,
const double& u,
const double& d,
const int& n){
vector< vector<double> > tree;
for (int i=1;i<=n;++i){
vector<double> r(i);
for (int j=0;j<i;++j){
r[j] = r0*pow(u,j)*pow(d,iв€’jв€’1);
};
tree.push back(r);
};
return tree;
};

Code 20.1: Building an interest rate tree

General references include Sundaresan .
Rendleman and Bartter  and Rendleman and Bartter  are the original references for building
standard binomial interest rate trees.

125
#include <cmath>
#include <algorithm>
#include <vector>
using namespace std;

double bond option price call zero american rendleman bartter(const double& X,
const double& option maturity,
const double& S,
const double& M, // term structure paramters
const double& interest, // current short interest rate
const double& bond maturity, // time to maturity for underlying bo
const double& maturity payment,
const int& no steps) {
double delta t = bond maturity/no steps;

double u=exp(S*sqrt(delta t));
double d=1/u;
double p up = (exp(M*delta t)в€’d)/(uв€’d);
double p down = 1.0в€’p up;

vector<double> r(no steps+1);
r=interest*pow(d,no steps);
double uu=u*u;
for (int i=1;i<=no steps;++i){ r[i]=r[iв€’1]*uu;};
vector<double> P(no steps+1);
for (int i=0;i<=no steps;++i){ P[i] = maturity payment; };
int no call steps=int(no steps*option maturity/bond maturity);
for (int curr step=no steps;curr step>no call steps;в€’в€’curr step) {
for (int i=0;i<curr step;i++) {
r[i] = r[i]*u;
P[i] = exp(в€’r[i]*delta t)*(p down*P[i]+p up*P[i+1]);
};
};
vector<double> C(no call steps+1);
for (int i=0;i<=no call steps;++i){ C[i]=max(0.0,P[i]в€’X); };
for (int curr step=no call steps;curr step>=0;в€’в€’curr step) {
for (int i=0;i<curr step;i++) {
r[i] = r[i]*u;
P[i] = exp(в€’r[i]*delta t)*(p down*P[i]+p up*P[i+1]);
C[i]=max(P[i]в€’X, exp(в€’r[i]*delta t)*(p up*C[i+1]+p down*C[i]));
};
};
return C;
};

Code 20.2: RB binomial model for European call on zero coupon bond

126
Chapter 21

Term Structure Derivatives
Contents
21.1 Vasicek bond option pricing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127

21.1 Vasicek bond option pricing

If the term structure model is VasicekвЂ™s model there is a solution for the price of an option on a zero coupon
bond, due to Jamshidan .
Under VacisekвЂ™s model the process for the short rate is assumed to follow.

dr = a(b в€’ r)dt + ПѓdZ

where a, b and Пѓ are constants. We have seen earlier how to calculate the discount factor in this case. We
now want to consider an European Call option in this setting.
Let P (t, s) be the time t price of a zero coupon bond with a payment of \$1 at time s (the discount factor).
The price at time t of a European call option maturing at time T on on a discount bond maturing at time
s is( See Jamshidan  and Hull )

P (t, s)N (h) в€’ XP (t, T )N (h в€’ ПѓP )

where
1 P (t, s) 1
h= ln + ПѓP
ПѓP P (t, T )X 2

ПѓP = v(t, T )B(T, s)

1 в€’ eв€’a(T в€’t)
B(t, T ) =
a

Пѓ 2 (1 в€’ eв€’a(T в€’t) )
2
v(t, T ) =
2a

In the case of a = 0,
в€љ
v(t, T ) = Пѓ T в€’ t
в€љ
ПѓP = Пѓ(s в€’ T ) T в€’ t

127
#include "normdist.h"
#include "fin_recipes.h"
#include <cmath>

using namespace std;

double bond option price call zero vasicek(const double& X, // exercise price
const double& r, // current interest rate
const double& option time to maturity,
const double& bond time to maturity,
const double& a, // parameters
const double& b,
const double& sigma){
double T t = option time to maturity;
double s t = bond time to maturity;
double T s = s tв€’T t;
double v t T;
double sigma P;
if (a==0.0) {
v t T = sigma * sqrt ( T t ) ;
sigma P = sigma*T s*sqrt(T t);
}
else {
v t T = sqrt (sigma*sigma*(1в€’exp(в€’2*a*T t))/(2*a));
double B T s = (1в€’exp(в€’a*T s))/a;
sigma P = v t T*B T s;
};
double h = (1.0/sigma P) * log (
term structure discount factor vasicek(s t,r,a,b,sigma)/
(term structure discount factor vasicek(T t,r,a,b,sigma)*X) )
+ sigma P/2.0;
double c =
term structure discount factor vasicek(s t,r,a,b,sigma)*N(h)
в€’X*term structure discount factor vasicek(T t,r,a,b,sigma)*N(hв€’sigma P);
return c;
};

Code 21.1:

128
Appendix A

Normal Distribution approximations.
Contents
A.1 The normal distribution function . . . . . . . . . . . . . . . . . . . . . . . . . . . 129
A.2 The cumulative normal distribution . . . . . . . . . . . . . . . . . . . . . . . . . 129
A.3 Multivariate normal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130
A.4 Calculating cumulative bivariate normal probabilities . . . . . . . . . . . . . . 130
A.5 Simulating random normal numbers . . . . . . . . . . . . . . . . . . . . . . . . . 132
A.6 Cumulative probabilities for general multivariate distributions . . . . . . . . . 132
A.7 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132
A.8 Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133

We will in general not go into detail about more standard numerical problems not connected to п¬Ѓnance,
there are a number of well known sources for such, but we show the example of calculations involving the
normal distribution.

A.1 The normal distribution function

The nurmal distribution function
x2
n(x) = eв€’ 2

is calculated as

#include <cmath> // c library of math functions
using namespace std; // which is part of the standard namespace

// most C compilers deп¬Ѓne PI, but just in case it doesnвЂ™t
#ifndef PI
#deп¬Ѓne PI 3.141592653589793238462643
#endif

double n(double z) { // normal distribution function
return (1.0/sqrt(2.0*PI))*exp(в€’0.5*z*z);
};

Code A.1: The normal distribution function

A.2 The cumulative normal distribution

The solution of a large number of option pricing formulas are written in terms of the cumulative normal
distribution. For a random variable x the cumulative probability is the probability that the outcome is lower
than a given value z. To calculate the probability that a normally distubuted random variable with mean 0
and unit variance is less than z, N (z), one have to evaluate the integral
z z
x2
eв€’
Prob(x в‰¤ z) = N (z) = n(x)dx = dx
2

в€’в€ћ в€’в€ћ

129
There is no explicit closed form solution for calculation of this integral, but a large number of well known
approximations exists. Abramowiz and Stegun  is a good source for these approximations. The
following is probably the most used such approximation, it being pretty accurate and relatively fast. The
arguments to the function are assumed normalized to a (0,1 ) distribution.

#include <cmath> // math functions.
using namespace std;

double N(double z) {
if (z > 6.0) { return 1.0; }; // this guards against overп¬‚ow
if (z < в€’6.0) { return 0.0; };

double b1 = 0.31938153;
b2 = в€’0.356563782;
double
double b3 = 1.781477937;
b4 = в€’1.821255978;
double
double b5 = 1.330274429;
double p = 0.2316419;
double c2 = 0.3989423;

double a=fabs(z);
double t = 1.0/(1.0+a*p);
double b = c2*exp((в€’z)*(z/2.0));
double n = ((((b5*t+b4)*t+b3)*t+b2)*t+b1)*t;
n = 1.0в€’b*n;
if ( z < 0.0 ) n = 1.0 в€’ n;
return n;
};

Code A.2: The cumulative normal

A.3 Multivariate normal

The normal distribution is also deп¬Ѓned for several random variables. We then characterise the vector of
random variables
пЈ® пЈ№
x1
пЈЇ x2 пЈє
X=пЈЇ . пЈє
пЈЇ пЈє
пЈ°.пЈ» .
xn

A.4 Calculating cumulative bivariate normal probabilities

The most used multivariate normal calculation is the bivariate case, where we let x and y be bivariate
normally distributed, each with mean 0 and variance 1, and assume the two variables have correlation of ПЃ.
By the deп¬Ѓnition of correlation ПЃ в€€ [в€’1, 1]. The cumulative probability distribution

P (x < a, y < b) = N (a, b, ПЃ)
a b
1 x2 в€’ 2ПЃxy + y 2
1
exp в€’
= dxdy
1 в€’ ПЃ2
2
1 в€’ ПЃ2
2ПЂ
в€’в€ћ в€’в€ћ

There are several approximations to this integral. We pick one such, discussed in [Hull, 1993, Ch 10], shown
in code A.3

130
#include <cmath> // include the standard library mathematics functions
using namespace std; // which are in the standard namespace

double N(double); // deп¬Ѓne the univariate cumulative normal distribution as a separate function

#ifndef PI
const double PI=3.141592653589793238462643;
#endif

inline double f(double x, double y, double aprime, double bprime, double rho) {
double r = aprime*(2*xв€’aprime) + bprime*(2*yв€’bprime) + 2*rho*(xв€’aprime)*(yв€’bprime);
return exp(r);
};

inline double sgn( double x) { // sign function
if (x>=0.0) return 1.0;
return в€’1.0;
};

double N(double a, double b, double rho) {
// Numerical approximation to the bivariate normal distribution,
// as described e.g. in Hulls book
if ( (a<=0.0) && (b<=0.0) && (rho<=0.0) ) {
double aprime = a/sqrt(2.0*(1.0в€’rho*rho));
double bprime = b/sqrt(2.0*(1.0в€’rho*rho));
double A={0.3253030, 0.4211071, 0.1334425, 0.006374323};
double B={0.1337764, 0.6243247, 1.3425378, 2.2626645 };
double sum = 0;
for (int i=0;i<4;i++) {
for (int j=0; j<4; j++) {
sum += A[i]*A[j]* f(B[i],B[j],aprime,bprime,rho);
};
};
sum = sum * ( sqrt(1.0в€’rho*rho)/PI);
return sum;
}
else if ( a * b * rho <= 0.0 ) {
if ( ( a<=0.0 ) && ( b>=0.0 ) && (rho>=0.0) ) {
return N(a) в€’ N(a, в€’b, в€’rho);
}
else if ( (a>=0.0) && (b<=0.0) && (rho>=0.0) ) {
return N(b) в€’ N(в€’a, b, в€’rho);
}
else if ( (a>=0.0) && (b>=0.0) && (rho<=0.0) ) {
return N(a) + N(b) в€’ 1.0 + N(в€’a, в€’b, rho);
};
}
else if ( a * b * rho >= 0.0 ) {
double denum = sqrt(a*a в€’ 2*rho*a*b + b*b);
double rho1 = ((rho * a в€’ b) * sgn(a))/denum;
double rho2 = ((rho * b в€’ a) * sgn(b))/denum;
double delta=(1.0в€’sgn(a)*sgn(b))/4.0;
return N(a,0.0,rho1) + N(b,0.0,rho2) в€’ delta;
};
 << Пред. стр. страница 17(всего 20)ОГЛАВЛЕНИЕ След. стр. >>