Chapter 1
1
CHAPTER 9
Monte Carlo Option Pricings
9.1. The Monte Carlo Method
The Monte Carlo method provides numerical solution to a variety of mathematical problems by
performing statistical samplings on a computer. In risk-neutral pricing of options, we are most
interesting in evaluating the expected value of a function g(x) under a random variable x as
E( g(x) )
dx g(x) (x) (9.1)
where (x) is the probability density function of x. In general, it would be difficult to derive an
analytical formula for (9.1) and a numerical estimation seems to be the only way out. Monte Carlo
simulation provides a simple and flexible numerical method to solve these types of problems1.
Consider the random sample { x1, x2, … , xn } generated based on the probability density function
(x). The estimates of the mean and variance of g(x) are given by
and (9.2)
According to the central limit theorem, the random variable defined as
(9.3)
tends to follow a standard normal distribution with increasing sample size n and be irrespective of the
distribution of g(x). Thus, the sample average m approaches a normal distribution with mean E( g(x) )
and standard deviation (s/√n). On this basis, the confidence interval in the estimation of E( g(x) ) can
be obtained as
(9.4)
with, for example, confidence level of 68.27% for z 1. We refer to the term (s/√n) in (9.4) as the
standard error in the estimation of E( g(x) ). To reduce the standard error by a factor of ten, the sample
size has to be increased one hundredfold. The method is thus computationally inefficient in its
fundamental form called the crude Monte Carlo simulation.
There are variance reduction techniques that focus on reducing instead the size of s in the standard
error. The common approach is known as the control variate method. It takes the analytic solution of a
similar but simpler problem to improve the accuracy in the estimated solution of a complex problem.
Suppose that the expected value E( h(x) ) can be evaluated analytically as H. In related to the original
function g(x), we can introduce a new function given by
g̃(x) g(x) h(x) (9.5)
through the control variate h(x) and rewrite (9.1) as
E( g(x) ) H
dx g̃(x) (x) (9.6)
1 See P. P. Boyle, “Options : A Monte Carlo Approach”, Journal of Financial Economics, Vol. 4, Issue 3 (1977)
p 323-338.
m E( g(x) )
s
√n
s
√n
E( g(x) ) m z
m
n
i 1 g(xi)
1
n
s2
n
i 1 ( g(xi) m )
2
1
n 1
2
Thus, we can determine the confidence interval in the estimation of E( g(x) ) based on instead the
estimates of the mean and variance of g̃ (x) given by
(9.7)
It can be shown that the variances of g(x) and g̃(x) can be related as
var( g̃ (x) ) var( g(x) ) var( h(x) ) 2 cov( g(x) , h(x) ) (9.8)
If g(x) and h(x) are similar problems, the covariance between them is positive. In (9.8), the variance of
g̃(x) will be less than the variance of g(x) as long as cov( g(x) , h(x) ) ½ var( h(x) ). It is therefore
possible to reduce the size of the standard error by identifying a highly correlated problem with
known analytic solution.
An alternative approach is known as the antithetic variate method. In case of standard normal
variable, it makes use of the symmetric property around zero in the density function. Again, we can
introduce a new function given by
ĝ(x) ½[ g(x) g(x) ] (9.9)
through antithetic variate of the form x. We can rewrite (9.1) using the symmetric property of the
standard normal variable x as
E( g(x) )
dx ĝ(x) (x) (9.10)
Similarly, we can determine the confidence interval in the estimation of E( g(x) ) based on the
estimates of the mean and variance of ĝ (x) given by
(9.11)
The variance of ĝ(x) is expected to be smaller than the variance of g(x) as it is already an average
quantity of two samples. It can be shown that the two variances can be related as
var( ĝ (x) ) ½ var( g(x) ) ½ cov( g(x) , g(x) ) (9.12)
If the covariance between g(x) and g(x) is negative, it is always efficient to consider the estimates for
ĝ(x) rather than doubling the size of independent samples.
Figure 9.1: Monte Carlo simulation for E( e
x
).
Figure 9.1 depicts, for example, the Monte Carlo estimation of the expected value E( e
x
) with x
taken to be a standard normal variable. In EXCEL, random samples of x ( 0 , 1 ) can be generated
ŝ
√n
E( g(x) ) m̂ z
s̃
√n
E( g(x) ) ( H m̃ ) z
3
very easily by calling the function NORMSINV(RAND()). They will be used to evaluate the sample
values of g(x) e
x
for which mean and variance can be estimated. In the figure, we compare the
confidence interval evaluated through the crude simulation with those based on the variance reduction
techniques. For the control variate method, we have adopted a highly correlated problem E( x ) with
known analytic solution of zero. Thus, we introduce g̃(x) e
x
x for which E( e
x
) E( e
x
x ). For
the antithetic variate method, it is clear that the covariance between e
x
and e
x
is negative. We thus
define ĝ(x) ½ ( e
x
e
x
) with E( e
x
) E( ½ ( e
x
e
x
) ). In both cases, we have shown that there are
significant reductions in the size of the standard errors.
9.2. Risk-Neutral Valuation
Under risk-neutral valuation, current price of an option can be defined based on the present value
of its average maturity payoff at time T as
f0 Ê( e
rT fT | S0 ) (9.13)
where r is the constant interest rate. Here, we are averaging over realized maturity payoffs of the
option fT in respect to sample underlying asset prices generated through its risk-neutral process that
initiated at current price S0. In stochastic model, asset price return during the time increment from t to
t t is assumed to follow a random normal process as
St / St t
t ( 0 , 1 ) (9.14)
where and are respectively the mean rate and volatility of return. For traded asset such as stock,
the risk-neutral process is simply given by (9.14) with replaced by r in the drift term. Practically, it
is convenient to consider the asset price movement based on the risk-neutral process. For constant
volatility of return, it is shown to follow an iterative equation with arbitrary time duration given by
St St exp( ( r ½2 )
( 0 , 1 ) ) (9.15)
In particular, we have
ST S0 exp( ( r ½2 )T
T ( 0 , 1 ) ) (9.16)
that generates the maturity price ST directly from S0 in a single step. It will be useful when there is no
intermediate boundary condition for the option.
As our first example, we consider a European call option written on a non-dividend paying stock
with strike price K. It should be noted that there exits analytic solution for this option known as the
Black-Scholes formula given by
f0
BS
S0 N( d ) K e
rT N( d
T ) , (9.17)
Here, we take this simple example to demonstrate the use of Monte Carlo procedure for option
pricing. In this case, the option payoff will depend solely on the underlying asset price at maturity
regardless of its intermediate values. We can simply use (9.16) to generate the maturity price of the
asset by a single random number from ( 0 , 1 ). The sample maturity price can then be used to
evaluate the sample maturity payoff of the option according to the function
d
ln(S0 /K) ( r ½
2 )T
T
4
fT()
max{ ST() K , 0 } (9.18)
For variance reduction, we can adopt the maturity price
ST itself as the control variate and consider the
new function
fT̃()
max{ ST() K , 0 }
ST()
(9.19)
The control variate has an analytic solution given by Ê( e
rT
ST | S0 ) S0 in which we can rewrite
(9.13) as
f0 S0 Ê( e
rT fT̃ | S0 ) (9.20)
Alternatively, we can simply take as the antithetic variate and introduce the new function
fT̂() ½ [
max{ ST() K , 0 }
max{ ST() K , 0 } ] (9.21)
Figure 9.2 depicts the Monte Carlo results for the current price of the European call option with
parameters defined in cells B2:B6. Random samples of are first generated in EXCEL spreadsheet,
and they will be used to generate the sample maturity prices ST() as well as their antithetic values
ST(). Sample values of the functions fT(), fT̃(), and f̂T() can then be evaluated for which mean and
variance can be estimated. As shown in the figure, there are significant reductions in the size of the
standard errors with the use of the variance reduction techniques. As reference, the Black-Scholes
pricing in (9.17) is 12.34 for the European call option parameterized by B2:B6.
Figure 9.2 : Monte Carlo simulation for the current price of a European call option.
An important application of Monte Carlo method is in the pricing of exotic options with
intermediate boundary conditions. For example, a barrier knock-out option will have the maturity
payoff depending on whether the underlying asset price has surpassed a predetermined barrier level
prior to its maturity. Thus, it is sometimes necessary to generate the entire path of asset prices in
discrete time increments { t1, … , tN T } during the life of the option. To this end, we can use (9.15)
to generate iteratively the risk-neutral asset price at ti 1 based on that at ti as
Sti 1 Sti exp( ( r ½
2 )( ti 1 ti )
ti 1 ti i 1 ) (9.22)
and trace out the path by running i from 0 to N 1 with random numbers { 1, … , N } from ( 0 , 1 ).
The sample path can then be used to evaluate the sample option payoff according to some path-
dependent conditions. It should be noted that we have defined in the iteration t0 0 and the choice of
time increments will depend on the boundary conditions of the option.
5
9.3. VBA Implementation
A crucial part of the Monte Carlo method is the generation of the standard normal random
numbers. The procedure starts with a pseudo-random number generator (PRNG) that “randomly”
produces a real number between 0 and 1 with uniform probability so that every number will have the
same chance of being generated. The numbers generated by PRNG are not truly random in that they
are completely determined by an initial seed integer. In fact, they are considered to be random only
subject to standard statistical tests for randomness2. Uniform random numbers can be transformed into
standard normal random numbers through Box-Muller algorithm as
1
2ln(u1) cos( 2u2 ) (9.23)
2
2ln(u1) sin( 2u2 )
In (9.23), u1 and u2 are two uniform random numbers generated independently by PRNG. They can be
transformed into a pair of standard normal random numbers 1 and 2 that are also independent. A
more efficient transformation is to rewrite (9.23) so as to avoid the costly evaluation of trigonometric
functions. This is known as the polar rejection algorithm given by
(9.24)
, if w ( 2u1 1 )
2 ( 2u2 1 )
2 1
The drawback in (9.24) is that the input uniform random numbers should stay inside the unit circle
around the origin.
In EXCEL spreadsheet, PRNG is given by the RAND() function that automatically reseeds based
on the system clock. The transformation of uniform random numbers into standard normal random
numbers are conducted by the NORMSINV(RAND()) function that takes RAND() as input. In VBA
implementation, it is highly inefficient to generate random numbers interactively through EXCEL as
the number of calls will be enormous in Monte Carlo simulation. It is therefore essential to develop a
VBA routine that implements for instance the polar rejection algorithm in (9.24). It can be defined
through an external function called StdNormNum( ) that returns a standard normal random number on
each request. The VBA code of this function is given by Code 9.1. In VBA, PRNG is given by the
Rnd() function in which the generated sequence is completely determined by the default seed value.
To prevent using repetitive sequences in Monte Carlo simulation, we can alter the seed based on the
system clock by calling Randomize once in the main routine. It is important to note that we should
not Randomize every time we call Rnd() as the random behavior can be badly skewed.
In both the Box-Muller and polar rejection algorithms, two usable standard normal random
numbers are generated at the same time. In order to obtain maximum efficiency, we should utilize
both of them by saving the unused random number for the next request. In Code 9.1, this can be done
through two static variables flagSave and snnSave that retain their latest values after termination of the
procedure. The static variable flagSave is initialized to 0 at the very first call to the function when
flagSave is still an empty cell checked by the IsEmpty function. When flagSave 0, the function
generates two random numbers and returns only the one under snnUse. The second random number is
saved under snnSave with flagSave set to 1. The two static variables will retain these values upon the
2 See Chapter 7 of W.H. Press, S.A. Teukolsky, W.T. Vetterling, and B.P. Flannery, Numerical Recipes in C :
The Art of Scientific Computing, 2nd Edition, Cambridge University Press, 1997.
1 ( 2u1 1 )
2ln(w)
w
2 ( 2u2 1 )
2ln(w)
w
6
next request where the function simply returns the random number under snnSave and resets flagSave
0.
__________________________________________________________________________________
Function StdNormNum() As Double
Dim v1 As Double, v2 As Double, w As Double, fac As Double
Dim snnUse As Double
Static flagSave As Integer: If IsEmpty(flagSave) Then flagSave = 0
Static snnSave As Double
If (flagSave = 0) Then
newtrial:
v1 = 2# * Rnd() – 1#
v2 = 2# * Rnd() – 1#
w = v1 ^ 2 + v2 ^ 2
If (w >= 1#) Then GoTo newtrial
fac = Sqr(-2# * Log(w) / w)
snnSave = fac * v1
snnUse = fac * v2
flagSave = 1
Else
snnUse = snnSave
flagSave = 0
End If
StdNormNum = snnUse
End Function
__________________________________________________________________________________
Code 9.1 : VBA code of the StdNormNum( ) function.
To generate a very long random sequence, it is advisable to use instead the standard algorithm for
PRNG known as ran0( ) proposed by Park and Miller2. This generator has passed all theoretical tests
so far and is definitely safer to use despite it is relatively slower than Rnd(). In Code 9.2 , we have
included the VBA code of ran0( ) and refer the reader to the original literature for detail discussion of
the algorithm. Like that in Rnd(), the seed value in ran0( ) should only be initiated once in the main
routine. It can be any nonzero integer such as seed 56789 and it must not be altered between
successive calls in a random sequence.
__________________________________________________________________________________
Public seed As Long
Function ran0() As Double
Dim IA As Long: IA = 16807
Dim IM As Long: IM = 2147483647
Dim IQ As Long: IQ = 127773
Dim IR As Long: IR = 2836
Dim MASK As Long: MASK = 123459876
Dim AM As Double: AM = 1# / IM
Dim k As Long
seed = seed Xor MASK
k = seed / IQ
seed = IA * (seed – k * IQ) – IR * k
If (seed < 0) Then seed = seed + IM
ran0 = AM * seed
seed = seed Xor MASK
End Function
__________________________________________________________________________________
Code 9.2 : VBA code of the ran0( ) function.
We consider again the pricing of a European call option written on a non-dividend paying stock as
discussed in previous section. Here, it is efficient and convenient to develop VBA routines that
separately perform the three different simulation tasks as demonstrated in Figure 9.2. The VBA codes
of the McEuropeanCall( ) routine for crude simulation, the CMcEuropeanCall( ) routine for control
variate method, and the AMcEuropeanCall( ) routine for antithetic variate method are given by Code
9.3. They all require the input parameters { S0 , K , r , , T , n } and return the estimation of the
current option price together with the standard error. In both routines, random sample of is generated
by calling the StdNormNum( ) function, and it will be used to generate sample maturity price ST()
7
using (9.16). For antithetic variate method, the same random number with a negative sign will also be
used to generate the antithetic price ST(). Sample value of the functions fT(), fT̃(), and fT̂() are then
evaluated respectively in different routines using (9.18), (9.19), and (9.21). The estimates of the mean
and variance of their present values can be obtained by summing up all sample values and their
squares as in (9.2) through a loop that runs from 1 to n. It is also convenient to implement an
alternative expression for the variance estimate in (9.2) as
(9.25)
In both routines, the mean estimate will be returned as the estimation of the current option price. For
control variate method, it should be reminded to include in the mean estimate the analytic solution of
S0 for the control variate as discussed in (9.20). The standard error in the estimation can be determined
by taking an additional factor of 1/√n in the variance estimate. Table 9.1 illustrates the relative
performance of the three routines in terms of the sizes of standard errors and computation times. We
adopt the same parameterization as in Figure 9.2 and consider one million samples in the Monte Carlo
simulation. In this example, the control variate method is shown to be highly effective, while there is
only mild improvement for the antithetic variate method.
Crude Simulation Control Variate Method Antithetic Variate Method
Mean 12.338 12.334 12.340
Standard Error 0.018 0.011 0.010
Relative CPU Time 1.00 1.02 1.33
Table 9.1 : Relative performance of different Monte Carlo routines for European call option.
As our second example, we consider a European call option written on a stock with expected
dividend payments { D1, … , DN } at { t1, … , tN T } prior to its maturity. Assume Sti represents the
asset price right after the dividend Di has been paid at ti. We can use (9.22) to first generate the asset
price at ti 1 using Sti and then incorporate the price drop due to the payment Di 1 as
Sti 1 Sti exp( ( r ½
2 )( ti 1 ti )
ti 1 ti i 1 ) Di 1 (9.26)
In this way, the maturity price ST of the asset can be generated iteratively by running i from 0 to N 1.
If at any time the R.H.S. of (9.26) is less than or equal zero, the iteration should stop and restart all
over again from i 0. The maturity price can be used to evaluate the option payoff according to the
function
fT( 1, … , N )
max{ ST( 1, … , N ) K , 0 } (9.27)
For variance reduction, it is effective to adopt the case of a non-dividend paying stock as the control
variate and consider the new function
fT̃( 1, … , N )
max{ ST( 1, … , N ) K , 0 }
max{ ST
(0)
( 1, … , N ) K , 0 } (9.28)
The control price ST
(0)
at maturity is generated simultaneously with ST using the same random sequence
but with zero dividend payments. The control variate has an analytic solution given by the Black-
Scholes formula in (9.17) for which (9.13) can be written as
f0 f0
BS
Ê( e
rT
fT̃ | S0 ) (9.29)
It is always straightforward to take as the antithetic variate and generate the antithetic value of ST
using the negative sequence { 1, … , N }. It is then convenient to introduce the new function
s2
n
i 1 g
2(xi) m
2
1
n 1
n
n 1
8
fT̂( 1, … , N ) ½ [
max{ ST( 1, … , N ) K , 0 }
max{ ST( 1, … , N ) K , 0 } ] (9.30)
as the two payoffs are negatively correlated. However, it has been shown below that the effectiveness
of the antithetic variate method in (9.30) is limited.
The pseudo code of the McEuropeanCallDiv( ) routine for crude simulation is given by Code 9.4.
The routine reads in the parameters { S0 , K , r , , n } as well as the details of dividend payments as
N, { t1, … , tN T }, and { D1, … , DN }. To initiate (9.26), the current time t0 0 by definition should
also be included in the input time array. Similarly, the routine returns the estimation of the current
option price together with the standard error. In the inner loop, sample maturity price ST of the asset is
generated iteratively by running i in (9.26) from 0 to N 1. The iteration starts off from S0 and moves
forward by calling the StdNormNum( ) function successively. It should restart from S0 at the
beginning of the loop if the asset price is less than or equal zero at any time. Sample payoff fT in
(9.27) and its present value can be evaluated using the asset price immediately after the iteration. The
estimates of the mean and variance can be obtained based on the sample sum and squared sum
accumulated through the outer loop that runs from sample number 1 to n. The VBA code of the
McEuropeanCallDiv( ) routine is given by Code 9.5.
It is easy to modify the algorithm for crude simulation and include the appropriate procedures for
variance reduction. The VBA code of the CMcEuropeanCallDiv( ) routine for the control variate
method is given by Code 9.6. In the inner loop, the iteration of the control price is performed
simultaneously with ST using the same random sequence and the same starting price but with no
dividend payment. Again, we should restart both iterations if either one of them is less than or equal
zero. In this way, sample value of the function fT̃ in (9.28) can be evaluated based on the iterated
prices. In the estimation of the option price, we have also included the Black-Scholes formula given
by the function BsEuropeanCall( ). The pseudo code of the AMcEuropeanCallDiv( ) routine for the
antithetic variate method is also given by Code 9.6. In this case, the antithetic price is generated
simultaneously with ST using however the negative sequence. Similarly, sample value of the function
fT̂ in (9.30) can be evaluated based on the prices immediately after the iteration. Table 9.2 illustrates
the relative performance of the three routines based on one million simulation samples. We adopt the
same parameterization as in Figure 9.2 and consider the dividend payments of D1,2,3,4 0.25 at
quarterly time intervals { t1 0.25 , t2 0.50 , t3 0.75 , t4 1.00 }. Again, the control variate method
is shown to be much more effective than the antithetic variate method.
Crude Simulation Control Variate Method Antithetic Variate Method
Mean 11.782 11.794 11.802
Standard Error 0.018 0.0005 0.010
Relative CPU Time 1.00 1.31 1.36
Table 9.2 : Relative performance of different Monte Carlo routines for European call option written
on dividend paying asset.
9.4. Exotic Options
In this section, we apply the Monte Carlo method in the pricing of exotic options with path-
dependent payoffs at maturity. For example, an Asian call option will have its payoff depending on
the average asset price over a set of predetermined forward time { t1, t2, … , tN T } given by
(9.31)
AT ( St1 … StN )
1
N
9
Again, we can generate iteratively the asset prices at every averaging time by running the index i from
0 to N 1 in (9.22) with random numbers { 1, … , N }. The average asset price can be used to
evaluate the option payoff according to the function
fT( 1, … , N )
max{ AT( 1, … , N ) K , 0 } (9.32)
For variance reduction, it is effective to adopt the case for the geometric average asset price as the
control variate and consider the new function
fT̃( 1, … , N )
max{ AT( 1, … , N ) K , 0 }
max{ GT( 1, … , N ) K , 0 } (9.33)
The geometric average asset price is defined as
GT
N
St1 … StN (9.34)
It is easy to show that GT is lognormally distributed with risk-neutral average given by
Ê( GT | S0 ) S0 e
r1 ½
2( 2 1
) , (9.35)
The control variate has an analytic solution given by the Black-Scholes formula as3
f0
BS
Ê( GT | S0 ) e
rT N( ) K e
rT
N(
2 ) , (9.36)
The pseudo code of the CMcArAsianCall( ) routine for the pricing of an Asian call option is given
by Code 9.7. The routine reads in the parameters { S0 , K , r , , n }, a set of forward averaging time {
t1, … , tN T }, and N. Similarly, the input time array should also include the current time defined as
t0 0 to initiate (9.22). Again, the routine returns the estimation of the current option price together
with the standard error. In the inner loop, asset prices at every averaging time are generated iteratively
for which their arithmetic and geometric averages can be evaluated. It should be noted that the
multiplication of a long sequence of prices could possibly overflow the numerical procedure. It is
safer to determine the geometric average by summing up the logarithmic prices. Sample payoff fT̃ in
(9.33) and its present value can be evaluated using the iterated average prices. The estimates of the
mean and variance can be obtained based on the sample sum and squared sum. We have also included
the Black-Scholes formula for the geometric case given by the function BsGeAsianCall( ). The VBA
code of the CMcArAsianCall( ) routine is given by Code 9.8. Table 9.3 demonstrates the effectiveness
of the use of geometric average price as control variate based on one million simulation samples. We
adopt the same parameterization as in Figure 9.2 and consider the set of averaging time { t1 0.50 , t2
0.75 , t3 1.00 }.
Crude Simulation Control Variate Method
Mean 9.691 9.6787
Standard Error 0.014 0.0004
Relative CPU Time 1.00 1.12
Table 9.3 : Relative performance of the CMcArAsianCall( ) routine as compared with
crude simulation.
As our second example on exotic option, we consider a double barrier knock out (DKO) European
call option with maturity payoff depending on the breaching condition of the asset price with respect
to two predefined barriers ( H S0 L ). In principle, the intermediate knock-out condition is
3 See A. Kemna and A. Vorst, “A pricing method for options based on average asset values”, Journal of
Banking and Finance, Vol. 14, Issue 1 (1990) p 113-129.
ln( Ê( GT | S0 )/K ) ½
22
2
N
i 1 ( N i 1 )
( ti ti 1 )
1
N
10
subjected to continuous asset price movement between time 0 and T. In (9.22), we can adopt equal
time interval of t T/N for large N such that the generated prices { St0 S0 , St1 , … , StN ST } at
discrete time ti it will be a good approximation of a continuous path. Thus, we can generate
iteratively the asset price movement for the entire life of the option by running i from 0 to N 1 as
Sti 1 Sti exp( ( r ½
2 )t
t i 1 ) (9.37)
with random numbers { 1, … , N }. If the asset price has hit either one of the barriers at any time
prior to the option maturity, the option is immediately knocked out and the payoff is set to zero or to a
constant amount of rebate c payable upon hitting. Otherwise, it will be evaluated according to the
maturity price based on a call-type payoff function as
c e
r( T th ) , if Sth( 1, … , h ) H or Sth( 1, … , h ) L for th T
fT( 1, … , N ) (9.38)
max{ ST( 1, … , N ) K , 0 } , otherwise
For variance reduction, it is difficult in this case to identify a control variate that works consistently
and effectively for a general barrier condition. We have also checked that the antithetic price
movement generated through the negative random sequence doesn’t work well for barrier option.
The pseudo code of the McDKOCall( ) routine for the crude pricing of a DKO European call
option is given by Code 9.9. The routine reads in the parameters { S0 , K , H , L , c , r , , N , T , n }
and returns the estimation of the option price together with the standard error. Here, the intermediate
knock-out condition is examined through an external procedure called DKOBoundary( ) that assigns
crossflag TRUE if the input price has surpassed either one of the barriers H and L. It is essential to
first examine the knock-out condition at current time and return the trivial solution f0 0 with no
estimation error. In the inner loop, sample asset price movement with time interval t is generated
iteratively using (9.37) starting from its current value S0. The knock-out condition is checked at every
time step by calling DKOBoundary( ) for the newly iterated price. If crossflag TRUE at any time,
the hitting time is tagged as th and the iteration should stop immediately. Otherwise, it should continue
through the entire loop and exit as the maturity price ST. In either case, sample payoff fT in (9.38) can
be evaluated according to the flag crossflag in conjunction with the hitting time or the maturity price.
The VBA code of the McDKOCall( ) routine is given by Code 9.10.
There are two sources of error in the Monte Carlo pricing of the DKO option. Firstly, the method
relies on statistical estimation and the statistical error is governed by the sample size. Secondly, the
statistical estimates are actually referring to the case with discrete time steps. Thus, there is systematic
error in the method as we are not estimating the true option price under continuous price movement.
The error would depend on the size of t in the configuration of the problem. It can be improved by
taking smaller time steps, but there is a tradeoff between computational time and accuracy. Table 9.4
demonstrates the significance of the size of t in the pricing based on one million simulation samples.
We adopt the same parameterization as in Figure 9.2 and consider the barrier levels of H 120 and L
80 with zero rebate c 0. It can be seen in the table that the effect is quite critical especially for the
mean estimate. We have shown in Figure 9.3 that the systematic error4 for the mean estimate will only
4 Consider the systematic error given by | mt mtrue | a(t)
b where a and b are constants. As mtrue is not known
a priori, we can choose mbest to be the best estimate in Table 9.4 and write
log10| mt mbest | b log10t log10a log10| 1 | , ( mbest mtrue )/( mt mtrue )
In the region where t is not too small, the term is negligible and b can be estimated by the slope of the plot
log10| mt mbest | versus log10t in Figure 9.3. It can be shown that b 0.5 in the linear fitting.
11
decrease like the square root of t. If the size of t is not small enough, the error will be quite
substantial.
t 0.1 t 0.01 t 0.001 t 0.0001
( N 10 ) ( N 100 ) ( N 1000 ) ( N 10000 )
Mean 1.2328 0.7533 0.6089 0.5661
Standard Error 0.0035 0.0027 0.0024 0.0023
Relative CPU Time 1.00 7.51 69.61 697.22
Table 9.4 : Relative performance of the McDKOCall( ) routine under different configurations of time
step.
Figure 9.3 : The plot of systematic error versus t for the pricing of the DKO call option.
9.5. American Options
One of the most important problems in option pricing is the valuation of American-style options
with early exercising features. For example, an American put option can be exercised at any time prior
to its maturity with intrinsic value evaluated according to the payoff function (S) max{ K S , 0 }.
At time t with underlying asset price St, the fair value of an American put option can be defined based
on the risk-neutral expectation of its discounted payoff as
F( St , t ) Ê( e
r( s t ) (Ss) | St ) (9.39)
where s is the first exercising time along a random asset price movement starting off from St. It can
occur at option’s maturity or at any earlier time when the intrinsic value is greater than the fair value
at that point. Similarly, it is optimal to early exercise the option at time t based on the criteria as
(St) F( St , t ) (9.40)
The price of the American put option should also include this feature and is given by
f( St , t ) max{ F( St , t ) , (St) } (9.41)
t
| mt mbest |
Slope 0.54
12
The equality condition in (9.40) can serve to define a critical price Sc(t) at time t for which the early
exercising criteria will be satisfied5 whenever St Sc(t). In this way, the criteria can be implemented
very easily by checking the asset price with respect to a critical price. It should be noted that the fair
value F( St , t ) can only be evaluated when all critical prices, or a boundary, at forward time of t have
been determined. This means that critical price Sc(t) can only be determined through backward
iteration starting off from the option’s maturity with Sc(T) K. To determine the current price of the
American put option, we need to first generate the entire critical boundary prior to its maturity and
then evaluate the risk-neutral expectation as
f0 max{ Ê( e
rs (Ss) | S0 ) , (S0) } (9.42)
Again, we can adopt equal time interval of t T/N for large N and consider the discrete time steps
ti it for the entire life of the option with i runs from 0 to N. At time ti , consider Sti x and suppose
the critical boundary at forward time is given by { Sc(ti 1) , … , Sc(tN) K }. The fair value of the
American put option can be estimated by Monte Carlo simulation as
F( Sti x , ti ) Ê( e
r( s ti ) (Ss) | Sti x ) (9.43)
In (9.43), we can generate iteratively the asset price movement for the remaining life of the option by
running j in (9.44) from i to N 1 starting from Sti x.
Stj 1 Stj exp( ( r ½
2 )t
t j ) (9.44)
If the generated price Stj 1 Sc(tj 1) at any forward time prior to the option’s maturity, the option
should be exercised immediately and s tj 1 in (9.43). Otherwise, it will be exercised at maturity for
positive payoff and s T. The critical price Sc(ti) at time ti can be determined by identifying the value
of x for which the mean estimate in (9.43) matches the intrinsic value. This is the same as solving the
root of the function given by
y(x) F( x , ti ) (x) (9.45)
where y(x) is subjected to the standard error in the estimation of F( x , ti ). It can be done by
performing a linear regression for y(x) x in the neighborhood6 of its trial value x Sc(ti 1) with
xroot /. It is essential for the size of y(xroot) to be less than the standard error in F( x , ti ).
Otherwise, a new xroot should be determined through another linear regression around the old one until
y(xroot) becomes zero within the statistical error. In this way, the entire critical boundary can be
determined by iterating the above procedure for all ti in a backward scheme starting from tN 1 to t0
with Sc(tN) K.
We first develop a routine called FvAmericanPut( ) capable of estimating the fair value of an
American put option at forward time based on Monte Carlo simulation. The pseudo code of
FvAmericanPut( ) is given by Code 9.11. It reads in the parameters { x , K , r , , N , T , n }, a time
pointer i for the reference forward time ti , and the critical prices { Sc(ti 1) , … , Sc(tN) } thereafter. The
routine returns the Monte Carlo estimation of the fair value together with the standard error. In the
inner loop, asset price movement between ti and T is generated iteratively using (9.44) with time
interval of t. At each time step, the early exercising criteria is examined through an external
procedure called APBoundary( ) that assigns exerciseflag TRUE if the iterated asset price is less
than the critical price. If exerciseflag TRUE at any time, the iteration should terminate immediately
with stopping time taken to be the updated s in the loop. Otherwise, it should continue until the
5 This is true as both (St) and F( St , t ) are decreasing function with St. Also, (St) vanishes for large St.
6 Avoid taking the out-of-money region of the option by virtue of the linear approximation for (9.45).
13
option matures at which s T. The discounted payoff in (9.43) can then be evaluated using the
stopping time and the terminal asset price in the iteration. The VBA code of FvAmericanPut( ) is
given by Code 9.12.
We then develop a routine called McAmericanPut( ) that performs the backward iteration based on
the FvAmericanPut( ) routine. The pseudo code of McAmericanPut( ) is given by Code 9.13. It reads
in the parameters { S0 , K , r , , N , T , n } and returns the estimation of the option price together
with the standard error. In the backward iteration, we can generate the critical boundary by running
the time pointer i in FvAmericanPut( ) backward from N 1 to 0 with Sc(tN) K being the starting
value. At time ti in the loop, the critical prices from ti 1 to tN are presumably calculated in previous
steps. It is then straightforward to estimate F( x , ti ) in (9.45) by calling FvAmericanPut( ) and the
main task is to determine the root of y(x) that corresponds to the critical price Sc(ti) at time ti. We
adopt Sc(ti 1) as the trial value of the root and consider a linear regression for y(x) in the
neighborhood of this value. In the inner loop, we generate h 100 data points for the regression with
values of x that spread across a small region of 1% below the trial value. The linear regression can
be performed through an external routine called Regn( ) and it should be noted that the errors in y(x)
are roughly constant7. The resulting and would provide a good approximation for xroot and the
corresponding y(xroot) should be calculated by calling FvAmericanPut( ). In the IF statement, it is
important to check explicitly that y(xroot) is zero within the standard error in the estimation. If so, we
can assign xroot to be the critical price Sc(ti). Otherwise, it can be used to update the trial value and the
regression should be repeated starting from the line labeled as “NewTrialValue”. The VBA code of
McAmericanPut( ) is given by Code 9.14.
Once we have determined the critical boundary in the backward iteration, the current price of the
American put option can be estimated by calling FvAmericanPut( ) with S0 and compares the fair
value with the intrinsic value according to (9.42). The systematic error here would depend on the
accuracy of the critical boundary determined under discrete time interval. Figure 9.4 illustrates the
deviation in the critical boundaries determined based on, for example, t 0.1 and t 0.001. Table
9.5 demonstrates the overall significance of the size of t in the pricing of the American put option.
We adopt the same parameterization as in Figure 9.2 and consider again one million samples in the
simulation. It has been show in Figure 9.5 that the systematic error for the mean estimate will only
decrease like the linear order of t. Thus, the error will be less substantial here than in the pricing of
DKO option.
t 0.1 t 0.01 t 0.001
( N 10 ) ( N 100 ) ( N 1000 )
Mean 7.909 7.982 7.985
Standard Error 0.009 0.009 0.009
Relative CPU Time 1.00 26.76 752.62
Table 9.5 : Relative performance of the McAmericanPut( ) routine under different configurations
of time step.
7 In the fitting of a set of h data points (xk , yk) to a straight line y x where y is subjected to a constant
measurement error, the best fit parameters in the linear model are given by
( xy xy ) / , ( x2y xxy ) / , x2 x2
The estimation errors of and are both in the order of 1 /√h. Thus, the size of h must be sufficiently large to
ensure stability in the iteration of the critical boundary.
14
Figure 9.4 : Critical boundary of the American put option.
Figure 9.5 : The plot of systematic error versus t for the pricing of the American put option.
t
t 0.1
t 0.001
t
| mt mbest |
slope 1.08
15
__________________________________________________________________________________
Sub McEuropeanCall(assetPrice As Double, strike As Double, riskFree As Double, sigma As Double, maturity As Double,
nsample As Long, ByRef optionPrice As Double, ByRef stdErr As Double)
Dim sum As Double, sum2 As Double, gen As Double
Dim mean As Double, sd As Double, Ls As Long
Dim St As Double, fT As Double, pV As Double
sum = 0
sum2 = 0
For Ls = 1 To nsample
gen = StdNormNum()
St = assetPrice * Exp((riskFree - 0.5 * sigma ^ 2) * maturity + sigma * Sqr(maturity) * gen)
fT = CallPayoff(strike, St)
pV = Exp(-riskFree * maturity) * fT
sum = sum + pV
sum2 = sum2 + pV * pV
Next Ls
mean = sum / nsample
sd = Sqr(sum2 / (nsample - 1) - (nsample / (nsample - 1)) * mean ^ 2)
optionPrice = mean
stdErr = sd / Sqr(nsample)
End Sub
__________________________________________________________________________________
Sub CMcEuropeanCall(assetPrice As Double, strike As Double, riskFree As Double, sigma As Double, maturity As Double,
nsample As Long, ByRef optionPrice As Double, ByRef stdErr As Double)
Dim sum As Double, sum2 As Double, gen As Double
Dim mean As Double, sd As Double, Ls As Long
Dim St As Double, fT As Double, pV As Double
sum = 0
sum2 = 0
For Ls = 1 To nsample
gen = StdNormNum()
St = assetPrice * Exp((riskFree - 0.5 * sigma ^ 2) * maturity + sigma * Sqr(maturity) * gen)
fT = CallPayoff(strike, St) - St
pV = Exp(-riskFree * maturity) * fT
sum = sum + pV
sum2 = sum2 + pV * pV
Next Ls
mean = sum / nsample
sd = Sqr(sum2 / (nsample - 1) - (nsample / (nsample - 1)) * mean ^ 2)
optionPrice = assetPrice + mean
stdErr = sd / Sqr(nsample)
End Sub
__________________________________________________________________________________
Sub AMcEuropeanCall(assetPrice As Double, strike As Double, riskFree As Double, sigma As Double, maturity As Double,
nsample As Long, ByRef optionPrice As Double, ByRef stdErr As Double)
Dim sum As Double, sum2 As Double, gen As Double
Dim mean As Double, sd As Double, Ls As Long
Dim St As Double, Sta As Double, fT As Double, pV As Double
sum = 0
sum2 = 0
For Ls = 1 To nsample
gen = StdNormNum()
St = assetPrice * Exp((riskFree - 0.5 * sigma ^ 2) * maturity + sigma * Sqr(maturity) * gen)
Sta = assetPrice * Exp((riskFree - 0.5 * sigma ^ 2) * maturity + sigma * Sqr(maturity) * (-gen))
fT = (CallPayoff(strike, St) + CallPayoff(strike, Sta)) / 2
pV = Exp(-riskFree * maturity) * fT
sum = sum + pV
sum2 = sum2 + pV * pV
Next Ls
mean = sum / nsample
sd = Sqr(sum2 / (nsample - 1) - (nsample / (nsample - 1)) * mean ^ 2)
optionPrice = mean
stdErr = sd / Sqr(nsample)
End Sub
__________________________________________________________________________________
16
Function CallPayoff(strike As Double, assetPrice As Double) As Double
CallPayoff = Max(assetPrice - strike, 0)
End Function
Function Max(x As Double, y As Double) As Double
If x > y Then Max = x Else Max = y
End Function
__________________________________________________________________________________
Code 9.3 : VBA codes of the McEuropeanCall( ), CMcEuropeanCall( ), and AMcEuropeanCall( )
routines.
17
__________________________________________________________________________________
McEuropeanCallDiv( S0 , K , r , , N , t(0 : N) , D(1 : N) , n , f0 , error )
# zeroize the sample sum and squared sum
sum 0 , sum2 0
For( Ls 1 to n ){
# initialize the asset price
St S0
# generate the asset price at each dividend date
For( i 0 to N 1 ){
StdNormNum( )
St St exp( ( r ½2 )[ t( i 1 ) t( i ) ]
t( i 1 ) t( i ) ) D( i 1 )
If( St 0 ) go back to the statement St S0 and restart all over again }
# evaluate the option payoff function as in (9.27)
fT CallPayoff( K , St )
PV e
rT
fT
# accumulate the sample sum and squared sum
sum sum PV
sum2 sum2 PV 2 }
# evaluate the estimates of mean and variance
m sum / n
# return the estimation of option price and standard error
f0 m
error s /√n
__________________________________________________________________________________
Code 9.4 : Pseudo code of the McEuropeanCallDiv( ) routine.
s sum2 m
2 1
n 1
n
n 1
18
__________________________________________________________________________________
Sub McEuropeanCallDiv(assetPrice As Double, strike As Double, riskFree As Double, sigma As Double, nDiv As Integer,
timeDiv() As Double, paymentDiv() As Double, nsample As Long, ByRef optionPrice As Double,
ByRef stdErr As Double)
Dim sum As Double, sum2 As Double, gen As Double
Dim mean As Double, sd As Double, Ls As Long
Dim St As Double, fT As Double, pV As Double, i As Integer
sum = 0
sum2 = 0
For Ls = 1 To nsample
NewTrial:
St = assetPrice
For i = 0 To nDiv – 1
gen = StdNormNum()
St = St * Exp((riskFree – 0.5 * sigma ^ 2) * (timeDiv(i + 1) – timeDiv(i)) + sigma * Sqr(timeDiv(i + 1) – timeDiv(i)) * gen) – paymentDiv(i + 1)
If (St <= 0) Then GoTo NewTrial
Next i
fT = CallPayoff(strike, St)
pV = Exp(-riskFree * timeDiv(nDiv)) * fT
sum = sum + pV
sum2 = sum2 + pV * pV
Next Ls
mean = sum / nsample
sd = Sqr(sum2 / (nsample - 1) - (nsample / (nsample - 1)) * mean ^ 2)
optionPrice = mean
stdErr = sd / Sqr(nsample)
End Sub
__________________________________________________________________________________
Code 9.5 : VBA code of the McEuropeanCallDiv( ) routine.
19
__________________________________________________________________________________
Sub CMcEuropeanCallDiv(assetPrice As Double, strike As Double, riskFree As Double, sigma As Double, nDiv As Integer,
timeDiv() As Double, paymentDiv() As Double, nsample As Long, ByRef optionPrice As Double,
ByRef stdErr As Double)
Dim sum As Double, sum2 As Double, gen As Double
Dim mean As Double, sd As Double, Ls As Long
Dim St As Double, St0 As Double, fT As Double, pV As Double, i As Integer
sum = 0
sum2 = 0
For Ls = 1 To nsample
NewTrial:
St = assetPrice
St0 = assetPrice
For i = 0 To nDiv - 1
gen = StdNormNum()
St = St * Exp((riskFree - 0.5 * sigma ^ 2) * (timeDiv(i + 1) - timeDiv(i)) + sigma * Sqr(timeDiv(i + 1) - timeDiv(i)) * gen) - paymentDiv(i + 1)
St0 = St0 * Exp((riskFree - 0.5 * sigma ^ 2) * (timeDiv(i + 1) - timeDiv(i)) + sigma * Sqr(timeDiv(i + 1) - timeDiv(i)) * gen)
If (St <= 0 Or St0 <= 0) Then GoTo NewTrial
Next i
fT = CallPayoff(strike, St) - CallPayoff(strike, St0)
pV = Exp(-riskFree * timeDiv(nDiv)) * fT
sum = sum + pV
sum2 = sum2 + pV * pV
Next Ls
mean = sum / nsample
sd = Sqr(sum2 / (nsample - 1) - (nsample / (nsample - 1)) * mean ^ 2)
optionPrice = BsEuropeanCall(assetPrice, strike, riskFree, sigma, timeDiv(nDiv)) + mean
stdErr = sd / Sqr(nsample)
End Sub
__________________________________________________________________________________
Sub AMcEuropeanCallDiv(assetPrice As Double, strike As Double, riskFree As Double, sigma As Double, nDiv As Integer,
timeDiv() As Double, paymentDiv() As Double, nsample As Long, ByRef optionPrice As Double,
ByRef stdErr As Double)
Dim sum As Double, sum2 As Double, gen As Double
Dim mean As Double, sd As Double, Ls As Long
Dim St As Double, Sta As Double, fT As Double, pV As Double, i As Integer
sum = 0
sum2 = 0
For Ls = 1 To nsample
NewTrial:
St = assetPrice
Sta = assetPrice
For i = 0 To nDiv - 1
gen = StdNormNum()
St = St * Exp((riskFree - 0.5 * sigma ^ 2) * (timeDiv(i + 1) - timeDiv(i)) + sigma * Sqr(timeDiv(i + 1) - timeDiv(i)) * gen) - paymentDiv(i + 1)
Sta = Sta * Exp((riskFree - 0.5 * sigma ^ 2) * (timeDiv(i + 1) - timeDiv(i)) + sigma * Sqr(timeDiv(i + 1) - timeDiv(i)) * (-gen)) - paymentDiv(i + 1)
If (St <= 0 Or Sta <= 0) Then GoTo NewTrial
Next i
fT = (CallPayoff(strike, St) + CallPayoff(strike, Sta)) / 2
pV = Exp(-riskFree * timeDiv(nDiv)) * fT
sum = sum + pV
sum2 = sum2 + pV * pV
Next Ls
mean = sum / nsample
sd = Sqr(sum2 / (nsample - 1) - (nsample / (nsample - 1)) * mean ^ 2)
optionPrice = mean
stdErr = sd / Sqr(nsample)
End Sub
__________________________________________________________________________________
Function BsEuropeanCall(assetPrice As Double, strike As Double, riskFree As Double, sigma As Double, maturity As Double) As Double
Dim d1 As Double, d2 As Double
d1 = (Log(assetPrice / strike) + (riskFree + 0.5 * sigma ^ 2) * maturity) / (sigma * Sqr(maturity))
d2 = d1 - sigma * Sqr(maturity)
With Application.WorksheetFunction
BsEuropeanCall = assetPrice * .NormSDist(d1) - strike * Exp(-riskFree * maturity) * .NormSDist(d2)
End With
End Function
__________________________________________________________________________________
Code 9.6 : VBA codes of the CMcEuropeanCallDiv( ), and AMcEuropeanCallDiv( ) routines.
20
__________________________________________________________________________________
CMcArAsianCall( S0 , K , r , , N , t(0 : N) , n , f0 , error )
# zeroize the sample sum and squared sum
sum 0 , sum2 0
For( Ls 1 to n ){
# initialize the asset price, sum of price and logarithmic sum of price
St S0 , S 0 , lnS 0
# generate the asset price at each averaging time
For( i 0 to N 1 ){ StdNormNum( )
St St exp( ( r ½2 )[ t( i 1 ) t( i ) ]
t( i 1 ) t( i ) )
S S St
lnS lnS ln( St ) }
# evaluate the arithmetic and geometric average asset prices
AT S / N
GT exp( lnS / N )
# evaluate the option payoff function as in (9.33)
fT CallPayoff( K , AT ) CallPayoff( K , GT )
PV e
rT
fT
# accumulate the sample sum and squared sum
sum sum PV
sum2 sum2 PV 2 }
# evaluate the estimates of mean and variance
m sum / n
# return the estimation of option price and standard error
f0 BsGeAsianCall( S0 , K , r , , N , t(0 : N) ) m
error s /√n
__________________________________________________________________________________
Code 9.7 : Pseudo code of the CMcArAsianCall( ) routine.
s sum2 m
2 1
n 1
n
n 1
21
__________________________________________________________________________________
Sub CMcArAsianCall(assetPrice As Double, strike As Double, riskFree As Double, sigma As Double, nAvg As Integer,
timeAvg() As Double, nsample As Long, ByRef optionPrice As Double, ByRef stdErr As Double)
Dim sum As Double, sum2 As Double, gen As Double
Dim mean As Double, sd As Double, Ls As Long
Dim St As Double, fT As Double, pV As Double, i As Integer
Dim sumPrice As Double, sumLnPrice As Double
Dim amAvg As Double, gmAvg As Double
sum = 0
sum2 = 0
For Ls = 1 To nsample
St = assetPrice
sumPrice = 0
sumLnPrice = 0
For i = 0 To nAvg - 1
gen = StdNormNum()
St = St * Exp((riskFree - 0.5 * sigma ^ 2) * (timeAvg(i + 1) - timeAvg(i)) + sigma * Sqr(timeAvg(i + 1) - timeAvg(i)) * gen)
sumPrice = sumPrice + St
sumLnPrice = sumLnPrice + Log(St)
Next i
amAvg = sumPrice / nAvg
gmAvg = Exp(sumLnPrice / nAvg)
fT = CallPayoff(strike, amAvg) - CallPayoff(strike, gmAvg)
pV = Exp(-riskFree * timeAvg(nAvg)) * fT
sum = sum + pV
sum2 = sum2 + pV * pV
Next Ls
mean = sum / nsample
sd = Sqr(sum2 / (nsample - 1) - (nsample / (nsample - 1)) * mean ^ 2)
optionPrice = BsGeAsianCall(assetPrice, strike, riskFree, sigma, nAvg, timeAvg) + mean
stdErr = sd / Sqr(nsample)
End Sub
__________________________________________________________________________________
Function BsGeAsianCall(assetPrice As Double, strike As Double, riskFree As Double, sigma As Double, nAvg As Integer,
timeAvg() As Double) As Double
Dim meanG As Double
Dim nu1 As Double, nu2 As Double, i As Integer
Dim d1 As Double, d2 As Double
nu1 = 0
nu2 = 0
For i = 1 To nAvg
nu1 = nu1 + (nAvg - i + 1) * (timeAvg(i) - timeAvg(i - 1))
nu2 = nu2 + (nAvg - i + 1) ^ 2 * (timeAvg(i) - timeAvg(i - 1))
Next i
nu1 = nu1 / nAvg
nu2 = nu2 / nAvg ^ 2
meanG = assetPrice * Exp(riskFree * nu1 + 0.5 * sigma ^ 2 * (nu2 - nu1))
d1 = (Log(meanG / strike) + 0.5 * sigma ^ 2 * nu2) / (sigma * Sqr(nu2))
d2 = d1 - sigma * Sqr(nu2)
With Application.WorksheetFunction
BsGeAsianCall = Exp(-riskFree * timeAvg(nAvg)) * (meanG * .NormSDist(d1) - strike * .NormSDist(d2))
End With
End Function
__________________________________________________________________________________
Code 9.8 : VBA code of the CMcAsianCall( ) routine.
22
__________________________________________________________________________________
McDKOCall( S0 , K , H , L, c , r , , N , T , n , f0 , error )
# check the knock-out condition for the current asset price
Call DKOBoundary( H , L , S0 , crossflag )
If( crossflag ) then
f0 0 , error 0
Exit subroutine
Endif
# define the size of time interval
t T / N
# zeroize the sample sum and squared sum
sum 0 , sum2 0
For( Ls 1 to n ){
# initialize the asset price
St S0
# generate the asset price at each intermediate time and check the knock-out condition
For( i 0 to N 1 ){ StdNormNum( )
St St exp( ( r ½2 ) t
t )
Call DKOBoundary( H , L , St , crossflag )
If( crossflag ) then
th it
exit i
Endif }
# evaluate the option payoff function as in (9.38)
If( crossflag ) then
fT c e
r( T th )
Else
fT CallPayoff( K , St )
Endif
PV e
rT
fT
# accumulate the sample sum and squared sum
sum sum PV
sum2 sum2 PV 2 }
# evaluate the estimates of mean and variance
m sum / n
# return the estimation of option price and standard error
s sum2 m
2 1
n 1
n
n 1
23
f0 m
error s /√n
__________________________________________________________________________________
Code 9.9 : Pseudo code of the McDKOCall( ) routine.
__________________________________________________________________________________
Sub McDKOCall(assetPrice As Double, strike As Double, upperBarrier As Double, lowerBarrier As Double, rebate As Double,
riskFree As Double, sigma As Double, nStep As Integer, maturity As Double, nsample As Long,
ByRef optionPrice As Double, ByRef stdErr As Double)
Dim sum As Double, sum2 As Double, gen As Double
Dim mean As Double, sd As Double, Ls As Long
Dim St As Double, fT As Double, pV As Double, i As Integer
Dim crossFlag As Boolean
Dim timeHit As Double
Dim dtime As Double: dtime = maturity / nStep
Call DKOBoundary(upperBarrier, lowerBarrier, assetPrice, crossFlag)
If (crossFlag) Then
optionPrice = 0
stdErr = 0
Exit Sub
End If
sum = 0
sum2 = 0
For Ls = 1 To nsample
St = assetPrice
For i = 0 To nStep - 1
gen = StdNormNum()
St = St * Exp((riskFree - 0.5 * sigma ^ 2) * dtime + sigma * Sqr(dtime) * gen)
Call DKOBoundary(upperBarrier, lowerBarrier, St, crossFlag)
If (crossFlag) Then
timeHit = i * dtime
Exit For
End If
Next i
If (crossFlag) Then
fT = rebate * Exp(riskFree * (maturity - timeHit))
Else
fT = CallPayoff(strike, St)
End If
pV = Exp(-riskFree * maturity) * fT
sum = sum + pV
sum2 = sum2 + pV * pV
Next Ls
mean = sum / nsample
sd = Sqr(sum2 / (nsample - 1) - (nsample / (nsample - 1)) * mean ^ 2)
optionPrice = mean
stdErr = sd / Sqr(nsample)
End Sub
__________________________________________________________________________________
Sub DKOBoundary(upperBarrier As Double, lowerBarrier As Double, assetPrice As Double, ByRef crossFlag As Boolean)
If (assetPrice >= upperBarrier Or assetPrice <= lowerBarrier) Then
crossFlag = True
Else
crossFlag = False
End If
End Sub
__________________________________________________________________________________
Code 9.10 : VBA codes of the McDKOCall( ) and DKOBoundary( ) routines.
24
__________________________________________________________________________________
FvAmericanPut( x , K , r , , N , T , i , Sc( i 1 : N ) , n , F , error )
# define the size of time interval and the reference forward time
t T / N , ti it
# zeroize the sample sum and squared sum
sum 0 , sum2 0
For( Ls 1 to n ){
# initialize the asset price
St x
# generate the asset price at each time step, check early exercising criteria, and update the latest time
For( j i to N 1 ){ StdNormNum( )
St St exp( ( r ½2 ) t
t )
Call APBoundary( St , Sc( j 1 ) , exerciseflag )
s ( j 1 )t
If( exerciseflag ) exit j }
# evaluate the discounted payoff in (9.43)
PV e
r( s ti ) PutPayoff( K , St )
# accumulate the sample sum and squared sum
sum sum PV
sum2 sum2 PV 2 }
# evaluate the estimates of mean and variance
m sum / n
# return the estimation of option price and standard error
F m
error s /√n
__________________________________________________________________________________
Code 9.11 : Pseudo code of the FvAmericanPut( ) routine.
s sum2 m
2 1
n 1
n
n 1
25
__________________________________________________________________________________
Sub FvAmericanPut(assetPrice As Double, strike As Double, riskFree As Double, sigma As Double, nStep As Integer, maturity As Double,
iptr As Integer, Scritical() As Double, nsample As Long, ByRef fairValue As Double, ByRef stdErr As Double)
Dim sum As Double, sum2 As Double, gen As Double
Dim mean As Double, sd As Double, Ls As Long
Dim St As Double, pV As Double, j As Integer
Dim exerciseFlag As Boolean
Dim timeStop As Double
Dim dtime As Double: dtime = maturity / nStep
Dim timeRef As Double: timeRef = iptr * dtime
sum = 0
sum2 = 0
For Ls = 1 To nsample
St = assetPrice
For j = iptr To nStep - 1
gen = StdNormNum()
St = St * Exp((riskFree - 0.5 * sigma ^ 2) * dtime + sigma * Sqr(dtime) * gen)
timeStop = (j + 1) * dtime
Call APBoundary(St, Scritical(j + 1), exerciseFlag)
If (exerciseFlag) Then Exit For
Next j
pV = Exp(-riskFree * (timeStop - timeRef)) * PutPayoff(strike, St)
sum = sum + pV
sum2 = sum2 + pV * pV
Next Ls
mean = sum / nsample
sd = Sqr(sum2 / (nsample - 1) - (nsample / (nsample - 1)) * mean ^ 2)
fairValue = mean
stdErr = sd / Sqr(nsample)
End Sub
__________________________________________________________________________________
Sub APBoundary(assetPrice As Double, criticalPrice As Double, exerciseFlag As Boolean)
If (assetPrice <= criticalPrice) Then
exerciseFlag = True
Else
exerciseFlag = False
End If
End Sub
__________________________________________________________________________________
Function PutPayoff(strike As Double, assetPrice As Double) As Double
PutPayoff = Max(strike - assetPrice, 0)
End Function
__________________________________________________________________________________
Code 9.12 : VBA code of the FvAmericanPut( ) routine.
26
__________________________________________________________________________________
McAmericanPut( S0 , K , r , , N , T , n , f0 , error )
# define the configuration of the linear regression
h 100 , 0.01
# define the critical price at option’s maturity
Sc( N ) K
# generate the critical boundary
For( i N 1 to 0, 1 ){
# define the trial value for the root of (9.45)
xroot Sc( i 1 )
NewTrialValue : inc xroot / h
# generate plotting points for (9.45) in the neighborhood of the trial value
For( k 1 to h ){ Xfit(k) xroot k inc
Call FvAmericanPut( Xfit(k) , K , r , , N , T , i , Sc( i 1 : N ) , n , F , error )
Yfit(k) F PutPayoff( K , Xfit(k) ) }
# perform a linear regression for (9.45) and obtain an approximation for its root
Call Regn( Xfit( 1 : h ) , Yfit( 1 : h ) , h , , )
xroot /
# check the validity of the approximated root and, if necessary, use it as a new trial value.
Call FvAmericanPut( xroot , K , r , , N , T , i , Sc( i 1 : N ) , n , F , error )
If( | F PutPayoff( K , xroot ) | error ) Then GoTo NewTrialValue
# assign the valid root of (9.45) as the critical price
Sc( i ) xroot }
# knowing the critical boundary, evaluate the option price and standard error according to (9.42)
Call FvAmericanPut( S0 , K , r , , N , T , i , Sc( 0 : N ) , n , F , error )
f0 max( F , PutPayoff( K , S0 ) )
__________________________________________________________________________________
Code 9.13 : Pseudo code of the McAmericanPut( ) routine.
27
__________________________________________________________________________________
Sub McAmericanPut(assetPrice As Double, strike As Double, riskFree As Double, sigma As Double, nStep As Integer, maturity As Double,
nsample As Long, ByRef optionPrice As Double, ByRef stdErr As Double)
Dim Scritical() As Double: ReDim Scritical(0 To nStep)
Dim iptr As Integer, k As Integer
Dim nFit As Integer: nFit = 100
Dim delta As Double: delta = 0.01
Dim inc As Double
Dim xroot As Double, fairValue As Double
Dim xFit() As Double: ReDim xFit(1 To nFit)
Dim yFit() As Double: ReDim yFit(1 To nFit)
Dim slope As Double, intercept As Double
Scritical(nStep) = strike
For iptr = nStep - 1 To 0 Step -1
xroot = Scritical(iptr + 1)
NewTrialValue: inc = delta * xroot / nFit
For k = 1 To nFit
xFit(k) = xroot - k * inc
Call FvAmericanPut(xFit(k), strike, riskFree, sigma, nStep, maturity, iptr, Scritical, nsample, fairValue, stdErr)
yFit(k) = fairValue - PutPayoff(strike, xFit(k))
Next k
Call Regn(xFit, yFit, nFit, slope, intercept)
xroot = -intercept / slope
Call FvAmericanPut(xroot, strike, riskFree, sigma, nStep, maturity, iptr, Scritical, nsample, fairValue, stdErr)
If (Abs(fairValue - PutPayoff(strike, xroot)) > stdErr) Then GoTo NewTrialValue
Scritical(iptr) = xroot
Range(“A1”).Offset(nStep – 1 – iptr, 0) = iptr * (maturity / nStep)
Range(“B1”).Offset(nStep – 1 – iptr, 0) = xroot
Application.StatusBar = “Done simulation i : ” & iptr
Next iptr
Call FvAmericanPut(assetPrice, strike, riskFree, sigma, nStep, maturity, 0, Scritical, nsample, fairValue, stdErr)
optionPrice = Max(fairValue, PutPayoff(strike, assetPrice))
End Sub
__________________________________________________________________________________
Sub Regn(xFit() As Double, yFit() As Double, nFit As Integer, ByRef slope As Double, ByRef intercept As Double)
Dim avgX As Double, avgY As Double, avgX2 As Double, avgY2 As Double, avgXY As Double
Dim i As Integer
avgX = 0
avgY = 0
avgX2 = 0
avgY2 = 0
avgXY = 0
For i = 1 To nFit
avgX = avgX + xFit(i)
avgY = avgY + yFit(i)
avgX2 = avgX2 + xFit(i) ^ 2
avgY2 = avgY2 + yFit(i) ^ 2
avgXY = avgXY + xFit(i) * yFit(i)
Next i
avgX = avgX / nFit
avgY = avgY / nFit
avgX2 = avgX2 / nFit
avgY2 = avgY2 / nFit
avgXY = avgXY / nFit
slope = (avgXY – avgX * avgY) / (avgX2 – avgX ^ 2)
intercept = (avgX2 * avgY – avgX * avgXY) / (avgX2 – avgX ^ 2)
End Sub
__________________________________________________________________________________
Code 9.14 : VBA code of the McAmericanPut( ) routine.