程序代写代做代考 scheme Excel flex algorithm finance Chapter 1

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

√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( 2u2 ) (9.23)

2  

2ln(u1) sin( 2u2 )

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 r1  ½ 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 )  ½ 22   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  it 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 | mt  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| mt  mbest |  b log10t  log10a  log10| 1   | ,   ( mbest  mtrue )/( mt  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| mt  mbest | versus log10t 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 ) (Ss) | 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 | mt  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 rs (Ss) | 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  it 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 ) (Ss) | 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  xy ) /  ,   ( x2y  xxy ) /  ,   x2  x2 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 | mt  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  it 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  it # 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.