MCS Framework FEegs

Download as pdf or txt
Download as pdf or txt
You are on page 1of 10

Monte Carlo Simulation: IEOR E4703 Fall 2004

c 2004 by Martin Haugh


The Monte Carlo Framework, Examples from
Finance and Generating Correlated Random
Variables
1 The Monte Carlo Framework
Suppose we wish to estimate some quantity, = E[h(X)], where X = {X
1
, . . . , X
n
} is a random vector in R
n
,
h() is a function from R
n
to R, and E[|h(X)|] < .
Note that X could represent the values of a stochastic process at dierent points in time. For example, X
i
might be the price of a particular stock at time i and h() might be given by
h(X) =
X
1
+. . . +X
n
n
so then is the expected average value of the stock price. To estimate we use the following algorithm:
Monte Carlo Algorithm
for i = 1 to n
generate X
i
set h
i
= h(X
i
)
set

n
=
h
1
+h
2
+...+h
n
n
Note: If n is large, it may be necessary to keep track of

i
h
i
within the for loop, so that we dont have to
store each value of h
i
.
Question: Why is

a good estimator?
Answer: As we saw previously, there are two reasons:
1.

n
is unbiased. That is
E[

n
] =
E[

n
i
h
i
]
n
=
E[

n
i
h(X
i
)]
n
=
n
n
=
2.

n
is consistent. That is

n
wp 1 as n .
This follows from the Strong Law of Large Numbers (SLLN).
Remark 1 We can also estimate probabilities this way by representing them as expectations. In particular, if
= P(X A), then = E[I
A
(X)] where
I
A
(X) =
_
1 if X A
0 otherwise
The Monte Carlo Framework, Examples from Finance and Generating Correlated Random Variables 2
2 Examples from Finance
Example 1 (Portfolio Evaluation)
Consider two stocks, A and B, and let S
a
(t) and S
b
(t) be the time t prices of A and B, respectively. At time
t = 0, I buy n
a
units of A and n
b
units of B so my initial wealth is W
0
= n
a
S
a
(0) +n
b
S
b
(0). Suppose my
investment horizon is T years after which my terminal wealth, W
T
, is given by
W
T
= n
a
S
a
(T) +n
b
S
b
(T).
Note that this means that I do not trade in [0, T]. Assume S
a
GBM(
a
,
a
), S
b
GBM(
b
,
b
), and that
S
a
and S
b
are independent. We would now like to estimate
P
_
W
T
W
0
.9
_
,
i.e., the probability that the value of my portfolio drops by more than 10%. Note that we may write
S
a
(T) = S
a
(0) exp
_
(
a

2
a
/2)T +
a
B
a
(T)
_
S
b
(T) = S
b
(0) exp
_
(
b

2
b
/2)T +
b
B
b
(T)
_
where B
a
and B
b
are independent SBMs.
Let L be the event that W
T
/W
0
.9 so that that the quantity of interest is := P(L) = E[I
L
]. The problem
of estimating therefore falls into our Monte Carlo framework. In this example, X = (S
a
(T), S
b
(T)) and
I
L
(X) =
_
_
_
1 if
n
a
S
a
(T)+n
b
S
b
(T)
n
a
S
a
(0)+n
b
S
b
(0)
0.9
0 otherwise
Lets assume the following parameter values:
T = .5 years,
a
= .15,
b
= .12,
a
= .2,
b
= .18, S
a
(0) = $100, S
b
(0) = $75 and n
a
= n
b
= 100.
This then implies W
0
= $17, 500.
We then have the following algorithm for estimating :
Monte Carlo Estimation of
for i = 1 to n
generate X
i
= (S
i
a
(T), S
i
b
(T))
compute I
L
(X
i
)
set

n
=
I
L
(X
1
)+...+I
L
(X
n
)
n
The Monte Carlo Framework, Examples from Finance and Generating Correlated Random Variables 3
Sample Matlab Code
> n=1000; T=0.5; na =100; nb=100;
> S0a=100; S0b=75; mua=.15; mub=.12; siga=.2; sigb=.18;
> W0 = na*S0a + nb*S0b;
> BT = sqrt(T)*randn(2,n);
> STa = S0a * exp((mua - (siga^2)/2)*T + siga* BT(1,:));
> STb = S0b * exp((mub - (sigb^2)/2)*T + sigb* BT(2,:));
> WT = na*STa + nb*STb;
> theta_n = sum(WT./W0 < .9) / n
2.1 Introduction to Security Pricing
We assume again that S
t
GBM(, ) so that
S
T
= S
0
e
(
2
/2)T+B
T
.
In addition, we will always assume that there exists a risk-free cash account so that if W
0
is invested in it at
t = 0, then it will be worth W
0
exp(rt) at time t. We therefore interpret r as the continuously compounded
interest rate. Suppose now that we would like to estimate the price of a security that pays h(X) at time T,
where X is a random quantity (variable, vector, etc.) possibly representing the stock price at dierent times in
[0, T]. The theory of asset pricing then implies that the time 0 value of this security is
h
0
= E
Q
[e
rT
h(X)]
where E
Q
[.] refers to expectation under the risk-neutral probability
1
measure.
Risk-Neutral Asset Pricing
In the GBM model for stock prices, using the risk-neutral probability measure is equivalent to assuming that
S
t
GBM(r, ).
Note that we are not saying that the true stock price process is a GBM(r, ). Instead, we are saying that for
the purposes of pricing securities, we pretend that the stock price process is a GBM(r, ).
Example 2 (Pricing a European Call Option)
Suppose we would like to estimate
2
C
0
, the price of a European call option with strike K and expiration T,
using simulation. Our risk-neutral pricing framework tells us that
C
0
= e
rT
E
_
max
_
0, S
0
e
(r
2
/2)T+B
T
K
__
and again, this falls into our simulation framework. We have the following algorithm.
1
The risk-neutral probability measure is often called the equivalent martingale measure or EMM.
2
C
0
can be calculated exactly using the Black-Scholes formula but we ignore this fact for now!
The Monte Carlo Framework, Examples from Finance and Generating Correlated Random Variables 4
Estimating the Black-Scholes Option Price
set sum = 0
for i = 1 to n
generate S
T
set sum = sum+ max (0, S
T
K)
set

C
0
= e
rT
sum/n
Example 3 ( Pricing Asian Options)
Let the time T payo of the Asian option be
h(X) = max
_
0,

m
i=1
SiT
m
m
K
_
so that X =
_
ST
m
, S2T
m
, . . . , S
T
_
.
We can then use the following Monte Carlo algorithm to estimate C
0
= E
Q
[e
rT
h(X)].
Estimating the Asian Option Price
set sum = 0
for i = 1 to n
generate ST
m
, S2T
m
, . . . , S
T
set sum = sum+ max
_
0,

m
i=1
S
iT
m
m
K
_
set

C
0
= e
rT
sum/n
3 Generating Correlated Normal Random Variables and
Brownian Motions
In the portfolio evaluation example of Lecture 4 we assumed that the two stock price returns were independent.
Of course this assumption is too strong since in practice, stock returns often exhibit a high degree of correlation.
We therefore need to be able to simulate correlated random returns. In the case of geometric Brownian motion,
(and other models based on Brownian motions) simulating correlated returns means simulating correlated
normal random variables. And instead of posing the problem in terms of only two stocks, we will pose the
problem more generally in terms of n stocks. First, we will briey review correlation and related concepts.
3.1 Review of Covariance and Correlation
Let X
1
and X
2
be two random variables. Then the covariance of X
1
and X
2
is dened to be
Cov(X
1
, X
2
) := E[X
1
X
2
] E[X
1
]E[X
2
]
The Monte Carlo Framework, Examples from Finance and Generating Correlated Random Variables 5
and the correlation of X
1
and X
2
is then dened to be
Corr(X
1
, X
2
) = (X
1
, X
2
) =
Cov(X
1
, X
2
)
_
Var(X
1
)Var(X
2
)
.
If X
1
and X
2
are independent, then = 0, though the converse is not true in general. It can be shown that
1 1.
Suppose now that X = (X
1
, . . . , X
n
) is a random vector. Then , the covariance matrix of X, is the (n n)
matrix that has (i, j)
th
element given by
i,j
:= Cov(X
i
, X
j
).
Properties of the Covariance Matrix
1. It is symmetric so that
T
=
2. The diagonal elements satisfy
i,i
0
3. It is positive semi-denite so that x
T
x 0 for all x R
n
.
We will now see how to simulate correlated normal random variables.
3.2 Generating Correlated Normal Random Variables
The problem then is to generate X = (X
1
, . . . , X
n
) where X MN(0, ). Note that it is then easy to handle
the case where E[X] = 0. By way of motivation, suppose Z
i
N(0, 1) and IID for i = 1, . . . , n. Then
c
1
Z
1
+. . . c
n
Z
n
N(0,
2
)
where
2
= c
2
1
+. . . +c
2
n
. That is, a linear combination of normal random variables is again normal.
More generally, let C be a (n m) matrix and let Z = (Z
1
Z
2
. . . Z
n
)
T
. Then
C
T
Z MN(0, C
T
C)
so our problem clearly reduces to nding C such that
C
T
C = .
Question: Why is this true?
Finding such a matrix, C, requires us to compute the Cholesky decomposition of .
3.3 The Cholesky Decomposition of a Symmetric Positive-Denite Matrix
A well known fact from linear algebra is that any symmetric positive-denite matrix, M, may be written as
M = U
T
DU
where U is an upper triangular matrix and D is a diagonal matrix with positive diagonal elements.
Since our variance-covariance matrix, , is symmetric positive-denite, we can therefore write
= U
T
DU
= (U
T

D)(

DU)
= (

DU)
T
(

DU).
The matrix C =

DU therefore satises C
T
C = . It is called the Cholesky Decomposition of .
The Monte Carlo Framework, Examples from Finance and Generating Correlated Random Variables 6
3.3.1 Cholesky Decomposition in Matlab
It is easy to compute the Cholesky decomposition of a symmetric positive-denite matrix in Matlab using the
chol command. This means it is also easy to simulate multivariate normal random vectors as well.
As before, let be an (n n) variance-covariance matrix and let C be its Cholesky decomposition. If
X MN(0, ) then we can generate random samples of X in Matlab as follows:
Sample Matlab Code
>> Sigma = [1.0 0.5 0.5;
0.5 2.0 0.3;
0.5 0.3 1.5];
>> C=chol(Sigma);
>> Z=randn(3,1000000);
>> X=C*Z;
>> cov(X)
ans =
0.9972 0.4969 0.4988
0.4969 1.9999 0.2998
0.4988 0.2998 1.4971
Remark 2 We must be very careful to premultiply Z by C
T
and not C.
Example 4 (A Faulty )
We must ensure that is a genuine variance-covariance matrix. Consider the following Matlab code.
Matlab Code
>> Sigma=[0.5 0.9 0.4;
0.9 0.7 0.9;
0.4 0.9 0.9];
>> C=chol(Sigma)
Question: What is the problem here?
More formally, we have the following algorithm for generating multivariate random vectors, X.
Generating Correlated Normal Random Variables
generate Z MN(0, I)
/ Now compute the Cholesky Decomposition /
compute C such that C
T
C =
set X = C
T
Z
The Monte Carlo Framework, Examples from Finance and Generating Correlated Random Variables 7
3.4 Generating Correlated Brownian Motions
Generating correlated Brownian motions is, of course, simply a matter of generating correlated normal random
variables.
Denition 1 We say B
a
t
and B
b
t
are correlated SBMs with correlation coecient if E[B
a
t
B
b
t
] = t.
For two such SBMs, we then have
Corr(B
a
t
, B
b
t
) =
Cov(B
a
t
, B
b
t
)
_
Var(B
a
t
)Var(B
b
t
)
=
E[B
a
t
B
b
t
] E[B
a
t
]E[B
b
t
]
t
= .
Now suppose S
a
t
and S
b
t
are stock prices that follow GBMs such that
Corr(B
a
t
, B
b
t
) =
where B
a
t
and B
b
t
are the standard SBMs driving S
a
t
and S
b
t
, respectively. Let r
a
and r
b
be the continuously
compounded returns of S
a
t
and S
b
t
, respectively, between times t and t +s. Then it is easy to see that
Corr(r
a
, r
b
) = .
This means that when we refer to the correlation of stock returns, we are at the same time referring to the
correlation of the SBMs that are driving the stock prices. We will now see by example how to generate
correlated SBMs and, by extension, GBMs.
Example 5 (Portfolio Evaluation Revisited)
Recalling the notation of Example 1, we assume again that S
a
GBM(
a
,
a
) and S
b
GBM(
b
,
b
).
However, we no longer assume that S
a
t
and S
b
t
are independent. In particular, we assume that
Corr(B
a
t
, B
b
t
) =
where B
a
t
and B
b
t
are the SBMs driving A and B, respectively. As mentioned earlier, this implies that the
correlation between the return on A and the return on B is equal to . We would like to estimate
P
_
W
T
W
0
.9
_
,
i.e., the probability that the value of the portfolio drops by more than 10%. Again, let L be the event that
W
T
/W
0
.9 so that that the quantity of interest is := P(L) = E[I
L
], where X = (S
a
T
, S
b
T
) and
I
L
(X) =
_

_
1 if
n
a
S
a
T
+n
b
S
b
T
n
a
S
a
0
+n
b
S
b
0
0.9
0 otherwise.
Note that we can write
S
a
t+s
= S
a
t
exp
_
(
a

2
a
/2)s +

s V
a
_
S
b
t+s
= S
b
t
exp
_
(
b

2
b
/2)s +

s V
b
_
where (V
a
, V
b
) MN(0, ) with
=
_

2
a

a

b

2
b
_
.
So to generate one sample value of (S
a
t+s
, S
b
t+s
), we need to generate one sample value of V = (V
a
, V
b
). We do
this by rst generating Z MN(0, I), and then setting V = C
T
Z, where C is the Cholesky decomposition of
. To estimate we then set t = 0, s = T and generate n samples of I
L
(X). The following Matlab function
accomplishes this.
The Monte Carlo Framework, Examples from Finance and Generating Correlated Random Variables 8
A Matlab Function
function[theta] = portfolio_evaluation(mua,mub,siga,sigb,n,T,rho,S0a,S0b,na,nb);
% This function estimates the probability that wealth of the portfolio falls
%by more than 10% n is the number of simulated values of W_T that we use
W0 = na*S0a + nb*S0b;
Sigma = [siga^2 siga*sigb*rho;
siga*sigb*rho sigb^2];
B = randn(2,n);
C=chol(Sigma);
V = C * B;
STa = S0a * exp((mua - (siga^2)/2)*T + sqrt(T)*V(1,:));
STb = S0b * exp((mub - (sigb^2)/2)*T + sqrt(T)*V(2,:));
WT = na*STa + nb*STb;
theta = mean(WT/W0 < .9);
The function portfolio evaluation.m can now be executed by typing portfolio evaluation at the Matlab prompt.
3.4.1 Generating Correlated Log-Normal Random Variables
Let X be a multivariate lognormal random vector with mean

and variance -covariance matrix


3

. Then we
can write X = (e
Y
1
, . . . , e
Y
n
) where
Y := (Y
1
, . . . , Y
n
) MN(, ).
Suppose now that we want to generate a value of the vector X. We can do this as follows:
1. Solve for and in terms of

and

.
2. Generate a value of Y
3. Take X = exp(Y)
Step 1 is straightforward and only involves a few lines of algebra. (See Law and Kelton for more details.) In
particular, we now also know how to generate multivariate log-normal random vectors.
3
For a given positive semi-denite matrix,

, it should be noted that it is not necessarily the case that a multivariate


lognormal random vector exists with variance -covariance matrix,

.
The Monte Carlo Framework, Examples from Finance and Generating Correlated Random Variables 9
4 Simulating Correlated Random Variables in General
In general, there are two distinct approaches to modelling with correlated random variables. The rst approach
is one where the joint distribution of the random variables is fully specied. In the second approach, the joint
distribution is not fully specied. Instead, only the marginal distributions and correlations between the variables
are specied.
4.1 When the Joint Distribution is Fully Specied
Suppose we wish to generate a random vector X = (X
1
, . . . , X
n
) with joint CDF
F
x
1
,...,x
n
(x
1
, . . . , x
n
) = P(X
1
x
1
, . . . , X
n
x
n
)
Sometimes we can use the method of conditional distributions to generate X. For example, suppose n = 2.
Then
F
x
1
,x
2
(x
1
, x
2
) = P(X
1
x
1
, X
2
x
2
)
= P(X
1
x
1
) P(X
2
x
2
|X
1
x
1
)
= F
x
1
(x
1
)F
x
2
|x
1
(x
2
|x
1
)
So to generate (X
1
, X
2
), rst generate X
1
from F
x
1
() and then generate X
2
independently from F
x
2
|x
1
().
This of course will work in general for any value of n, assuming we can compute and simulate from all the
necessary conditional distributions. In practice the method is somewhat limited because we often nd that we
cannot compute and / or simulate from the distributions.
We do, however, have exibility with regards to the order in which we simulate the variables. Regardless, for
reasonable values of n this is often impractical. However, one very common situation where this method is
feasible is when we wish to simulate stochastic processes. In fact we use precisely this method when we simulate
Brownian motions and Poisson processes. More generally, we can use it to simulate other stochastic processes
including Markov processes and time series models among others.
Another very important method for generating correlated random variables is the Markov Chain Monte Carlo
(MCMC) method
4
.
4.2 When the Joint Distribution is not Fully Specied
Sometimes we do not wish to specify the full joint distribution of the random vector that we wish to simulate.
Instead, we may only specify the marginal distributions and the correlations between the variables. Of course
such a problem is then not fully specied since in general, there will be many dierent joint probability
distributions that have the same set of marginal distributions and correlations. Nevertheless, this situation often
arises in modelling situations when there is only enough data to estimate marginal distributions and correlations
between variables. In this section we will briey mention some of the issues and methods that are used to solve
such problems.
Again let X = (X
1
, . . . , X
n
) be a random vector that we wish to simulate.
Now, however, we do not specify the joint distribution F(x
1
, . . . , x
n
)
Instead, we specify the marginal distributions F
x
i
(x) and the covariance matrix R
4
See Ross for an introduction to MCMC.
The Monte Carlo Framework, Examples from Finance and Generating Correlated Random Variables 10
Note again that in general, the marginal distributions and R are not enough to uniquely specify the joint
distribution
Even now we should mention that potential diculties already exist. In particular, we might have a consistency
problem. That is, the particular R we have specied may not be consistent with the specied marginal
distributions. In this case, there is no joint CDF with the desired marginals and correlation structure.
Assuming that we do not have such a consistency problem, then one possible approach, based on the
multivariate normal distribution, is as follows.
1. Let (Z
1
, . . . , Z
n
) MN(0, ) where is a covariance matrix with 1s on the diagonal. Therefore
Z
i
N(0, 1) for i = 1, . . . , n.
2. Let () and () be the PDF and CDF, respectively, of a standard normal random variable. It can then
be seen that ((Z
1
), . . . , (Z
n
)) is a random vector where the marginal distribution of each (Z
i
) is
uniform. How would you show this?
3. Since the Z
i
s are correlated, the uniformly distributed (Z
i
)s will also be correlated. Let

denote the
variance-covariance matrix of the (Z
i
)s.
4. We now set X
i
= F
1
x
i
((Z
i
)). Then X
i
has the desired marginal distribution, F
x
i
(). Why?
5. Now since the (Z
i
)s are correlated, the X
i
s will also be correlated. Let

denote the
variance-covariance matrix of the X
i
s.
6. Recall that we want X to have marginal distributions F
x
i
() and covariance matrix R. We have satised
the rst condition and so to satisfy the second condition, we must have

= R. So we must choose the


original in such a way that

= R. (1)
7. In general, choosing appropriately is not trivial and requires numerical work. It is also true that there
does not always exist a such that (1) holds.

You might also like