. 12
( 14)



1992 1993 1994 1995 1996 1997

3 month interest rates

1992 1993 1994 1995 1996 1997

Figure 15.4. Interest rates time series: German (thick dotted line),
Japanese (thin dotted line), American (thick straight line), Thai (thin
straight line). XFGbasket.xpl

self-¬nancing strategy which produces a positive expected payo¬:

• at time t
“ borrow the portfolio (1+rj )’1 E(θj,t+h |Ft )Xj,t from the countries
composing the basket,
“ lend (1 + r0 )’1 Yt to Thailand,
“ invest the di¬erence (1 + rj )’1 E(θj,t+h |Ft )Xj,t ’ (1 + r0 )’1 Yt in
the numeraire currency at the risk-free rate r1 ,
• at time t + h
“ withdraw the amount Yt+h from Thailand,
15.2 Estimating the coe¬cients of an exchange rate basket 335

E(θj,t+h |Ft )Xj,t+h ,
“ pay back the loan of
“ keep the di¬erence.

The expression for the pro¬t and for its expected value are:
= Yt+h ’ E(θj,t+h |Ft )Xj,t+h
« 
(1 + rj )’1 E(θj,t+h |Ft )Xj,t ’ (1 + r0 )’1 Yt 
+(1 + r1 ) 
« 
(1 + rj )’1 E(θj,t+h |Ft )Xj,t ’ (1 + r0 )’1 Yt  .
E(Πt+h |Ft ) = (1 + r1 ) 

15.2.2 Estimation results

For the implementation of the investment strategy described above one needs
the estimate of the, possibly time-varying, basket weights. The precision of
the estimation has a direct impact on the economic result of the investment.
Therefore, we compare four di¬erent estimators of the basket weights: the
adaptive, the recursive, the window and the Kalman ¬lter estimator using
economic criteria for a one month and for a three month investment horizon.
In particular we compute the average expected pro¬t and the average realized
The adaptive estimation procedure requires three parameters: m, » and µ. The
choice of m0 does not in¬‚uence the results very much and it can be reasonably
set to 30. This value represents the minimal amount of data which are used
for the estimation, and in the case of a structural break, the minimal delay
before having the chance of detecting the change point. The selection of » and
µ is more critical. These two values determine the sensitivity of the algorithm.
Small values would imply a fast reaction to changes in the regressor coe¬cients,
but but they would also lead to the selection of intervals of homogeneity which
are possibly too small. Large values would imply a slower reaction and con-
sequently the selection of intervals which can be too large. To overcome this
problem we suggest the following approach.
The main idea is that small changes in the values of » and µ should not a¬ect
the estimation results. Therefore we restrict our attention on a set S of possible
336 15 Locally time homogeneous time series modeling

USD weight


1994 1995 1996 1997

DEM weight

1994 1995 1996 1997

JPY weight

1994 1995 1996 1997
Figure 15.5. Estimated exchange rate basket weights: adaptive
(straight line), recursive (thine dotted line), window (thick dotted line).

pairs (», µ). In the present context we chose all the even number between 2
and 8:
S = {(», µ)| », µ ∈ {2, 4, 6, 8}}
15.2 Estimating the coe¬cients of an exchange rate basket 337

Then we compare the 16 pairs with the following criterion at each time t:
t’1 d
(»— , µ— ) = arg min Ys ’ θj,s|s’h Xj,s  .
(»,µ)∈S s=t’200 j=1

Finally, we estimate the value of θt+h|t with the selected pair (»— , µ— ). The
appeal of the above selection criterion consists of the fact that it leads to the
choice of the pair (», µ) which has provided the least quadratic hedging costs
over the past trading periods. Notice that in general we have di¬erent results
depending on the length of the forecasting horizon: here one and three month.
Figure 15.5 shows the results for the three month horizon. It is interesting to
see that the adaptive estimate tends to coincide with the recursive estimate
during the ¬rst half of the sample, more or less, while during the second half
of the sample it tends to follow the rolling estimate.
We remark that the problem of selecting free parameters is not speci¬c to the
adaptive estimator. The window estimator requires the choice of the length
of the window: k, while the Kalman ¬lter needs the speci¬cation of the data
generating process of θt and the determination of Σ and σ. In this application
k is set equal to 250, Σ and σ are estimated recursively from the data using the
OLS, while θ0|0 and P0|0 are initialized using the ¬rst 350 observations which
are then discarded. We remark that this choice is consistent with the one of
Christo¬ersen and Giorgianni (2000).
Table 15.2 shows the result of the simulated investment. The investments are
normalized such that at each trading day we take a short position of 100 USD
in the optimal portfolio of the hard currencies. The result refers to the period
April 9 1993 to February 12 1997 for the one month horizon investment and
June 7 1993 to February 12 1997 for the three month horizon investment. Notice
¬rst that the average realized pro¬ts are positive and, as far as the three month
investment horizon is concerned, they are signi¬cantly larger than zero among
all methods. This provides a clear evidence for the fact that arbitrage pro¬ts
were possible with in the framework of the Thai Bath basket for the period
under study. The comparison of the estimator also show the importance of
properly accounting for the time variability of the parameters. The recursive
estimator shows modest result as far as the realized pro¬ts are concerned and
the largest bias between expected the realized pro¬t. On one side, the bias is
reduced by the window estimator and by the Kalman ¬lter, but on the other
side these two methods provide a worse performance as far as the realized pro¬t
are concerned. Finally, the adaptive estimator appears to be the best one, its
338 15 Locally time homogeneous time series modeling

ONE MONTH HORIZON Recursive Window KF Adaptive

Average Expected Pro¬ts .772 .565 .505 .553

Average Realized Pro¬t .403 .401 .389 .420

Standard errors (.305) (.305) (.330) (.333)

THREE MONTH HORIZON Recursive Window KF Adaptive

Average Expected Pro¬ts 1.627 1.467 1.375 1.455

Average Realized Pro¬t 1.166 1.141 1.147 1.182

Standard errors (.464) (.513) (.475) (.438)

Table 15.2. Summary statistics of the pro¬ts.

bias is much smaller than the one of the recursive estimator and it delivers the
largest realized pro¬ts for both investment horizons.

15.3 Estimating the volatility of ¬nancial time
The locally time homogeneous approach appears to be also appropriate for the
estimation of the volatility of ¬nancial time series. In order to provide some
motivation we ¬rst describe the stylized facts of ¬nancial time series. Let St
de¬ne the price process of a ¬nancial asset such as stocks or exchange rates,
then the returns are de¬ned as follows:

Rt = ln St ’ ln St’1 .

Stylized facts of ¬nancial asset returns are: a leptokurtic density, variance
clustering and highly persistent autocorrelation of square and absolute returns
(see Figure 15.6). Further details and examples on this topic can be found in
Taylor (1986) and in Franke, H¨rdle and Hafner (2001).
15.3 Estimating the volatility of ¬nancial time series 339


0 50 100 150 200 250 300

ACF of the absolute returns

0 10 20 30 40 50
Figure 15.6. JPY/USD returns XFGretacf.xpl

15.3.1 The standard approach

The returns of ¬nancial time series are usually modeled by the following equa-
Rt = σ t µt
Where σt is a strictly positive process, which describes the dynamics of the vari-
ance of Rt , and ξt has a standard normal distribution: ξt ∼ N (0, 1). Standard
parametric models of the volatility are of (G)ARCH type:
2 2 2
σt = ω + ±Rt’1 + βσt’1 ,
340 15 Locally time homogeneous time series modeling

like in Engle (1995) and Bollerslev (1995), and of stochastic volatility type:
2 2
ln σt = θ0 + θ1 ln σt’1 + νt ,

as described by Harvey, Ruiz and Shephard (1995). These models have been
expanded in order to incorporate other characteristics of the ¬nancial return
time series: TARCH, EGARCH and QARCH explicitly assume an asymmet-
ric reaction of the volatility process to the sign of the observed returns, while
IGARCH and FIGARCH model the long memory structure of the autocorre-
lations of the square returns.

15.3.2 The locally time homogeneous approach

A common feature to all the models which have been cited in the previous
section is that they completely describe the volatility process by a ¬nite set of
parameters. The availability of very large samples of ¬nancial data has given
the possibility of constructing models which display quite complicated param-
eterizations in order to explain all the observed stylized facts. Obviously those
models rely on the assumption that the parametric structure of the process
remains constant through the whole sample. This is a nontrivial and possi-
bly dangerous assumption in particular as far as forecasting is concerned as
pointed out in Clements and Hendry (1998). Furthermore checking for param-
eter instability becomes quite di¬cult if the model is nonlinear, and/or the
number of parameters is large. Whereby those characteristics of the returns
which are often explained by the long memory and (fractal) integrated nature
of the volatility process, could also depend on the parameters being time vary-
ing. We want to suggest an alternative approach which relies on a locally time
homogeneous parameterization, i.e. we assume that the volatility σ follows a
jump process and is constant over some unknown interval of time homogeneity.
The adaptive algorithm, which has been presented in the previous sections,
also applies in this case; its aim consists in the data-driven estimation of the
interval of time homogeneity, after which the estimate of the volatility can be
simply obtained by local averaging.

15.3.3 Modeling volatility via power transformation

Let St be an observed asset process in discrete time, t = 1, 2, . . . , „ and Rt
are the corresponding returns: Rt = log(St /St’1 ) . We model this process via
15.3 Estimating the volatility of ¬nancial time series 341

the conditional heteroscedasticity assumption

Rt = σt µt , (15.10)

where µt , t ≥ 1 , is a sequence of independent standard Gaussian random
variables and σt is the volatility process which is in general a predictable
random process, that is, σt is measurable with respect to Ft’1 with Ft’1 =
σ(R1 , . . . , Rt’1 ) .
The model equation (15.10) links the volatility σt with the observations Rt
via the multiplicative errors µt . In order to apply the theory presented in
Section 15.1 we need a regression like model with additive errors. For this
reason we consider the power transformation, which leads to a regression with
additive noise and so that the noise is close to a Gaussian one, see Carroll and
Ruppert (1988). Due to (15.10) the random variable Rt is conditionally on
Ft’1 Gaussian and it holds
2 2
E Rt |Ft’1 = σt .

Similarly, for every γ > 0 ,
γ γ γ
E Rt |Ft’1 = σt E (|ξ|γ |Ft’1 ) = Cγ σt ,
γ 2
γ 2γ 2γ 2
= σt E (|ξ|γ ’ Cγ ) = σt Dγ
’ Cγ σt |Ft’1
E Rt

where ξ denotes a standard Gaussian r.v., Cγ = E|ξ|γ and Dγ = Var|ξ|γ .

Therefore, the process |Rt |γ allows for the representation
γ γ
|Rt |γ = Cγ σt + Dγ σt ζt , (15.11)

where ζt is equal (|ξ|γ ’ Cγ ) /Dγ . A suitable choice of the value of γ provides
that the distribution of
(|ξ|γ ’ Cγ ) /Dγ
is close to the normal. In particular the value of γ = 0.5 appears to be almost
optimal, see Figure 15.7.

15.3.4 Adaptive estimation under local time-homogeneity

The assumption of local time homogeneity means that the function σt is
constant within an interval I = [„ ’ m, „ ] , and the process Rt follows the
342 15 Locally time homogeneous time series modeling

Normal and powertransformed densities


-1 0 1 2
Figure 15.7. Normal and power transformed densities for γ = 0.5.

regression-like equation (15.11) with the constant trend θI = Cγ σI which can
be estimated by averaging over this interval I :
|Rt |γ .
θI = (15.12)

By (15.11)

Cγ Dγ 1 sγ
γ γ
θI = σt + σt ζt = θt + θt ζt (15.13)
|I| |I| |I| |I|
t∈I t∈I t∈I t∈I

with sγ = Dγ /Cγ so that

EθI = E θt , (15.14)
s2 s2
γ γ 2
E θt ζt = E θt . (15.15)
|I|2 |I|2
t∈I t∈I

De¬ne also
2 2
vI =2 θt .
15.3 Estimating the volatility of ¬nancial time series 343

In view of (15.15) this value is called the conditional variance of θI . Under
local homogeneity it holds θt is constantly equal to θI for t ∈ I , and hence,

EθI = θI
s2 θI
vI = Var θI = .

A probability bound analogous to the one in Section 15.1 holds also in this
case. Let the volatility coe¬cient σt satisfy the condition b ¤ σt ¤ bB with
some constants b > 0, B > 1 . Then there exists aγ > 0 such that it holds for
every » ≥ 0


P |θI ’ θ„ | > »vI ¤ 4 e»(1 + log B) exp ’ . (15.16)

The proof of the statement above and some related theoretical results can be
found in Mercurio and Spokoiny (2000).
For practical application one has to substitute the unknown conditional stan-
dard deviation with its estimate: vI = sγ θI |I|’1/2 . Under the assumption of
time homogeneity within an interval I = [„ ’ m, „ ] equation (15.16) allows
to bound |θI ’ θJ | by »vI + µvJ for any J ‚ I, provided that » and µ
are su¬ciently large. Therefore we can apply the same algorithm described in
Section 15.1 in order to estimate the largest interval of time homogeneity and
the related value of θ„ . Here, as in the previous section, we are faced with
the choice of three tuning parameters: m0 , », and µ. Simulation studies and
repeated trying on real data by Mercurio and Spokoiny (2000) have shown that
the choice of m0 is not particularly critical and it can be selected between 10
and 50 without a¬ecting the overall results of the procedure.
As described in Section 15.2.2, the choice of » and µ is more delicate. The
in¬‚uence of » and µ is similar to the one of the smoothing parameters in
the nonparametric regression. The likelihood of rejecting a time homogeneous
interval decreases with increasing » and/or µ. This is clear from equation
(15.6). Therefore if » and µ are too large this would make the algorithm too
conservative, increasing the bias of the estimator, while too small values of
» and µ would lead to a frequent rejection and to a high variability of the
estimate. Once again, a way of choosing the optimal values of » and µ can be
made through the minimization of the squared forecast error. One has to de¬ne
a ¬nite set S of the admissible pair of » and µ. Then for each pair belonging
344 15 Locally time homogeneous time series modeling

to S one can compute the corresponding estimate: θt and then select the
optimal pair and the corresponding estimate by the following criterion:
T 2
|Rt |γ ’ θt
(», µ) = min .

Figure 15.8 shows the result of the on-line estimation of the locally time homo-
geneous volatility model for the JPY/USD exchange rate. The bottom plot, in
particular, shows the estimated length of the interval of time homogeneity: m,
at each time point.

15.4 Technical appendix
In this section we give the precise conditions under which the bound (15.4)
holds. De¬ne:
WI = VI’1 ,
VI = σ ’2 Xt Xt

furthermore let wij,I denote the√ elements of WI . For some positive constants
b > 0, B > 1, ρ < 1, r ≥ 1, » > 2 and for i = 1 . . . p consider the random set
were the following conditions are ful¬lled:
± 
b ¤ wjj,I ¤ bB
 
 
 
 
 
 
||VI µ||2 wjj,I ¤ r
Ai,I =
 {µ∈RK :||µ||=1} 
 
 
 
 
|wji,I /wjj,I | ¤ ρ ∀ i = 1, . . . , p
 

Let (Y1 , X1 ) . . . (Y„ , X„ ) obey (15.1), where the regressors are possibly stochas-
tic, then it holds holds for the estimate θI :

P |θi,I ’ θi,„ | > » wii,I ; Ai,I

¤ 4e ln(4B)(1 + 2ρ r(d ’ 1)»)p’1 » exp(’»2 /2), i = 1, . . . , p.

A proof of this statement can be found in Liptser and Spokoiny (1999). For
a further generalization, where the hypothesis of local time homogeneity holds
only approximatively, see H¨rdle et al. (2000).
15.4 Technical appendix 345

YEN/DM returns

0 5 10 15 20


0 5 10 15 20

Interval of homogeneity

0 5 10 15 20
Figure 15.8. From the top: returns, estimated locally time homoge-
neous volatility and estimated length of the interval of time homogene-
ity. XFGlochom.xpl

Bollerslev, T. (1995). Generalised autoregressive conditional heteroskedasticity,
in Engle (1995).
346 15 Locally time homogeneous time series modeling

Carroll, R. and Ruppert, D. (1988). Transformation and Weighting in Regres-
sion, Chapman and Hall, New York.

Christo¬ersen, P. and Giorgianni, L. (2000). Interest rate in currency basket:
Forecasting weights and measuring risk, Journal of Business and Eco-
nomic Statistics 18: 321“335.
Chui, C. and Chen, G. (1998). Kalman Filtering, Information Sciences, third
edn, Springer-Verlag, Berlin.
Clements, M. P. and Hendry, D. F. (1998). Forecastng Economic Time Series,
Cambridge University Press, Cambridge.
Cooley, T. F. and Prescott, E. C. (1973). An adaptive regression model, Inter-
national Economic Review 14: 364“371.
Eichengreen, B., Masson, P., Savastano, M. and Sharma, S. (1999). Transition
Strategies and Nominal Anchors on the Road to Greater Exchange Rate
Flexibility, number 213 in Essays in International Finance, Princeton Uni-
versity Press.

Elliot, R. J., Aggoun, L. and Moore, J. B. (1995). Hidden Markov Models,
Springer-Verlag, Berlin.
Engle, R. F. (ed.) (1995). ARCH, selected readings, Oxford University Press,
Franke, J., H¨rdle, W. and Hafner, C. (2001). Einf¨hrung in die Statistik der
a u
Finanzm¨rkte, Springer, Berlin.

H¨rdle, W., Herwartz, H. and Spokoiny, V. (2001). Time inhomogeous multiple
volatility modelling. Discussion Paper 7, Sonderforschungsbereich 373,
Humboldt-Universit¨t zu Berlin. To appear in Financial Econometrics.
H¨rdle, W., Spokoiny, V. and Teyssi`re, G. (2000). Adaptive estimation for
a e
a time inhomogeneouse stochastic volatility model. Discussion Paper 6,
Sonderforschungsbereich 373, Humboldt-Universit¨t zu Berlin.
Harvey, A., Ruiz, E. and Shephard, N. (1995). Multivariate stochastic variance
models, in Engle (1995).
Lepski, O. (1990). One problem of adaptive estimation in gaussian white noise,
Theory Probab. Appl. 35: 459“470.
15.4 Technical appendix 347

Lepski, O. and Spokoiny, V. (1997). Optimal pointwise adaptive methods in
nonparametric estimation, Annals of Statistics 25: 2512“2546.

Liptser, R. and Spokoiny, V. (1999). Deviation probability bound for mar-
tingales with applications to statistical estimation, Stat. & Prob. Letter
46: 347“357.
Mercurio, D. and Spokoiny, V. (2000). Statistical inference for time-
inhomogeneous volatility models. Discussion Paper 583, Weierstrass In-
stitute for Applied Analysis and Stochastic, Berlin.
Mercurio, D. and Torricelli, C. (2001). Estimation and arbitrage opportunities
for exchange rate baskets. Discussion Paper 37, Sonderforschungsbereich
373, Humboldt-Universit¨t zu Berlin.
Musiela, M. and Rutkowski, M. (1997). Martingale Methods in Financial Mod-
elling, number 36 in Application of Mathemathics. Stochastic Modelling
and Applied Probability, Springer, New York.
Spokoiny, V. (1998). Estimation of a function with discontinuities via lo-
cal polynomial ¬t with an adaptive window choice, Annals of Statistics
26: 1356“1378.
Taylor, S. J. (1986). Modelling Financial Time Series, Wiley, Chichester.
16 Simulation based Option Pricing
Jens L¨ssem and J¨rgen Schumacher
u u

16.1 Simulation techniques for option pricing
We introduce Monte Carlo techniques and Quasi Monte Carlo techniques for
option pricing. First, we give an idea how to use simulation techniques to
determine option prices, then - using the developed basic methods - we give
examples how to price more complex i.e. exotic options even on more than one
underlying. Finally we present a short guideline how to price exotic options
with the proposed techniques.
First, we take a look at a European put on one underlying stock, a pricing
problem which can be solved analytically e.g. by using the Black-Scholes for-
mula. We start with this problem not only because it has become a kind of
”standard problem” but also to have the possibility to compare the results of
our approximation with an analytical solution. At the same time we look at
the time-complexity of the used simulation technique. Next, we show how to
price path dependent options with Monte Carlo methods. Afterwards, we show
how to price a stock option on several underlyings. This implies that we have
to solve a multi-dimensional simulation problem.

16.1.1 Introduction to simulation techniques

The idea behind randomized algorithms is that a random sample from a pop-
ulation (of input variables) is representative for the whole population. As a
consequence, a randomized algorithm can be interpreted as a probability dis-
tribution on a set of deterministic algorithms.
We will see that there are three main advantages to randomized algorithms:
1. Performance: For many problems, it can be shown that randomized algo-
350 16 Simulation based Option Pricing

rithms run faster than the best known deterministic algorithm. 2. Simplicity:
Randomized algorithms are easier to describe and implement than comparable
deterministic algorithms. 3. Flexibility: Randomized algorithms can be easily
In general one distinguishes two types of randomized algorithms. Las Vegas
algorithms are randomized algorithms that always give correct results with
only the variation from one run to another being its running time. Monte
Carlo algorithms are randomized algorithms that may produce an incorrect
solution for which one can bound the probability of occurrence. The quality of
the solution can be seen as a random variable.
Within this chapter, we focus on Monte Carlo algorithms calculating the value
of the following integral
f (x)dx (16.1)

by evaluation of f (x) for independent uniform distributed random vectors
X1 , X2 , . . . , Xn , Xi ∈ [0, 1]d .
The arithmetic mean of the values f (Xi ) can be seen as a guess for the expected
value of the random variable f (Xi ) and therefore can be interpreted as an
approximation for the value of the integral. According to the strong law of
large numbers the estimator for the expected value (the arithmetic mean of
the random function values) is converging to the expected value (the value of
the integral) with an increasing sample size. The probability that the absolute
error of the approximation result exceeds a ¬xed positive value is limited and
decreases to zero with an increasing sample size if the variance of f is ¬nite.

16.1.2 Pricing path independent European options on one

For the case of a European option on one underlying we have to approximate
the following integral via Monte Carlo simulation:

r(T ’t)
CT (ST )g(ST |St , r, σ, T ’ t)dST (16.2)
e E [CT (ST )|St ] =

CT {f (x, St , r, σ, T ’ t)}dx
= (16.3)
16.1 Simulation techniques for option pricing 351

)(T ’t)))2
exp ’ (log ST ’(log St ’(r’0.5σ
2σ 2 (T ’t)
g(ST |St , r, σ, T ’ t) =
2πσ 2 (T ’ t)ST

is the risk neutral density function of the Black Scholes model with parameters:

ST : price of the underlying at maturity
St : price of the underlying at time t
r : risk free interest rate
σ : volatility of log returns of the underlying
T ’t : time to maturity

ST = f (x, St , r, σ, T ’ t) = St exp{(r ’ σ 2 )(T ’ t) + σ T ’ tF ’1 (x)}
transforms the uniform distributed values x in g(ST |St , r, σ, T ’ t) distributed
underlying values ST . F ’1 (x) is the inverse of the cumulative normal distri-
bution function and CT (y) is the payo¬ function of the option.
The Monte Carlo simulation calculates the value of the integral in the following

1 n
1. n independent random underlying values ST . . . ST are generated by com-
puting f (x, St , r, σ, T ’t) for a set of uniformly distributed pseudo random
numbers X1 , . . . , Xn .
i i
2. the option payo¬ CT (ST ) is calculated for each ST .
3. the value of the integral in (16.3) is then approximated by the arithmetic
mean of the option payo¬s:
¯ i
C= CT (ST )
n i=1

We will now derive an estimate of the approximation error of the arithmetic
1 n
mean. We assume that ST . . . ST are independent random underlying samples
of the g(ST |St , r, σ, T ’ t) density. Using this assumption we can conclude that
C is a random variable with expected value
¯ = er(T ’t) Ct (St )
352 16 Simulation based Option Pricing

Additionally we have to assume that the variance of the option payo¬s CT (ST )
is given by:

CT (ST )2 g(ST |St , r, σ, T ’ t)dST ’ E [CT (ST )]
Var [CT (ST )] = (16.4)

exists. Then we get:
1 1
¯ i
Var C = 2 Var CT (ST ) = Var [CT (ST )] (16.5)
n n

1 n
because of the independence of ST , . . . , ST .
The expected value of the random variable C equals the value of the inte-
gral er(T ’t) Ct (St ) and its variance converges to zero with increasing n. The
probability that the approximation error is greater than a ¬xed positive value
decreases to 0 with an increasing number n. A ¬rst estimation of the error is
given by the Chebychev inequality for C,

n Var [CT (ST )]
P |C ’ er(T ’t) Ct (St )| ≥ a ¤
The bound given by this equation is rather imprecise since we do not make any
assumptions on the distribution of the random variable. Only the expected
value and the variance are used in the previous equation. According to the
central limit theorem the distribution of C converges to a normal distribution
for n ’ ∞. It follows that the di¬erence between the approximation and the
integral, C ’ er(T ’t) Ct (St ) is approximately normally distributed with mean 0
and standard deviation

Var [CT (ST )]
σC = (16.7)
for large n. According to Boyle (1977) a value of n > 1000 is su¬ciently large
in order to use the normal distribution for error estimation purposes.

We get the following equation if we assume that C ’ er(T ’t) Ct (St ) is normal

C ’ er(T ’t) Ct (St ) ¤ =√ exp ’
P du (16.8)
2π ¯

16.1 Simulation techniques for option pricing 353

If we choose k as a multiple of the standard deviation σC of C, then we get:

C ’ er(T ’t) Ct (St )
C ’ er(T ’t) Ct (St ) ¤ kσC ¤k
P =P
=√ exp ’ du
2π ’k
=p (16.9)

Given a ¬xed probability level p, the error converges to zero with O(1/ n).
The error interval holds for k = 1, 2, 3 with the respective probabilities p =
0.682, 0.955, 0.997
The con¬dence intervals for a given probability level depend on the standard
deviation of the payo¬ function CT (ST ):

σ CT = Var [CT (ST )] . (16.10)

In general, this standard deviation cannot be calculated with analytical meth-
ods. Therefore one calculates the empirical standard deviation σ and uses it
as a proxy for the error bounds:

1 2
CT (ST ) ’ C
¯ . (16.11)

Figure 16.1 shows the evolution of the absolute error of the price for a European
call option calculated by Monte Carlo methods compared with √ analyticthe
solution. One can observe that the error tends to zero with O (1/ n).
We would like to give some of the main properties of algorithms using Monte
Carlo techniques. First from (16.9) it follows that the error bound tends to zero

with O (1/ n) for a ¬xed probability level p. Second, the probability that a

¬xed error bound holds converges to 1 with O (1/ n), Mavin H. Kalos (1986).
Since these results hold independent of the dimension of the problem, which
a¬ects only the variance of the payo¬ function with respect to the Black-Scholes
risk neutral density, the Monte Carlo method is especially well suited for the
evaluation of option prices in multidimensional settings. Competing pricing
methods e.g ¬nite di¬erences have exponential growing computational costs in
354 16 Simulation based Option Pricing

Errors in MC Simulation

0.03 0.06 0.09 0.12 0.15
absolute error

0 500000 1000000
number of iterations
Figure 16.1. Absolute error of a European Call option price calculated
by Monte Carlo simulations vs. n’1/2

the dimension of the problem. Another advantage of the Monte Carlo pricing
method is the error estimate given by the empirical standard deviation which
can be computed with a small additional e¬ort.
The two most important drawbacks of the Monte Carlo simulation, mentioned
in literature are its small convergence speed compared to other techniques for
options on few underlyings and the di¬culties occurring for options with early
exercise possibilities. For example, American options giving the investor the
possibility to exercise the option at any time before and at maturity, are di¬cult
to price. To evaluate an American option means to ¬nd an optimal exercise
strategy which leads - using only basic Monte Carlo techniques - to a recursive
algorithm with exponential time-complexity. But more advanced techniques
using importance sampling methods show that Monte Carlo simulations can
be applied to evaluate American contracts, Broadie (2000).

16.1.3 Pricing path dependent European options on one

There are two categories of options. Path-independent options are options
whose payo¬ depends only on the underlying prices at maturity. Path-
16.1 Simulation techniques for option pricing 355

dependent options are options whose payo¬ depends on underlying price out-
comes St1 , . . ., Stm at several time points t1 ¤ . . . ¤ tm within the lifetime of
the respective option.
Within the group of path-dependent options one can distinguish options with
a payo¬ function depending on a continuously de¬ned path variable and op-
tions with a payo¬ function depending on a ¬xed number of underlying values.
The price of an option with many - usually equally spaced - exercise dates is
often approximated by the price of an option with a continuously de¬ned path
variable and vice versa.
Examples for path-dependent options are barrier options, lookback options,
and Asian options. The latter have a payo¬ function which is linked to the
average value of the underlying on a speci¬c set of dates during the life of
the option. One distinguishes two basic forms of Asian options: options on
the geometric mean (for which the price can be calculated with standard tech-
niques) and options on the arithmetic mean (for which the price can not be
determined using standard approaches). Asian options are frequently used in
commodity markets. The volatility of the underlying prices of the commodities
is usually very high so that prices for vanilla options are more expensive than
for comparable Asian-style options.

16.1.4 Pricing options on multiple underlyings

In this section we show how to extend the Monte Carlo simulation technique
to higher dimensions. The problem is not only that one has to deal with higher
dimensional integrals, but also that one has to incorporate the underlying cor-
relation structure between the considered securities. In our framework we need
the covariance matrix of the log returns on an annual basis.
In general, a basket option is an option on several underlyings (a basket of
underlyings). Basket options can be European-, American or even Asian-style
options. Normally, the average of the underlying prices is taken to calculate
the price of the basket option, but sometimes other functions are used.
The advantage of the usage of basket options instead of a series of one dimen-
sional options is that the correlation between securities is taken into account.
This may lead to better portfolio hedges. We will look at a basket option on
¬ve underlyings where the underlying price of the best security in the basket
is taken to calculate the option price.
356 16 Simulation based Option Pricing

16.2 Quasi Monte Carlo (QMC) techniques for
option pricing

16.2.1 Introduction to Quasi Monte Carlo techniques

QMC methods can be considered as an alternative to Monte Carlo simulation.
Instead of (pseudo) random numbers, Quasi Monte Carlo algorithms use the
elements of low discrepancy sequences to simulate underlying values.
The discrepancy of a set of points P ‚ [0, 1]s measures how evenly these points
are distributed in the unit cube. The general measure of discrepancy is given

A(B; P )
’ »s (B)
Dn (B; P ) := sup (16.12)

where A(B; P ) is the number of points in P belonging to B, »s (B) is the
Lebesgue measure of the set B, B is a family of Lebesgue measurable subsets
of [0, 1]s , and n is the number of elements in P .
The discrepancy of a set is the largest di¬erence between the number of points
in a subset and the measure of the subset. If we de¬ne B to be the family J of
subintervals i=1 [0, ui ), then we get a special measure, the star-discrepancy:

Dn (P ) := Dn (J ; P ) (16.13)

16.2.2 Error bounds

For the star-discrepancy measure and reasonable assumption on the nature of
the function that has to be integrated an upper bound on the error is given by
the following theorem:

THEOREM 16.1 (Koksma-Hlawka) If the function f is of ¬nite variation
V (f ) in the sense of Hardy and Krause, then the following equation holds for
all sets of points {x1 , . . . , xn } ‚ I s = [0, 1]s
1 —
f (xi ) ’ f (u)du ¤ V (f )Dn (x1 , . . . , xn ) (16.14)
n Is
16.2 Quasi Monte Carlo (QMC) techniques for option pricing 357

A proof is given in Niederreiter (1992).

This means that the error is bounded from above by the product of the
variation V (f ), which in our case is model and payo¬ dependent and the star-
discrepancy of the sequence. The bound cannot be used for an automatic error
estimation since the variation and the star-discrepancy cannot be computed
easily. It has been shown though that sequences exist with a star-discrepancy
of the order O(n’1 (ln n)s ). All sequences with this asymptotic upper bound
are called low-discrepancy sequences Niederreiter (1992). One particular
low-discrepancy sequence is the Halton sequence.

16.2.3 Construction of the Halton sequence

We start with the construction of the one-dimensional Halton sequence within
the interval [0, 1]. An element of this sequence is calculated by using the fol-
lowing equation:

nk,i p’k’1
xi = (16.15)
with i > 0, p = 2 and nk,i determined by the following equation:

nk,i pk ; 0 ¤ nk,i < p; nk,i ∈ N
i= (16.16)

Note that with the above equation nk,i is a function of i and takes values only
in {0; 1}. To illustrate the algorithm we calculate the ¬rst three points.

i = 1: n0,1 = 1, nk,1 = 0 for every k > 0
i = 2: n1,2 = 1, nk,2 = 0 for every k = 1
i = 3: n0,3 = n1,3 = 1, nk,3 = 0 for every k > 1

Therefore we get the sequence 1/2, 1/4, 3/4, 1/8, 5/8, .... The extension of this
construction scheme to higher dimensions is straightforward. For every dimen-
sion j = 1, . . . , d we de¬ne xj by

xj nk,i (j)p’k’1
= (16.17)
i j
358 16 Simulation based Option Pricing

with pj is the jth smallest prime number and nk,i (j) is calculated as follows:

nk,i (j)pk ; 0 ¤ nk,i (j) < pj ; nk,i (j) ∈ N ∀j
i= (16.18)

By using p1 = 2, p2 = 3 we get the following two-dimensional Halton sequence:
(1/2; 1/3), (1/4; 2/3), .... In contrast to grid discretization schemes like i/n i =
1, ..., n low-discrepancy sequences ¬ll the integration space in an incremental
way avoiding the exponential growth of grid points of conventional schemes.
XploRe provides quantlets to generate pseudo random numbers and low dis-
crepancy sequences. For the generation of the pseudo random numbers we use

erg = randomnumbers (seqnum,d,n)
generates n pseudo random vectors of dimension d

where seqnum is the number of the random generator according to Table 16.1,
d is the dimension of the random vector and n the number of vectors generated.

0 Park and Miller with Bays-Durham shu¬„e
1 L™Ecuyer with Bays-Durham shu¬„e
2 Knuth
3 generator from G. Marsaglia et al. Marsaglia (1993)
4 random number generator of your system
5 generator from ACM TOMS 17:98-111
6 multiply with carry gen. (Marsaglia) Marsaglia (1993)
Table 16.1. Random generator that can be used in XploRe

The generation of low discrepancy sequences is provided by

erg = lowdiscrepancy (seqnum,d,n)
generates the ¬rst n low discrepancy sequence vectors of dimen-
sion d

where seqnum is the number of the low discrepancy sequence according to Table
16.2 Quasi Monte Carlo (QMC) techniques for option pricing 359

0 Halton sequence
1 Sobol sequence
2 Faure sequence
3 Niederreiter sequence
Table 16.2. Low-discrepancy sequences available in XploRe,
(Niederreiter, 1992) .

16.2.4 Experimental results

Figure 16.2 shows that two dimensional Halton points are much more equally
spaced than pseudo random points. This leads to a smaller error at least for
“smooth” functions.

First 1000 random points First 1000 Halton points
dimension 2

dimension 2



0 0.5 1 0 0.5 1
dimension 1 dimension 1
Figure 16.2. 1000 two-dimensional pseudo random points
vs. 1000 Halton points XFGSOPRandomNumbers.xpl,

The positive e¬ect of using more evenly spread points for the simulation task
is shown in Figure 16.3. The points of a low-discrepancy sequence are designed
in order to ¬ll the space evenly without any restrictions on the independence
of sequence points where as the pseudo random points are designed to show no
statistically signi¬cant deviation from the independence assumption. Because
of the construction of the low discrepancy sequences one cannot calculate an
360 16 Simulation based Option Pricing

Errors in QMC vs. MC Simulation

log absolute error

0 500000 1000000
number of iterations
Figure 16.3. Absolute error of a random sequence and the Halton
sequence for a put option

empirical standard deviation of the estimator like for Monte Carlo methods
and derive an error approximation for the estimation. One possible way out
of this dilemma is the randomization of the low-discrepancy sequences using
pseudo random numbers i.e. to shift the original quasi random numbers with
pseudo random numbers Tu¬n (1996). If x1 , . . . , xn are scalar elements of a
low-discrepancy sequence X then we can de¬ne a new low discrepancy sequence

xi + xi + <= 1
W ( ) = {y1 , . . . , yn } with yi = (16.19)
(xi + ) ’ 1 otherwise

for a uniformly distributed value . Then we can calculate an empirical stan-
dard deviation of the price estimates for di¬erent sequences W ( ) for indepen-
dent values which can be used as a measure for the error. Experiments with
payo¬ functions for European options show that this randomization technique
reduces the convergence rate proportionally.
The advantage of the Quasi Monte Carlo simulation compared to the Monte
Carlo simulation vanishes if the dimension increases. Especially the compo-
nents with a high index number of the ¬rst elements in low-discrepancy se-
quences are not evenly distributed Niederreiter (1992). Figure 16.4 shows that
the 49th and 50th component of the ¬rst 1000 Halton points are not evenly
distributed. But the result for the ¬rst 10000 points of the sequence shows that
the points become more evenly spread if the number of points increases.
16.3 Pricing options with simulation techniques - a guideline 361

First 1000 Halton Points First 10000 Halton Points

dimension 50

dimension 50



0 0.5 1 0 0.5 1
dimension 49 dimension 49
Figure 16.4. The ¬rst 1000 and 10000 Halton points of dimension 49
and 50 XFGSOPLowDiscrepancy.xpl

However by using the Brownian Bridge path construction method we can limit
the e¬ect of the high dimensional components on a simulated underlying path
and the corresponding path variable for the most common path dependent
options, Moroko¬ (1996). This method start with an empty path with known
start value and calculates at each step the underlying value for a time point
with maximum time distance to all other time points with known underlying
value until the whole path is computed. Experimental results show that we
can still get a faster convergence of the QMC simulation for options up to 50
time points if we apply this path construction method.

16.3 Pricing options with simulation techniques -
a guideline
In this section we would like to give a short guideline how to price exotic
options with Monte Carlo and Quasi Monte Carlo simulation techniques
within the framework described above. Furthermore we give some indications
about the limits of these techniques.
362 16 Simulation based Option Pricing

16.3.1 Construction of the payo¬ function

As a ¬rst step we have to de¬ne the payo¬ function corresponding to our
option product. Within the methods de¬ned in the quantlib finance we have
to consider three di¬erent cases.

One underlying + path independent
In this case the payo¬ function is called by the pricing routine with the sim-
ulated underlying value at maturity as the single argument. It calculates the
corresponding payo¬ and returns this value. We have de¬ned the payo¬ func-
tion for a put option with strike price 100 as an example for a one dimensional
payo¬ function.

Several underlying + path independent
For options whose payo¬ depends on the underlying values of several assets at
maturity, we have to de¬ne a payo¬ function on the vector of the underlying
values at maturity. An example for such an option is an exchange option that
permits to swap a de¬ned share with the best performing share in a basket. Its
payo¬ function is given by:

CT ((ST , . . . , ST )) = max{0, ±i (ST ’ K i ) + 55 ’ ST |i = 1, .., 5}
1 5 i 3

One underlying + path dependent
The third category of option types that are captured are path dependent op-
tions on one underlying. The payo¬ function of these options depends on the
underlying values at several ¬xed time points during the lifetime of the option.
Payo¬ functions for these contracts are called with a vector of underlying val-
ues whose ith element is the underlying value at the time ti which has to be
speci¬ed in the model.

16.3.2 Integration of the payo¬ function in the simulation

After de¬ning the payo¬ function in XploRe we can start to calculate a price
estimate with the help of the appropriate simulation routine. In the one di-
mensional case we just have to call
16.3 Pricing options with simulation techniques - a guideline 363

erg = BlackScholesPathIndependent1D (s0,r,vola,dt,opt,
MC estimation of the option price for a path independent option.
erg = BlackScholesPathIndependent1DQMC (s0,r,vola,dt,opt,
QMC estimation of the option price for a path independent

to get a price estimate and for the Monte Carlo case an empirical standard
deviation with respect to a start price of s0, a continuous risk free interest
rate of r, a volatility vola, a time to maturity of dt years, the payo¬ function
opt, sample size itr and the random/low-discrepancy generator with number
gen. Table 16.1 shows the random number generators and table 16.2 the low-
discrepancy generators that can be used. An application of these routines for
a Put option can be found in XFGSOP1DPut.xpl.
Pricing path-dependent options is only slightly more complicated. Here we
have to de¬ne the vector of time points for which underlying prices have to
be generated. This vector replaces the time to maturity used to price path
independent options. Then we can apply one of the following methods to
compute a price estimate for the path dependent option

erg = BlackScholesPathDependent1D (s0,r,vola,times,opt,
MC estimation of the option price for path-dependent options.
erg = BlackScholesPathDependent1DQMC (s0,r,vola,times,opt,
QMC estimation of the option price for path-dependent options,

with respect to the start price s0, the continuous risk free interest rate r, the
volatility vola, the time scheme times, the payo¬ function opt, sample size
itr and the random/low-discrepancy generator with number gen, as given in
Tables 16.1 and 16.2. Using the above quantlets, we calculate the price of an
Asian call option in XFGSOP1DAsian.xpl.
364 16 Simulation based Option Pricing

In the case of multidimensional options we have to de¬ne a start price vector
and a covariance matrix instead of a single underlying price and volatility value.
Then we can call one of the multi-dimensional simulation routines:

erg = BlackScholesPathIndependentMD (s0,r,vola,dt,opt
MC estimation of the option price in the multidimensional Black
Scholes model

erg = BlackScholesPathIndependentMDQMC (s0,r,vola,dt,opt
QMC estimation of the option price in the multidimensional
Black Scholes model

with respect to the m dimensional start price vector s0, the continuous risk free
interest rate r, the m—m covariance matrix vola, the time to maturity dt, the
payo¬ function opt, the number of iterations itr and the generator number
gen according to the generators in Tables 16.1 and 16.2. Both quantlets are
illustrated in XFGSOPMD.xpl.
If in addition a dividend is paid during the time to maturity, we can use the
following two quantlets to calculate the option prices.

erg = BlackScholesPathIndependentMDDiv (s0,r,div,vola
MC estimation of the option price in the multidimensional Black
Scholes model
erg = BlackScholesPathIndependentMDDivQMC (s0,r,div,vola
QMC estimation of the option price in the multidimensional
Black Scholes model

The additional argument div is a m dimensional vector of the continuously paid
dividends. An application of these functions for our basket option is provided
in XFGSOPMDDiv.xpl.


. 12
( 14)