gives 4.212.)
4 As mentioned in Section 12.13, o* is in theory slightly greater than o, the volatility of S. In practice, the use ofimplied volatilities avoids the need for analysts to distinguish between o and o*.Numerical Procedures 405
At each node:
Upper vaiue = Underlying Asset PriceLower value = Option Price
Shading indicates where option is exercised
Strike price = 50Discount factor per step = 0.9917
Time step, dt = 0.0833 years, 30.42 daysGrowth factor per step, a = 1.0084Probability of up move, p = 0.5073Up step size, u = 1.1224Down step size, d = 0.8909
         
 
 
 
Node Time:0.0000 0.0833 0.1667 0.2500 0.3333 0.4167
 
18.4 EXTENSIONS OF THE BASIC TREE APPROACH
 
 
We now explain some ways the basic tree methodology we have described can be extended.
Time-Dependent Interest Rates
The usual assumption when American options are being valued is that interest rates are constant.When the term structure is steeply upward or downward sloping, this may not be a satisfactoryassumption. It is more appropriate to assume that the interest rate for a future period of length dtequals the current forward interest rate for that period. We can accommodate this assumption bysetting
a= ef (18.12)406 CHAPTER 18
for nodes at time t, where f(t) is the forward rate between times t and ¢ + d¢. This does not changethe geometry of the tree because u and d do not depend on a. The probabilities on the branchesemanating from nodes at time ¢ are?
efO8t gy
p=—— = and 1-p=
y — efl5
u—d
The rest of the way that we use the tree is the same as before, except that when discounting betweentunes ¢ and (+ dt we use f(t). A similar modification of the basic tree can be used to value indexoptions, foreign exchange options, and futures options. The dividend yield on an index or a foreignrisk-free rate can be made a function of time by using a similar approach to that just described.
Control Variate Technique
A technique known as the control variate technique can be used for the evaluation of an Americanoption.® This involves using the same tree to calculate both the value of the American option, fy,and the value of the corresponding European option, fg. We also calculate the Black-Scholes priceof the European option, fgg. The error given by the tree in the pricing of the European option isassumed equal to that given by the tree in the pricing of the American option. This gives theestimate of the price of the American option to be
fat fos — Se
To illustrate this approach, Figure 18.10 values the option in Figure {8.3 on the assumption that itis European. The price obtained is $4.32. From the Black-Scholes formula, the true Europeanprice of the option is $4.08. The estimate of the American price in Figure 18.3 is $4.49. The controlvariate estimate of the American price is therefore
4.49 + 4.08 — 4.32 = 4.25
A good estimate of the American price, calculated using L00 steps, is 4.278. The controi variateapproach, therefore, does produce a considerable improvement over the basic tree estimate of 4.49in this case.
The control variate technique in effect involves using the tree to calculate the difference betweenthe European and the American price rather than the American price itself. We give a furtherapplication of the control variate technique when we discuss Monte Carlo: simulation later in thechapter.
18.5 ALTERNATIVE PROCEDURES FOR CONSTRUCTING TREES
 
 
The original Cox, Ross, and Rubinstein approach is not the only way of building a binomial tree.Instead of imposing the assumption u = 1/d on equations (18.2) and (18.3), we can set p = 0.5.
* For a sufficiently large number of time steps, these probabilities are always positive.
® See J. Hull and A. White, “The Use of the Control Variate Technique in Option Pricing,” Journal of Financial andQuantitative Analysis, 23 (September 1988), 237-51.Numerical Procedures 407
At each node:
Upper value = Underlying Asset PriceLower value = Option Price
Shading indicates where option is exercised
Strike price = 50Discount factor per step = 0.9917
Time step, dt = 0.0833 years, 30.42 daysGrowth factor per step, a = 1.0084Probability of up move, p = 0.5073Up step size, u = 1.1224Down step size, d = 0.8909
         
 
 
 
 
 
 
 
 
 
Node Time:
aacno nn 9aaq n4Aanz nOnnD n3200 OATVUE VU.UO09 VU.100F V.ZOUU V.Id00 UF10F
Figure18.10 Tree. as produced by DerivaGem, for the European version of the option in Figure 18.3.At each node, the upper number is the stock price, and the lower number is the option price
A solution to the equations when terms of order higher than 6¢ are ignored is then
v= po 12) bio bi d= ot-8 281-01When the stock provides a dividend yield at rate g, the variable r becomes r — q in these formulas.This allows trees with p = 0.5 to be built for options on indives, foreign exchange, and futures.This alternative tree-building procedure has the advantage over the Cox, Ross, and Rubinsteinapproach that the probabilities are always 0.5 regardless of the value of o or the number of timesteps.’ Its disadvantage is that it is not as easy to calculate delta, gamma, and rho from the treebecause the value of the underlying asset at the central node at time 2 dt is no longer the same as attime zero.
7 When time steps are so large that o < |(r — q)Vét|, the Cox, Ross, and Rubinstein tree gives negative probabilities.The alternative procedure described here does not have that drawback.408 CHAPTER 18
Example18.6 Consider a nine-month American call option on the Canadian dollar. The currentexchange rate is 0.7900, the strike price is 0.7950, the U.S. risk-free interest rate is 6% per annum, theCanadian risk-free interest rate is 10% per annum, and the volatility of the exchange rate is 4% perannum. In this case, Sy = 0.79, K = 0.795, r = 0.06, r7 = 0.10, o = 0.04, and T = 0.75. We divide thelife of the option into three-month periods for the purposes of constructing the tree, so thatbt = 0.25, We set the probabilities on each branch to 0.5 and
 
1 = el". 06-0-10-0.0016/2)0.2540.04V0.25 __ 1 Qggg
d= el!98-0.10-0.0016/2)0.25-0.04V0.25 _ 4 9793
The tree for the exchange rate is shown in Figure 18.11 and tree gives the value of the optionas 30.0026.
Trinomial Trees
Trinomial trees can be used as an alternative to binomial trees. The general form of the tree is asshown in Figure 18.12. Suppose that p,, p,,, and p, are the probabilities of up, middle, and down
At each node:
Upper value = Underlying Asset PriceLower value = Option Price
Shading indicates where option is exercised
Strike price = 0.795Discount factor per step = 0.9851Time step, dt = 0.2500 years, 91.25 days
Probability of up move, p = 0.5000
          
8(eSa10
 
0.000
 
Node Time:0.0000 0.2500 0.5000 0.7500
Figure 18.11 Binomial tree for an American call option on the Canadian dollar. At each node, theupper number is spot exchange rate, and the Jower number is option price. All probabilities are 0.5Numerical Procedures 409
reso
Squ2° Sou
X
Sou
 
 
Sod
Sod?
/\
Sod3Figure 18.12 Trinomial stock price tree
movements at each node and St is the length of the time step. For a non-dividend-paying stock,parameter values that match the mean and standard deviation of price changes when terms oforder higher than 8t are ignored are
!u=e® 38r d=-u
ét btPa=nV TG SOF Pm Pua 30) +5
For a stock paying a dividend yield at rate g, we replace r by r—q in these equations.Calculations for a trinomial tree are analogous to those for a binomial tree. The trinomial treeapproach proves to be equivalent to the explicit finite difference method, which will be describedin Section 18.8.
 
The Adaptive Mesh Model
Figlewski and Gao have proposed a method, which they call the adaptive mesh model, for buildingtrees where a high-resolution (small-ét) tree is grafted on to a low-resolution (large-8t) tree.® Itturns out that, for an American option, we need high resolution close to the strike price and close
5 See §. Figlewski and B. Gao, “The Adaptive Mesh Model: A New Approach to Efficient Option Pricing,” Journalof Financial Economics, 53 (1999), 313-51.410 CHAPTER 18
<— Strike Price
 
Figure 18.13 Adaptive mesh model for an American-style option
to maturity. Figure 18.13 indicates how we can achieve this. In this figure a trinomial tree with atime step of 6t/4 has been grafted on to a trinomial tree with a time step of dt. The calculation ofprobabilities on branches and the rollback procedure is analogous to that for a regular tree. Thetree in Figure 18.13 leads to a significant improvement in numerical efficiency over a binomial ortrinomiai tree.
18.6 MONTE CARLO SIMULATION
 
 
We now move on to discuss Monte Carlo simulation. This uses the risk-neutral valuation result.The expected payoff in a risk-neutral world is calculated using a sampling procedure. It is thendiscounted at the risk-free interest rate.
Consider a derivative dependent on a single market variable S that provides a payoff at time T.Assuming that interest rates are constant, we can value the derivative as follows:”
1. Sample a random path for S$ in a risk-neutral world,2. Calculate the payoff from the derivative.
° We discuss how Monte Carlo simulation can be used with stochastic interest rates in Section 21.4.Numerical Procedures 411
3. Repeat steps | and 2 to get many sample values of the payoff from the derivative in a risk-neutral world.
4. Calculate the mean of the sample payoffs to get an estimate of the expected payoff in a risk-neutral world.
5. Discount the expected payoff at the risk-free rate to get an estimate of the value of thederivative.
Suppose that the process followed by the‘underlying market variable in a risk-neutral world isdS = fiSdt + oS dz (18.13) ©
where dz is a Weiner process, /i is the expected return in a risk-neutral world, and o is thevolatility.'° To simulate the path followed by S, we divide the life of the derivative into N shortintervals of length ét and approximate equation (18.13) as
S(t + dt) — S(t) = S(t) 8t + oS(t)eV5E (18.14)
where S(t) denotes the value of S at time t, and ¢ is a random sample from a normal distributionwith mean zero and standard deviation 1.0. This enables the value of S at time dt to be calculatedfrom the initial value of S, the value at time 2 5¢ to be calculated from the value at time dt, and soon. An illustration of the procedure is given in Section 11.4. One simulation trial involvesconstrucling a complete path for S using N random samples from a normal distribution.
In practice, it is usually more accurate to simulate In $ rather than S. From [t6’s lemma theprocess followed by In S is
dinS= (a-S)atoas
so that2In S(t-+ 8¢) — In S(t) = (« - 5) &t + oeJ/ 51or equivalently
S(t + 5t) = S(t) LG - 2) ot+ aevii] (18.15)
This equation is used to construct a path for S$ in a similar way to equation (18.14). Whereasequation (18.14) is true only in the limit as d¢ tends to zero, equation (18.15) is exactly true for all 5t.
The advantage of Monte Carlo simulation is that it can be used when the payoff depends on thepath followed by the underlying variable S as well as when it depends only on the final value of S.Payoffs can occur at several times during the life of the derivative rather than all at the end. Anystochastic process for S can be accommodated. As will be shown shortly, the procedure can also beextended to accommodate situations where the payoff from the derivative depends on severalunderlying market variables. The drawbacks of Monte Carlo simulation are that it is computa-tionally very time-consuming and cannot easily handle situations where there are early exerciseopportunities."!
10 Tf § is the price of a non-dividend-paying stock, ji = r; if it is an exchange rate, fi =r — 7; and so on. Note thatthe volatility is the same in a risk-neutral world as in the real world, as shown in Section 10.7.
‘| A number of researchers have suggested ways Monte Carlo simulation can be extended to value American options(see Section 20.9).412 CHAPTER 18
 
In the particular case where a derivative provides a payoff at time T dependent only on the valueof S$ at that time, it is not necessary to sample a whole path for S. Instead we can jump straightfrom the value of S at time zero to its value at time 7. When the process in equation (18.13) isassumed,
S(T) = S(0) exp] ( - 5) T+ acv7|
Derivatives Dependent on More than One Market Variable
Consider the situation where the payoff from a derivative depends on x variables, 6, (1 <i <2).Define s; as the volatility of 6,, #, as the expected growth rate of 6; in a risk-neutral world, and px,as the instantaneous correlation between 6, and 6,.'" As in the single-variable case, the life of thederivative must be divided into N subintervals of length 5¢. The discrete version of the process for9; is then
6,(t + St) — 6,(t) = m,0,(t) bt + 9,0,(t)e,Vt (18.16)
where ¢; is a random sample from a standard normal distribution. The coefficient of correlationbetween €; and e& is py for 1 <i,k <n. One simulation trial involves obtaining N samples of thee, (1 <i <a) from a multivariate standardized normal distribution. These are substituted intoequation (18.16) to produce simulated paths for each 0,, thereby enabling a sample value for thederivative to be calculated.
Generating the Random Samples from Normal Distributions
Most programming languages incorporate routines for sampling a random number between zeroand one. An approximate sample from a univariate standardized normal distribution can beobtained from the formula
12e=)°R,-6 (18.17)i=l
where the R; (1 <i < 12) are independent random numbers between zero and one, and e is therequired sample from (0, 1). This approximation is satisfactory for most purposes.
When two correlated samples €, and €) from standard normal distributions are required, anappropriate procedure is as follows. Independent samples x, and x. from a univariate standardizednormal distribution are obtained as just described. The required samples €, and ¢€, are thencalculated as follows:
eax€2 = px, ty V1 -
where p is the coefficient of correlation.
Consider next the situation where we require correlated samples from normal distributionswhere the coefficient of correlation between sample i and sample j is p;;. We first sample nindependent variables x; (1 <i<z), from univariate standardized normal distributions. The
"2 Note that sj, Mj, and py are not necessarily constant; they may depend on the 6;.Numerical Procedures 413
required samples are €,; (I <i <a”), where
k=
For ¢; to have the correct variance and the correct correlation with the €, (1 < j < i), we must have
L2y Oi =k=l
and, for all j <i,iSe einer = Pik=l
The first sample, €,, is set equal to x,. These equations for the @’s can be solved so that €, iscalculated from x, and x2, €3 is calculated from x), x» and x3, and so on.'? The procedure is knownas the Cholesky decomposition.
Number of Trials
The number of simulation trials carried out depends on the accuracy required. If M independenttrials are carried out as described above, it is usual to calculate the standard deviation as well as themean of the discounted payoffs given by the simulation trials for the derivative. Denote the meanby u and the standard deviation by . The variable y is the simulation’s estimate of the value ofthe derivative. The standard error of the estimate is
2VMA 95% confidence interval for the price f of the derivative is therefore given by1.96@ 1.96@bom we foe t+ ——
This shows that our uncertainty about the value of the derivative is inversely proportional to thesquare root of the number of trials. To double the accuracy of a simulation, we must quadruple thenumber of trials; to increase the accuracy by a factor of 10, the number of trials must increase by afactor of 100; and so on.
Bo Nm ke neAPpiuicatigns
Monte Carlo simulation tends to be numerically more efficient than other procedures when thereare three or more stochastic variables. This is because the time taken to carry out a Monte Carlosimulation increases approximately linearly with the number of variables, whereas the time takenfor most other procedures increases exponentially with the number of variables. One advantage ofMonte Carlo simulation is that it can provide a standard error for the estimates that it makes.Another is that it is an approach that can accommodate complex payoffs and complex stochasticprocesses. It can also be used when the payoff depends on some function of the whole path
'8 Tf the equations for the «’s do not have real solutions, the assumed correlation structure is internally inconsistent,as explained in Section 17.7.414 CHAPTER 18
followed by a variable, not just its terminal value. As already noted, a limitation of the MonteCarlo sunulaiion approach is that it is difficult to use it for non-European-siyle derivatives.
Calculating the Greek Letters
The Greek letters discussed in Chapter 14 can be calculated using Monte Carlo simulation.Suppose that we are interested in the partial derivative of f with x, where f is the value of thederivative and x is the value of an underlying variable or a parameter. First, Monte Carlosimulation is used in the usual way to calenlate an estimate, f, for the value of the derivative.A small increase dx is then made in the value of x, and a new value for the derivative, f*, icalculated in the same way as f. An estimate for the hedge parameter is given by
 
 
f-fOx
 
In order to minimize the standard error of the estimate, the number of time intervals N, therandom number streams, and the number of trials M should be the same for calculating both fand f*.
Sampling through a Tree
Instead of implementing Monte Carlo simulation by randomly sampling from the stochasticprocess for an underlying variable, we can sample paths for the underlying variable using abinomial tree. Suppose we have a binomial tree where the probability of an “‘up”” movement is 0.6.The procedure for sampling a random path through the tree is as follows. At each node, we samplea random number between zero and one. If the number is less than 0.4, we take the down path. If itis greater than 0.4, we take the up path. Once we have a complete path from the initial node to theend of the tree, we can calculate a payoff. This completes the first trial. A similar procedure is usedto complete more trials. The mean of the payoffs is discounted at the risk-free rate to get anestimate of the value of the derivative.
148.7 VARIANCE REDUCTION UR
PROCEDURES
 
Tf the simulation is carried out as described so far, a very large value of M is usually necessary toestimate f with reasonable accuracy. This is very expensive in terms of computation time. In thissection, we examine a number of variance reduction procedures that can lead to dramatic savingsin computation time.
Antithetic Variable Technique
In the antithetic variable technique, a simulation trial involves calculating two values of thederivative. The first value, f;, is calculated in the usual way; the second value, fy, is calculatedby changing the sign of all the random samples from standard normal distributions. (If € is asample used to calculate f,, —¢ is the corresponding sample used to calculate f,.) The sample valueof the derivative calculated from a simulation trial is the average of f, and f). This works wellbecause when one value is above the true value, the other tends to be below, and vice versa.Numerical Procedures 415
Denote f as the average of f, and fp:fit2
ie
f=
 
The final estimate of the value of the derivative is the average of the f’s. If @ is the standarddeviation of the /’s, and M is the number of simulation trials (.e., the number of pairs of valuescalculated), the standard error of the estimate is @//M. This is generally much less than thestandard error calculated using 2M random trials.
Control Variate Technique
We have already given one example of the control variate technique in connection with the use oftrees to value American options (see Section 18.4). The control variate technique is applicablewhen there are two similar derivatives, A and B, Derivative A is the security under consideration;derivative B is similar to derivative A and has an analytic solution available. Two simulationsusing the same random number streams and the same 6¢ are carried out in parallel. The first is usedto obtain an estimate, fx, of the value of A; the second is used to obtain an estimate, fy, of thevalue of B. A better estimate, f,, of the value of A is then obtained using the formula
In=fa- fet fo (18.18)
where fg is the known true value of B. Hull and White provide an example of the use of thecontrol variate technique when evaluating the effect of stochastic volatility on the price of aEuropean call option.'* In this case, f, is the estimated value of the option assuming stochasticvolatility and fg is its estimated value assuming constant volatility.
Importance Sampling
Importance sampling is best explained with an example. Suppose that we wish to calculate theprice of a deep-out-of-the-money European call option with strike price K and maturity T. If wesample values for the underlying asset price at time T in the usual way, most of the paths will leadto zero payoff. This is a waste of computation time because the zero-payoff paths contribute verylittle to the determination of the value of the option. We therefore try to choose only importantpaths, that is, paths where the stock price is above K at maturity.
Suppose that F is the unconditional probability distribution function for the stock price attime T and that the probability, k, of the stock price being greater than K at maturity is knownanalytically. Then G = F/k is the probability distribution of the stock price conditional on thestock price being greater than K. To implement importance sampling, we sample from G ratherthan F. The estimate of the value of the option is the average discounted payoff multiplied by k.
Stratified Sampling
Stratified sampling is a way of sampling from the probability distribution of a market variable at afuture time. It involves dividing the distribution into ranges or intervals and sampling from eachinterval according to its probability. Suppose, for example, that there are ten equally likelyintervals. We choose a sampling scheme that ensures that exactly 10% of our samples are from
14 See J. Hull and A. White, “The Pricing of Options on Assets with Stochastic Volatilities,” Journal of Finance,42 (June 1987), 281-300.416 CHAPTER 18
 
the first interval, 10% from the second interval, and so on. If a Jarge number of intervals are used,we can regard either the mean or the median value of the variable, conditional on it being in aninterval, as a representative value for the interval. When sampling from an interval, we then alwayspick the representative value.
Curran shows results from using this procedure to value both European call options and path-dependent options.'> In the case of a standard normal distribution when there are n intervals, wecan calculate the representative value for the ith interval as
wifi 0s)n)
where N' is the inverse cumulative normal distribution. Kor example, when n = 4 the representativevalues corresponding to the four intervals are N_'(0.125), N7'(0.375), N~1(0.625), N”1(0.875). Thefunction N~' can be calculated iteratively using the approximation to the N given in Section 12.9 (orthe NORMSINV function in Excel). An alternative approach is suggested by Moro.'°
 
Moment Matching
Moment matching involves adjusting the samples taken from a standardized normal distributionso that the first, second, and possibly higher moments are matched. Suppose that the normaldistribution samples used to calculate the change in the value of a particular variable over aparticular time period are €; (1 < i <n). To match the first two moments, we calculate the mean ofthe samples, m, and the standard deviation of the samples, ». We then define adjusted samples
yi 1 <i<a)yas€é;-—m
 
y=SsThese adjusted samples have the correct mean of zero and the correct standard deviation of 1.0.We use the adjusted samples for all calculations.
Moment matching saves computation time, but can lead to memory problems because everynumber sampled must be stored until the end of the simulation. Moment matching is sometimestermed quadratic resampling. It is often used in conjunction with the antithetic variable technique.Because the latter automatically matches all odd moments, the goal of moment matching then
becomes that of matching the second moment and, possibly, the fourth moment.
Using Quasi-Random Sequences
A quasi-random sequence (also called a low-discrepancy sequence) is a sequence of representativesamples from a probability distribution.'’ Descriptions of the use of quasi-random sequences havebeen given by Brotherton-Ratcliffe and by Press et a/.'* Quasi-random sequences can have thedesirable property that they lead to the standard error of an estimate being proportional to 1/M
rather than 1 /VM, where M is the sample size.
'S See M. Curran, “Strata Gems,” RISK, March 1994, pp. 70-71.‘6 See B. Moro, “The Full Monte.” R/SK, February 1985, pp. 57-58.? The term quasi-random is a misnomer. A quasi-random sequence is totally deterministic.
'8 See R. Brotherton-Ratcliffe, “Monte Carlo Motoring,” RISK, December 1994, pp. 53-58; W. H. Press, S. A.Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in C: The Art of Scientific Computing, 2nd edn.,Cambridge University Press, Cambridge, 1992.Numerical Procedures 417
 
 
  
 
 
 
 
 
 
 
 
 
LOT yr MOET TT UTETTTT ETT LTTost” + oosh[> * | Ee0.6 7 0.6-04F 4 04h.E fo Piet02h, 4 02F -OCs tr rte tT ghee0 02 04 06 0.8 1.0 0 0.2 04 06 08 10Points 1 to 128 Points 129 to 512
 
 
TIT TT ETT,
ee
     
 
     
q
 
 
 
 
Fie bi Coes 0 “Le0 02 #04 06 08 1.0 0 0.6
Points 513 to 1024 Points 1 to 1024
 
 
  
“04
Figure 18.14 First 1024 points of a Soboi’ sequence
Quasi-random sampling is similar to stratified sampling. The objective is to sample representa-tive values for the underlying variables. In stratified sampling, it is assumed that we know inadvance how many samples will be taken. A quasi-random sampling scheme is more flexible. Thesamples are taken in such a way that we are always “filling in” gaps between existing samples. Ateach stage of the simulation, the sampled points are roughly evenly spaced throughout theprobability space.
Figure 18.14 shows points generated in two dimensions using a procedure suggested by Sobol’.!?It can he seen that successive points do tend to fill in the gaps left by previous points.
Representative Sampling through a Tree
As explained previously, we can implement Monte Carlo simulation by sampling paths through atree. In the spirit of stratified sampling, we can choose representative paths through the tree instead
'? See I. M. Sobol’, USSR Computational Mathematics and Mathematical Physics, 7, no. 4 (1967), 86-112. Adescription of the Sobol’ procedure is given by W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flanncry,Numerical Recipes in C: The Art of Scientific Computing, 2nd edn., Cambridge University Press, Cambridge, 1992.418 CHAPTER 18
of random paths. A set of paths is representative if the proportion of paths going through anygiven node equals (or is close to) the probability of that node being reached. The construction ofrepresentative paths is discussed by Mintz.
The first step is to choose the total number of sample paths. We can then calculate the expectednumber of sample paths going through each node. (This is the probability of the node beingreached times the total number of sample paths.) Suppose that w,; is the expected number of sample
paths going through node (/, j). Typically the u,, are not integers. It is necessary to develop an
algorithm for the “integerization” of each uj (i.e., the decision as to whether each w;; is rounded upor down). The result of the algorithm should be that the correct number of sample paths leaves theinitial node and that the number of sample paths entering each subsequent node equals the numberof sample paths leaving the node. Once we have determined the number of sample paths that willpass along each branch, a “sampling without replacement” procedure can be used to define thesample paths. . ‘
18.8 FINITE DIFFERENCE METHODS
 
 
Finite difference methods value a derivative by solving the differential equation that the derivativesatisfies. The differential equation is