IEOR E4602: Quantitative Risk Management Spring 2016

c 2016 by Martin Haugh

Monte-Carlo Methods for Risk Management

In these lecture notes we discuss Monte-Carlo (MC) techniques that are particularly useful in risk-management
applications. We focus on importance sampling and stratified sampling, both of which are variance reduction
techniques that can be very useful in estimating risk measures associated with rare-events. We consider several
applications of these techniques to risk management but there are many more. For example, Chapter 9 of
Glasserman’s Monte-Carlo Methods in Financial Engineering (2004) contains applications of importance and
stratified sampling to credit risk and the estimation of risk measures in both light- and heavy-tailed settings. We
also discuss the concept of credit valuation adjustment (CVA) which reflects the price adjustment that must be
made to derivative prices when we account for the possibility of counter-party defaults. Finally we briefly discuss
the use of Monte-Carlo for estimation of the Greeks.
We assume that readers are already familiar with Monte-Carlo simulation and know, in particular, how to
generate random variables and analyze simulation output. We also assume that readers have had some exposure
to variance reduction methods.

1 Variance Reduction and Importance Sampling

We begin with a seemingly contrived example but one that nonetheless gained some notoriety at the beginning
of the financial crisis in August 2007 when Goldman Sachs Asset Management’s (GSAM) Global Alpha fund
incurred steep losses. In explaining these losses, the CFO of Goldman Sachs claimed they had seen 25 standard
deviation moves several days in a row. Not surprisingly this claim resulted in considerable criticism.

Example 1 (Just How Unlucky is a 25 Standard Deviation (Negative) Return?)

Suppose we wish to estimate θ := P(X ≥ 25) = E[I{X≥25} ] where X ∼ N(0, 1). The usual Monte-Carlo
approach to this problem proceeds as follows:

1. Generate X1 , . . . , Xn IID N(0, 1)

2. Set Ij = I{Xj ≥25} for j = 1, . . . , n
Pn
3. Set θbn = j=1 Ij /n
√
4. Compute an approximate 95% CI as θbn ± 1.96 × σ
bn / n where σ
bn is the standard deviation of the Ij ’s.

For this problem, however, the usual approach would be completely inadequate since approximating θ to any
reasonable degree of accuracy would require n to be inordinately large. For example, we will soon see that on
average we would have to set n ≈ 3.26 × 10137 in order to obtain just one non-zero value of I. Clearly this is
impractical and a much smaller value of n would have to be used. Using a much smaller value of n, however,
would almost inevitably result in an estimate, θbn = 0, and an approximate confidence interval [L, U ] = [0, 0]! So
the naive Monte-Carlo approach does not work here.

Before proceeding further, it is not unreasonable to ask why such a problem would be important. After all, if
you want to estimate θ = P(X ≥ 25), isn’t it enough to know that θ is very close to 0? Put another way, do we
care whether θ = 10−10 or θ = 10−20 ? For many problems, this is a valid objection, as we may care just how
small θ as long as we know that it is indeed “small”. However, for many other other problems it is very
Monte-Carlo Methods for Risk Management 2

important to know θ to a much greater level of accuracy. For example, suppose we are designing a nuclear
power plant and we want to estimate the probability, θ, of a meltdown occurring sometime in the next 100
years. We would expect θ to be very small, even for a poorly designed power plant. However, this is not enough.
Should a meltdown occur, then clearly the consequences could be catastrophic and so we would like to know θ
to a very high degree of accuracy.
For another example, suppose we want to price a deep-out-of-the-money option using simulation. Then the price
of the option will be very small, perhaps lying between .1 cents and 10 cents. Clearly a bank is not going to
suffer if it misprices this option and sells it for .1 cents when the correct value is 10 cents. But what if the bank
sells 1 million of these options? And what if the bank makes similar trades several times a week? Then it
becomes very important to price the option correctly. A particularly rich source of examples can be found in
risk-management where we seek to estimate risk measures such as the VaRα or ESα of a given portfolio.

1.1 Introduction and Main Results

Note that these examples require estimating the probability of a rare event. Even though the events are rare,
they are very important because when they do occur their impact can be very significant. We now consider
importance sampling, a variance reduction technique that can be invaluable when estimating rare event
probabilities and expectations.
Suppose we wish to estimate θ = Ef [h(X)] where X has PDF f (or PMF, if X is a discrete random variable).
Let g be another PDF with the property that g(x) 6= 0 whenever f (x) 6= 0. That is, g has the same support as
f . Then
Z  
f (x) h(X)f (X)
θ = Ef [h(X)] = h(x) g(x) dx = Eg
g(x) g(X)
where Eg [ · ] denotes an expectation with respect to the density g. This has very important implications for
estimating
P θ. The original simulation method is to generate n samples of X from the density, f , and set
θbn = h(Xj )/n. An alternative method, however, is to generate n values of X from the density, g, and set
n
X h(Xj )f (Xj )
θbn,is = .
j=1
ng(Xj )

θbn,is is then an importance sampling estimator of θ. We often define h∗ (X) := h(X)f (X)/g(X) so that
θ = Eg [h∗ (X)]. We refer to f and g as the original and importance sampling densities, respectively. We also
refer to f /g as the likelihood ratio.

Example 2 (Revisiting the 25 Standard Deviation Example)

Consider again the problem where we want to estimate θ = P(X ≥ 25) = E[I{X≥25} ] when X ∼ N(0, 1). We
may then write
Z ∞
1 x2
θ = E[I{X≥25} ] = I{X≥25} √ e− 2 dx
−∞ 2π
x2
 
Z ∞ √1 e− 2
 √1 e− 2 dx
(x−µ)2
2π
= I{X≥25}  (x−µ)2
−∞ √1 e− 2 2π
2π
h 2
i
= Eµ I{X≥25} e−µX+µ /2
2
where now X ∼ N(µ, 1) and e−µX+µ /2 is the likelihood ratio. If we set µ = 25, for example, and use n = 10
million samples, then we find an approximate 95% confidence interval for θ is given by [3.053, 3.074] × 10−138 .

We can of course also estimate expectations using importance sampling.

Monte-Carlo Methods for Risk Management 3
2
Example 3 Suppose we wish to estimate θ = E[X 4 eX /4 I{X≥2} ] where X ∼ N(0, 1). Then the same
2 2
argument as before implies that we may also write θ = Eµ [X 4 eX /4 e−µX+µ /2 I{X≥2} ] where now X ∼ N(µ, 1).

The General Formulation

Let X = (X1 , . . . , Xn ) be a random vector with joint PDF f (x1 , . . . , xn ) and suppose we wish to estimate
θ = Ef [h(X)]. Let g(x1 , . . . , xn ) be another PDF such that g(x) 6= 0 whenever f (x) 6= 0. Then we easily obtain

θ = Ef [h(X)]
= Eg [h∗ (X)]

where h∗ (X) := h(X)f (X)/g(X). Again we have two methods for estimating θ: the original method where we
simulate with respect to the density function, f , and the importance sampling method where we simulate with
respect to the density, g.

Pn 2


Example 4 Suppose we wish to estimate θ = P i=1 Xi ≥ 50 where the Xi ’s are IID N(0, 1). Then
θ = E[h(X)] where h(X) := I{P X 2 ≥50} and X := (X1 , . . . , Xn ). We could estimate θ using importance
i
sampling as follows.
2 2
e−x1 /2 e−xn /2 P
Z Z
θ = E[h(X)] = ... √ ... √ I dx1 . . . dxn
2π { Xi ≥50}
2
x1 xn 2π
2 2
! 2 2 2 2
e−x1 /2 e−xn /2 e−x1 /2σ e−xn /2σ
Z Z
n
= σ ... . . . −x2 /2σ2 √ ... √ I P X 2 ≥50} dx1 . . . dxn
2πσ 2 {
2 2
x1 xn e−x1 /2σ e n 2πσ 2 i

2 2 2 2
e−x1 /2σ e−xn /2σ
 
x2
Z Z
x2
1 (1−1/σ 2 ) n (1−1/σ 2 )
= σn ... e− 2 . . . e− 2 √ ... √ I P X 2 ≥50} dx1 . . . dxn
x1 xn 2πσ 2 2πσ 2 { i

  2 2
 
X1 Xn
(1− σ12 ) (1− σ12 )
= Eg σ n e− 2 . . . e− 2 I{P X 2 ≥50}
i

where Eg [.] denotes expectation under a multivariate normal distribution where the X ∼ MVN(0, σ 2 In ).
Thus far we have not addressed the issue of how to choose a good sampling density, g, so that we obtain a
variance reduction when we sample from g instead of f . We will now address this question.

Obtaining a Variance Reduction

As before, suppose we wish to estimate θ = Ef [h(X)] where X is a random vector with joint PDF, f . Without
loss of generality we will assume that h(X) ≥ 0. Now let g be another density with support equal to that of f .
Then we know
θ = Ef [h(X)] = Eg [h∗ (X)]
and this gives rise to two estimators:
1. h(X) where X ∼ f and
2. h∗ (X) where X ∼ g
The variance of the importance sampling estimator is given by
Z
Varg (h∗ (X)) = h∗ (x)2 g(x) dx − θ2

h(x)2 f (x)
Z
= f (x) dx − θ2
g(x)
Monte-Carlo Methods for Risk Management 4

while the variance of the original estimator is given by Varf (h(X)) = h(x)2 f (x) dx − θ2 . So the reduction
R

in variance is then given by

Z  
f (x)
Varf (h(X)) − Varg (h∗ (X)) = h(x)2 1 − f (x) dx. (1)
g(x)
In order to achieve a variance reduction, the integral in (1) should be positive. For this to happen, we would like
1. f (x)/g(x) > 1 when h(x)2 f (x) is small and
2. f (x)/g(x) < 1 when h(x)2 f (x) is large.

Now the important part of the density, f , could plausibly be defined to be that region, A say, in the support of
f where h(x)f (x) is large. But by the above observation, we would like to choose g so that f (x)/g(x) is small
whenever x is in A. That is, we would like a density, g, that puts more weight on A: hence the name
importance sampling. Note that when h involves a rare event so that h(x) = 0 over “most” of the state space,
it can then be particularly valuable to choose g so that we sample often from that part of the state space where
h(x) 6= 0. This is why importance sampling is most useful for simulating rare events. Further guidance on how
to choose g is obtained from the following observation.
As we are free to choose g, let’s suppose we choose g(x) = h(x)f (x)/θ. Then it is easy to see that
Varg (h∗ (X)) = θ2 − θ2 = 0
so that we have a zero variance estimator! This means that if we sample with respect to this particular choice of
g, then we would only need one sample and this sample would equal θ with probability1 one. Of course this is
not feasible in practice since we don’t know θ and therefore don’t know g either. However, all is not lost and this
observation can often guide us towards excellent choices of g that lead to extremely large variance reductions.

How to Choose a Good Sampling Distribution

We saw above that if we could choose g(x) = h(x)f (x)/θ, then we would obtain the best possible estimator of
θ, that is, one that has zero variance. In general, we cannot do this, but it does suggest that if we could choose
g(·) so that it is similar to h(·)f (·), then we might reasonably expect to obtain a large variance reduction.
What does the phrase “similar” mean? One obvious thing to do would be to choose g(·) so that it has a similar
shape to h(·)f (·). In particular, we could try to choose g so that g(x) and h(x)f (x) both take on their
maximum values at the same value, x∗ , say. When we choose g this way, we say that we are using the
maximum principle. Of course this only partially defines g since there are infinitely many density functions
that could take their maximum value at x∗ . Nevertheless, this is often enough to obtain a significant variance
reduction and in practice, we often take g to be from the same family of distributions as f . For example, if f is
multivariate normal, then we might also take g to be multivariate normal but with a different mean vector and /
or variance-covariance matrix.2

2
Example 5 Returning to Example 3, recall that we wished to estimate θ = E[h(X)] = E[X 4 eX /4 I{X≥2} ]
where X ∼ N(0, 1). If we sample from a PDF, g, that is also normal with variance 1 but mean µ, then we know
that g takes it maximum value at x = µ. Therefore, a good choice of µ might be
2 √
µ = arg max h(x)f (x) = arg max x4 e−x /4 = 8.
x x≥2
√
4 X 2 /4 − 8X+4
√
Then θ = Eg [h∗ (X)] = Eg [X e e I{X≥2} ] where g(·) denotes the N( 8, 1) PDF.

1 With this choice of g we have h∗ (x) = h(x)f (x)/g(x) = θ. Note also that this choice of g is valid since g(x) dx = 1 and
R
we have assumed h is non-negative.
2 We note that it is not necessary that f and g come from the same family of distributions. In fact sometimes it is necessary

to choose g from a different family of distributions. This might occur, for example, if it is difficult or inefficient to simulate
from the family of distributions to which f belongs. In that case, our reason for using importance sampling in the first place is
so that we can simulate from an ‘easier’ distribution, g.
Monte-Carlo Methods for Risk Management 5

Example 6 (Pricing an Asian Call Option)

For the purpose of option pricing, we assume that St ∼ GBM (r, σ 2 ), where St is the price of the stock at time
t and r is the risk-free interest rate. Suppose now that we want to price an Asian call option whose payoff at
time t is given by  Pm 
i=1 SiT /m
h(S) := max 0, −K (2)
m
where S := {SiT /m : i = 1, . . . , m}, T is the expiration date and K is the strike price. The price of this option
is then given by Ca = E[e−rT h(S)]. Now we can write
2 iT
√T
SiT /m = S0 e(r−σ /2) m +σ m (X1 +...+Xi )

where the Xi ’s are IID N(0, 1). This means that if f is the joint PDF of X = (X1 , . . . , Xm ), then (with a mild
abuse of notation) we may write
Ca = Ef [h(X1 , . . . , Xn )].
Now if K is very large relative to S0 so that the option is deep out-of-the-money then pricing the option using
simulation amounts to performing a rare event simulation. As a result, estimating Ca using importance sampling
will often result in a very large variance reduction. In order to apply importance sampling, we need to choose the
sampling density, g. For this, we could take g to be the multivariate normal PDF with variance-covariance
matrix equal to the identity, Im , and mean vector, µ∗ . As before, a good possible value of µ∗ might be
µ∗ = arg maxx h(x)f (x) which can be found using numerical methods.

Potential Problems with the Maximum Principle

Sometimes applying the maximum principle to choose g will be difficult. For example, it may be the case that
there are multiple or even infinitely many solutions to µ∗ = arg maxx h(x)f (x). Even when there is a unique
solution, it may be the case that finding it is very difficult. In such circumstances, an alternative method for
choosing g is to scale f . We will demonstrate this by example.

Example 7 (Using Scaling to Select g)

Assume in Example 4 that n = 2 so that θ = P X12 + X22 ≥ 50 = E[I{X12 +X22 ≥50} ] where X1 , X2 are IID
N(0, 1). Then
2 2
e−(x1 +x2 )/2
h(x)f (x) = I{x21 +x22 ≥50}
2π
so that h(x)f (x) = 0 inside the circle x21 + x22 ≤ 50 and h(x)f (x) takes on its maximum value at every point
on the circle x21 + x22 = 50. As a result, it is not possible to apply the maximum principle. Before choosing a
sampling density, g, recall that we would like g to put more weight on those parts of the sample space where
h(x)f (x) is large. One way to achieve this is by scaling the density of X = (X1 , X2 ) so that X is more
dispersed. For example, we could take g to be multivariate normal with mean vector 0 and variance-covariance
matrix  2 
σ 0
Σg =
0 σ2
where σ 2 > 1. Note that this simply means that under g, X1 and X2 are IID N(0, σ 2 ). Furthermore, when
σ 2 > 1, then more probability mass is given to the region X12 + X22 ≥ 50 as desired. We could choose the value
of σ using heuristic methods. One method would be to choose σ so that Eg [X12 + X22 ] = 50 which in this case
would imply that σ = 5. Why? We then have

X12 X22
   
2 2 2
θ = E[I{X12 +X22 ≥50} ] = Eg σ exp − (1 − 1/σ ) − (1 − 1/σ ) I{X12 +X22 ≥50} .
2 2
Pn
For the more general case where n > 2, we could proceed by again choosing σ so that Eg [ i=1 Xi2 ] = 50.
Monte-Carlo Methods for Risk Management 6

Difficulties with Importance Sampling

The most difficult aspect to importance sampling is in choosing a good sampling density, g. In general, one
needs to be very careful for it is possible to choose g according to some good heuristic such as the maximum
principle, but to then find that g results in a variance increase. In fact it is possible to choose a g that results in
an importance sampling estimator that has an infinite variance! This situation would typically occur when g puts
too little weight relative to f on the tails of the distribution. In more sophisticated applications of importance
sampling it is desirable to have (or prove) some guarantee that the importance sampling variance will be finite.
As an example of such a guarantee see Example 9 below.

1.2 Using Tilted Densities to Obtain a Good Sampling Distribution

Suppose f is light-tailed so that it has a moment generating function (MGF). Then a very useful way of
generating the sampling density, g, from the original density, f , is to use the MGF of f . We use Mx (t) to
denote the MGF and it is defined by
Mx (t) = Ef [etX ].
Then a tilted density of f is given by
etx f (x)
ft (x) =
Mx (t)
for −∞ < t < ∞. The tilted densities are useful since a random variable with density ft (·) tends to be larger
than one with density f when t > 0, and smaller when t < 0. This means, for example, that if we want to
sample more often from the region where X tends to be large, we might want to use a tilted density with t > 0
as our sampling density. Similarly, if we want to sample more often from the region where X tends to be small,
then we might use a tilted density with t < 0.

Example 8 Suppose X is an exponential random variable with mean 1/λ. Then f (x) = λe−λx for x ≥ 0, and
it is easy to see that ft (x) = Ce−(λ−t)x where C is the constant that makes the density integrate to 1.

Example 9 (The Probability that a Sum of Random Variables Will Exceed SomePValue)
n
Suppose X1 , . . . , Xn are independent random variables, where Xi has density fi (·). Let Sn := i=1 Xi and
suppose we want to estimate θ := P(Sn ≥ a) for some constant, a. If a is large so that we are dealing with a
rare event we should use importance sampling to estimate θ. Since Sn is large when the Xi ’s are large it makes
sense to sample each Xi from its tilted density function, fi,t (·) for some value of t > 0. We may then write
" n
# " n
! #
Y fi (Xi ) Y
θ = E[I{Sn ≥a} ] = Et I{Sn ≥a} = Et I{Sn ≥a} Mi (t) e−tSn
i=1
fi,t (Xi ) i=1

where Et [.] denotes expectation with respect to the Xi ’s Q

under the tilted densities, fi,t (·), and Mi (t) is the
n
moment generating function of Xi . If we write M (t) := i=1 Mi (t), then it is easy to see that the importance
sampling estimator, θbn,i , satisfies
0 ≤ θbn,i ≤ M (t)e−ta . (3)
Therefore a good choice of t would be that value that minimizes the bound in (3). We can minimize this by
minimizing log(M (t)e−ta ) = log(M (t)) − ta. It is straightforward to check that the minimizing value of t
satisfies µt = a where µt := Et [Sn ]. This can easily be found numerically.

Applications From Insurance: Estimating Ruin Probabilities

Continuing on from Example 9, if we define the stopping time τa := min{n ≥ 0 : Sn ≥ a}, then P (τa < ∞)
is the probability that Sn ever exceeds a. If E[X1 ] > 0 and the Xi ’s are IID with MGF, MX (t), then this
Monte-Carlo Methods for Risk Management 7

probability equals one. The case of interest is then when E[X1 ] ≤ 0. A similar argument to that of Example 9
yields
"∞ # ∞
X X  
θ = E[I{τa <∞} ] = E 1{τa =n} = E 1{τa =n}
n=1 n=1
X∞
n
Et 1{τa =n} (MX (t)) e−tSn
 
=
n=1
X∞
τ
Et 1{τa =n} (MX (t)) a e−tSτa
 
=
n=1
h i
= Et I{τa <∞} e−tSτa +τa ψ(t)

where ψ(t) := log(MX (t)) is the

 cumulant generating function. Note that if Et [X1 ] > 0 then τa < ∞ almost
surely and so we obtain θ = Et e−tSτa +τa ψ(t) . In fact, this is an example where importance sampling can be


used to ensure that the simulation stops almost surely. It is possible to choose a good value of t based on the
cumulant generating function.
Note that this example has direct applications to the estimation of ruin probabilities in the context of insurance
risk. For example, suppose Xi := Yi − cTi where Yi is the size of the ith claim, Ti is the inter-arrival time
between claims, c is the premium received per unit time and a is the initial reserve. Then θ is the probability
that the insurance company ever goes bankrupt. Only in very simple models is it possible to calculate θ
analytically. In general, Monte-Carlo approaches are required.

1.3 Estimating Conditional Expectations

Importance sampling can also be very useful for computing conditional expectations when the event being
conditioned upon is a rare event. For example, suppose we wish to estimate θ = E[h(X)|X ∈ A] where A is a
rare event and X is a random vector with PDF, f . Then the density of X, given that X ∈ A, is
f (x)
f (x|x ∈ A) = for x ∈ A
P(X ∈ A)
so
E[h(X)I{X∈A} ]
θ= . (4)
E[I{X∈A} ]
Now since A is a rare event we would be much better off if we could simulate using a sampling density, g, that
makes A more likely to occur. Then, as usual, we would have
Eg [h(X)I{X∈A} f (X)/g(X)]
θ= .
Eg [I{X∈A} f (X)/g(X)]

So to estimate θ using importance sampling, we would generate X1 , . . . , Xn with density g(·), and set
Pn
h(Xi )I{Xi ∈A} f (Xi )/g(Xi )
θbn,i = i=1
Pn .
i=1 I{Xi ∈A} f (Xi )/g(Xi )

In contrast to our usual estimators, θbn,i is no longer an average of n IID random variables but instead, it is the
ratio of two such averages. This has implications for computing approximate confidence intervals for θ. In
particular, confidence intervals should now be estimated using bootstrapping techniques. An obvious application
of this methodology in risk management is in the estimation of quantities similar to ES or CVaR.

Exercise 1 How does the practical problem of estimating the α-CVaR of a loss distribution, i.e. estimating
θ := E[L | L ≥ VaRα ], generally differ from the problem of estimating θ in (4)?
Monte-Carlo Methods for Risk Management 8

2 An Application of Importance Sampling to Credit Risk

Pm
We consider3 a portfolio loss of the form L = i=1 ei Yi where ei is the deterministic and positive exposure to
the ith credit and Yi is the default indicator with corresponding default probability, pi . We assume that Y
follows a Bernoulli mixture model which we now define.
Definition 1 Let p < m and let Ψ = (Ψ1 , . . . , Ψp )> be a p-dimensional random vector. Then we say the
random vector Y = (Y1 , . . . , Ym )> follows a Bernoulli mixture model with factor vector Ψ if there are
functions pi : Rp → [0, 1], 1 ≤ i ≤ m, such that conditional on Ψ the components of Y are independent
Bernoulli random variables satisfying P (Yi = 1 | Ψ = ψ) = pi (ψ).
We are interested in the problem of estimating θ := P (L ≥ c) where c is substantially larger than E[L]. Note
that a good importance sampling distribution for θ should also yield a good importance sampling distribution for
computing risk measures associated with the α-tail of the loss distribution where qα (L) ≈ c. We begin with the
case where the default indicators are independent.

2.1 Independent Default Indicators

We define Ω to be the state space of Y so that Ω = {0, 1}m . Then
m
Y
P ({y}) = pyi i (1 − pi )1−yi , y∈Ω
i=1

so that
m
Y m
Y
ML (t) = E[etL ] = E[etei Yi ] = pi etei + 1 − pi .


i=1 i=1
Let Qt be the corresponding tilted probability measure so that
Pm m
et i=1 ei y i Y etei yi
Qt ({y}) = P ({y}) = pyi (1 − pi )1−yi
tei + 1 − p ) i
ML (t) i=1
(pi e i
m
Y yi
= qt,i (1 − qt,i )1−yi
i=1

where qt,i := pi e /(pi e + 1 − pi ) is the Qt probability of the ith credit defaulting. Note that the default
tei tei

Pmqt,i → 1 as t → ∞ and qt,i → 0 as

indicators remain independent Bernoulli random variables under Qt . Since
t → −∞ it is clear that we can shift the mean of L to any value in (0, i=1 ei ). The same argument
Pm that was
used at the end of Example 9 suggests that we take t equal to that value that solves Et [L] = i=1 qi,t ei = c.
This value can be found easily using numerical methods.

2.2 Dependent Default Indicators

Suppose now that there is a p-dimensional factor vector, Ψ, such that the default indicators are independent
with default probabilities pi (ψ) conditional on Ψ = ψ. Suppose also that Ψ ∼ MVNp (0, Σ). The Monte-Carlo
scheme for estimating θ is to first simulate Ψ and to then simulate Y conditional on Ψ. We can apply
importance sampling to the second step using our discussion of independent default indicators above. However,
we can also apply importance sampling to the first step. A natural way to do this is to simulate Ψ from the
MVNp (µ, Σ) distribution for some µ ∈ Rp . The corresponding likelihood ratio, rµ (Ψ) say, (f /g in our earlier
notation) is given by the ratio of the two multivariate normal densities and satisfies
exp − 12 Ψ> Σ−1 Ψ


1
rµ (Ψ) =
= exp(−µ> Σ−1 Ψ + µ> Σ−1 µ).
exp − 12 (Ψ − µ)> Σ−1 (Ψ − µ) 2
3 This application is described in Section 8.5 of MFE, but is based on the original paper of Glasserman and Li (2005).
Monte-Carlo Methods for Risk Management 9

How Do We Choose µ?
Recall that the quantity of interest is θ := P (L ≥ c) = E[P (L ≥ c | Ψ)]. We know from our earlier discussion of
importance sampling that we would like to choose the importance sampling density, g ∗ (Ψ) say, so that
1
g ∗ (Ψ) ∝ P (L ≥ c | Ψ) exp(− Ψ> Σ−1 Ψ). (5)
2
Of course this is not possible since we do not know P (L ≥ c | Ψ), the very quantity that we wish to estimate.
The maximum principle applied to the MVNp (µ, Σ) distribution would then suggest taking µ equal to the value
of Ψ which maximizes the right-hand-side of (5). Again it is not possible to solve this problem exactly as we do
not know P (L ≥ c | Ψ) but numerical methods can be used to find good approximate solutions. See
Glasserman and Li (2005) for further details.

The Importance Sampling Algorithm for Estimating θ = P (L ≥ c)

1. Generate Ψ1 , . . . , Ψn independently from the MVNp (µ, Σ) distribution
2. For each Ψi estimate P (L ≥ c | Ψ = Ψi ) using the importance sampling distribution that we described in
our discussion of independent default indicators. Let θbnIS1 (Ψi ) be the corresponding estimator based on n1
samples.
3. The full importance sampling estimator is then given by
n
1X
θbnIS = rµ (Ψi ) θbnIS1 (Ψi ).
n i=1

3 Variance Reduction and Stratified Sampling

Consider a game show where contestants first pick a ball at random from an urn and then receive a payoff, Y.
The payoff is random and depends on the color of the selected ball so that if the color is c then Y is drawn from
the PDF, fc . The urn contains red, green, blue and yellow balls, and each of the four colors is equally likely to
be chosen. The producer of the game show would like to know how much a contestant will win on average when
he plays the game. To answer this question, she decides to simulate the payoffs of n contestants and take their
average payoff as her estimate. The payoff, Y , of each contestant is simulated as follows:
1. Simulate a random variable, I, where I is equally likely to take any of the four values r, g, b and y

2. Simulate Y from the density fI (y).

The average payoff, θ := E[Y ], is then estimated by
Pn
j=1 Yj
θbn := .
n
Now suppose n = 1000, and that a red ball was chosen 246 times, a green ball 270 times, a blue ball 226 times
and a yellow ball 258 times.
Question: Would this influence your confidence in θbn ? What if fg tended to produce very high payoffs and fb
tended to produce very low payoffs?
Is there anything that we could have done to avoid this type of problem occurring? The answer is yes. We know
that each ball color should be selected 1/4 of the time so we could force this to be true by conducting four
separate simulations, one each to estimate E[X|I = c] for c = r, g, b, y. Note that
1 1 1 1
E[Y ] = E[E[Y |I]] = E[Y |I = r] + E[Y |I = g] + E[Y |I = b] + E[Y |I = y]
4 4 4 4
Monte-Carlo Methods for Risk Management 10

so that an unbiased estimator of θ is obtained by setting

1 1 1 1
θbst,n := θbr,nr + θbg,ng + θbb,nb + θby,ny (6)
4 4 4 4

where θc := E[Y |I = c] for c = r, g, b, y.4 How does the variance of θbst,n compare with the variance of θbn , the
original raw simulation estimator? To answer this question, assume for now that nc = n/4 for each c, and that
Yc is a sample from the density, fc . Then a fair comparison of Var(θbn ) with Var(θbst,n ) should compare

Var(Y1 + Y2 + Y3 + Y4 ) with Var(Yr + Yg + Yb + Yy ) (7)

where Y1 , Y2 , Y3 and Y4 are IID samples from the original simulation algorithm, i.e. where we first select the
ball randomly and then receive the payoff, and the Yc ’s are independent with density fc (·), for c = r, g, b, y.
Now recall the conditional variance formula which states

Var(Y ) = E[Var(Y |I)] + Var(E[Y |I]). (8)

Each term in the right-hand-side of (8) is non-negative so this implies

Var(Y ) ≥ E[Var(Y |I)]

1 1 1 1
= Var(Y |I = r) + Var(Y |I = g) + Var(Y |I = b) + Var(Y |I = y)
4 4 4 4
Var(Yr + Yg + Yb + Yy )
=
4
which implies
Var(Y1 + Y2 + Y3 + Y4 ) = 4Var(Y ) ≥ Var(Yr + Yg + Yb + Yy ). (9)
As a result, we may conclude that using θbst,n instead of θbn leads to a variance reduction. This variance
reduction will be substantial if I accounts for a large fraction of the variance of Y . Note also that the
computational requirements for computing θbst,n are similar5 to those required for computing θbn . We call θbst,n a
stratified sampling estimator of θ, and we say that I is the stratification variable.

3.1 The Stratified Sampling Algorithm

We will now formally describe the stratified sampling algorithm. Suppose as usual that we wish to estimate
θ := E[Y ] where Y is a random variable. Let W be another random variable that satisfies the following two
conditions:
Condition 1: For any ∆ ⊆ R, P(W ∈ ∆) can be easily computed.
Condition 2: It is easy to generate (Y |W ∈ ∆), i.e., Y given W ∈ ∆.

In order to achieve a variance reduction it must also be the case that Y and W are dependent.
Pm Now divide R
into m non-overlapping subintervals, ∆1 , . . . , ∆m , such that pj := P(W ∈ ∆j ) > 0 and j=1 pj = 1. Note
that if W can take any value in R, then the first interval should be [−∞, a], while the final interval should be
[b, ∞] for some finite a and b.

We will use the following notation:

1. Let θj := E[Y |W ∈ ∆j ] and σj2 := Var(Y |W ∈ ∆j ).
2. We define the random variable I by setting I := j if W ∈ ∆j .
4θbc,n is an estimate of θc using nc samples. θbst,n is an estimate of θ using n samples, so it is implicitly assumed in (6) that
c
nr + ng + nb + ny = n.
5 For this example, the stratified estimator will actually require less work, but it is also possible in general for it to require

more work.
Monte-Carlo Methods for Risk Management 11

3. Let Y (j) denote a random variable with the same distribution as (Y |W ∈ ∆j ) ≡ (Y |I = j).
Our notation then implies θj = E[Y |I = j] = E[Y (j) ] and σj2 = Var(Y |I = j) = Var(Y (j) ). In particular we
have

θ = E[Y ] = E[E[Y |I]] = p1 E[Y |I = 1] + . . . + pm E[Y |I = m]

= p1 θ1 + . . . + pm θm .

Note that to estimate θ we only need to estimate the θi ’s since by condition 1 above, the pi ’s are easily
computed. Furthermore, we know how to estimate the θi ’s by condition 2. If we use ni samples to estimate θi ,
then an estimate of θ is given by
θbst,n = p1 θb1,n1 + . . . + pm θbm,nm . (10)
It is clear that θbst,n will be unbiased if for each i, θbi,ni is an unbiased estimate of θi .

Obtaining a Variance Reduction

How does the stratification estimator compare with the usual raw simulation estimator? As was the case with
the game show example, to answer this question we would like to compare Var(θbn ) with Var(θbst,n ). First we
need to choose n1 , . . . , nm such that n1 + . . . + nm = n. That is, we need to determine the number of samples,
ni , that will be used to estimate each θi , but in such a way that the total number of samples is equal to n.
Clearly, the optimal approach would be to choose the ni ’s so as to minimize Var(θbst,n ). Consider for now,
however, the sub-optimal allocation where we set nj := npj for j = 1, . . . , m. Then

Var(θbst,n ) = Var(p1 θb1,n1 + . . . + pm θbm,nm )

σ12 σ2
= p21 + . . . + p2m m
n1 nm
Pm 2
j=1 pj σj
= .
n
On the other hand,
Pm the usual simulation estimator has variance σ 2 /n where σ 2 := Var(Y ). Therefore, we need
only show that j=1 pj σj < σ 2 to prove that the non-optimized stratification estimator has a lower6 variance
2
Pm
than the usual raw estimator. The proof that j=1 pj σj2 < σ 2 is precisely the same as that used for the game
Pm
show example. In particular, equation (8) implies σ 2 = Var(Y ) ≥ E[Var(Y |I)] = j=1 pj σj2 and the proof is
complete!

Optimizing the Stratified Estimator

We know from (10) that
Pn1 (1) Pnm (m)
i=1 Yi i=1 Yi
θbst,n = p1 + . . . + pm
n1 nm
(j)
where for a fixed j, the Yi ’s are IID ∼ Y (j) . This then implies
m
σ12 σ2 X p2j σj2
Var(θbst,n ) = p21 + . . . + p2m m = . (11)
n1 nm j=1
nj

Therefore, to minimize Var(θbst,n ) we must solve the following constrained optimization problem:
m
X p2j σj2
min subject to n1 + . . . + nm = n. (12)
nj
j=1
nj
6 The optimized stratification estimator would then of course achieve an even greater variance reduction.
Monte-Carlo Methods for Risk Management 12

We can easily solve (12) using a Lagrange multiplier. The optimal solution is given by
!
∗ p j σj
nj = Pm n (13)
j=1 pj σj

P 2
m
and the minimized variance is given by Var(θbst,n∗ ) = p
j=1 j jσ /n. Note that the solution in (13) makes
intuitive sense: if pj is large, then other things being equal, it makes sense to expend more effort simulating
from stratum j, i.e., the region where Wj ∈ ∆j . Similarly, if σj2 is large then, other things again being equal, it
makes sense to simulate more often from stratum j so as to get a more accurate estimate of θj .

Remark 1 It is interesting to note the following connection between stratified sampling and importance
sampling. We know that when we importance sample, we like to sample more often from the important region.
The choice of nj in (13) also means that we simulate more often from the important region (in this case the
regions with large σj ’s) when we use optimized stratified sampling.

Remark 2 Note also the connection between stratified sampling and conditional Monte-Carlo. Both methods
rely on the conditional variance formula to prove that they lead to a variance reduction. The difference between
the two methods can best be explained as follows. Suppose we wish to estimate θ := E[Y ] using simulation and
we do this by first generating random variable, W , and then generating Y given W . In the conditional
expectation method, we simulate W first, but then compute E[Y |W ] analytically. In the stratified sampling
method, we effectively generate W analytically, and then simulate Y given W .

Advantages and Disadvantages of Stratified Sampling

The obvious advantage of stratified sampling is that it leads to a variance reduction which can be very
substantial if the stratification variable, W , accounts for a large fraction of the variance of Y . The main
disadvantage of stratified sampling is that typically we do not know the σj2 ’s so it is impossible to compute the
optimal nj ’s exactly. Of course we can overcome this problem by first doing m pilot simulations to estimate
each σj . If we let Np denote the total number of pilot simulations, then a good heuristic is to use Np /m runs
for each individual pilot simulation. In order to obtain a reasonably good estimate of σj2 , a useful rule-of-thumb
is that Np /m should be greater than 30. If m is large however, and each simulation run is computationally
expensive, then it may be the case that a lot of effort is expended in trying to estimate the optimal nj ’s.
One method of overcoming this problem is to abandon the pilot simulations and simply use the sub-optimal
allocation where nj = npj . We saw earlier that this allocation still results in a variance reduction which
sometimes can be substantial. In practice, both methods are used. The decision to conduct pilot simulations
should depend on the problem at hand. For example, if you have reason to believe that the σj ’s will not vary
too much then it should be the case that the optimal allocation and the sub-optimal allocation will be very
similar. In this case, it is probably not worth doing the pilot simulations. On the other hand, if the σj ’s vary
considerably, then conducting the pilot runs may be worthwhile. Of course, a combination of the two is also
possible where a only a subset of the pilot simulations is conducted.

The stratified simulation algorithm is given below. We assume that the pilot simulations have already been
completed, or it has been decided not to conduct them at all; either way, the nj ’s have been pre-computed. We
2
also show how the estimate, θbn,st , and the estimated variance, σbn,st , can be computed without having to store
P (j) 2 P (j)
all the generated samples. That is, we simply keep track of Yi and Yi for each j since these
7 2
quantities are all that is required to compute θn,st and σ
b bn,st .
 Pnj (j)
  Pnj (j)

i=1 Yi i=1 Yi
that θbn,st = m
7 Note
Pm
p2j . Both quantities can be estimated
P
j=1 nj
pj and Var(θbn,st ) = j=1 Var nj
P nj (j) P nj (j) 2
knowing just i=1 Yi and i=1 Yi . Any simulation study that requires a large number of samples should only keep track
of these quantities, thereby avoiding the need to store every sample.
Monte-Carlo Methods for Risk Management 13

Stratification Simulation Algorithm for Estimating θ

2
set θbn,st = 0; σ
bn,st = 0;
for j = 1 to m
set sumj = 0; sum squaresj = 0;
for i = 1 to nj
(j)
generate Yi
(j)
set sumj = sumj + Yi
(j) 2
set sum squaresj = sum squaresj + Yi
end for
set θj = sumj /nj

set σbj2 = sum squaresj − sum2j /nj /(nj − 1)
set θbn,st = θbn,st + pj θj
2 2
set σbn,st =σ bn,st +σbj2 p2j /nj
end for
set approx. 100(1 − α) % CI = θbn,st ± z1−α/2 σ
bn,st

3.2 Some Applications of Stratified Sampling to Option Pricing

Our examples relate to option pricing in a geometric Brownian motion (GBM) setting. Since GBM is a very
poor model of asset price dynamics, however, these examples are not directly applicable in practice. However,
the specific techniques and results that we discuss are certainly applicable. For example, we know a multivariate
t random vector can be generated using a multivariate normal random vector together with an independent
chi-squared random variable. Therefore the techniques we discuss below for multivariate normal random vectors
can also be used8 for multivariate t random vectors. They are also applicable more generally to normal-mixture
distributions, Gaussian and t copulas, as well as the simulation of SDE’s. As stated earlier, Chapter 9 of
Glasserman (2004) should be sufficient to convince any reader of the general applicability of stratified sampling
(as well as importance sampling) to risk management.

Example 10 (Pricing a European Call Option)

Suppose we wish to price a European call option where St ∼ GBM (r, σ 2 ). Then
C0 = E e−rT max(0, ST − K) = E[Y ]
 

 2
√ 
where Y = h(X) := e−rT max 0, S0 e(r−σ /2)T +σ T X − K for X ∼ N(0, 1). While we know how to
compute C0 analytically, it is worthwhile seeing how we could estimate it using stratified simulation. Let
W = X be our stratification variable. To see that we can stratify using this choice of W note that:

(1) Computing P(W ∈ ∆)

For ∆ ⊆ R, P(W ∈ ∆) can easily be computed. Indeed, if ∆ = [a, b], then P(W ∈ ∆) = Φ(b) − Φ(a), where
Φ(.) is the CDF of a standard normal random variable.

(2) Generating (Y |W ∈ ∆)
(h(X)|X ∈ ∆) can easily be generated. We do this by first generating X̃ := (X|X ∈ ∆) and then take h(X̃).
We generate X̃ as follows. First note that if X ∼ N(0, 1), then we can generate an X using the inverse
8 See Chapter 9 of Glasserman (2004) for a risk management application in a multivariate t setting.
Monte-Carlo Methods for Risk Management 14

transform method by setting X = Φ−1 (U ). The problem with such an X is that it may not lie in ∆ = [a, b].
However, we can overcome this problem by simply generating Ũ ∼ U (Φ(a), Φ(b)) and then setting
X̃ = Φ−1 (Ũ ). It is then straightforward to check that X̃ ∼ (X|X ∈ [a, b]).
It is therefore clear that we can estimate C0 using X as a stratification variable.

Example 11 (Pricing an Asian Call Option)

Recall that the discounted payoff of an Asian call option is given by
 Pm 
−rT i=1 SiT /m
Y := e max 0, −K (14)
m

and that it’s price is given by Ca = E[Y ] where as before we assume St ∼ GBM (r, σ 2 ). Now each SiT /m may
be expressed as r !
2 iT T
SiT /m = S0 exp (r − σ /2) +σ (X1 + . . . + Xi ) (15)
m m

where the Xi ’s are IID N(0, 1). This means that we may then write Ca = E [h(X1 , . . . , Xm )] where the function
h(.) is given implicitly by (14) and (15). To estimate Ca using our standard simulation algorithm, we would
simply generate sample values of h(X1 , . . . , Xm ) and take their average as our estimate. We can also, however,
estimate Ca using stratified9 sampling.
To do so, we must first choose a stratification variable, W . One possible choice would be to set W = Xj for
some j. However, thisP is unlikely to capture much of the variability of h(X1 , . . . , Xm ). A much better choice
m
would be to set W = j=1 Xj . Of course, we need to show that such a choice is possible. That is, we need to
show that P(W ∈ ∆) is easily computed, and that (Y |W ∈ ∆) is easily generated.

(1) Computing P(W ∈ ∆)

Since X1 , . . . , Xm are IID N(0, 1), we immediately have that W ∼ N(0, m). If ∆ = [a, b] then

P(W ∈ ∆) = P (N(0, m) ∈ ∆) P (a ≤ N(0, m) ≤ b)

=
 
a b
= P √ ≤ N(0, 1) ≤ √
m m
   
b a
= Φ √ −Φ √ .
m m
   
Similarly, if ∆ = [b, ∞), then P(W ∈ ∆) = 1 − Φ √bm , and if ∆ = (−∞, a], then P(W ∈ ∆) = Φ √am .

(2) Generating (Y |W ∈ ∆)
We need two results from the theory of multivariate normal random variables. The first result is well known to
us:

1. Suppose X = (X1 , . . . , Xm ) ∼ MVN(0, Σ). If we wish to generate a sample vector X, we first generate
Z ∼ MVN(0, Im ) and then set
X = CT Z (16)
where CT C = Σ. One possibility of course is to let C be the Cholesky decomposition of Σ, but in fact
any matrix C that satisfies CT C = Σ will do.

9 The method we now describe is also useful for pricing other path dependent options. See Glasserman (2004) for further

details.
Monte-Carlo Methods for Risk Management 15
p
2. Let a = (a1 a2 . . . am ) satisfy ||a|| = 1, i.e. a21 + . . . + a2m = 1, and let
Z = (Z1 , . . . , Zm ) ∼ MVN(0, Im ). Then
( m
)
 X
(Z1 , . . . , Zm )  ai Zi = w ∼ MVN(wa> , Im − a> a).

i=1

Pm
Therefore, to generate {(Z1 , . . . , Zm ) | i=1 ai Zi = w} we just need to generate a vector, V, where
V ∼ MVN(wa> , Im − a> a) = wa> + MVN(0, Im − a> a).
>
Generating such a V is very easy since (Im − a> a) (Im − a> a) = Im − a> a. That is, Σ> Σ = Σ where
Σ = Im − a> a, so we can take C = Σ in (16).

We can now return to the problem of generating (Y | W ∈ P ∆). Since Y = h(X1 , . . . , Xm ), we can clearly
m
generate (Y | W ∈ ∆) if we can generate [(X1 , . . . , Xm ) | i=1 Xi ∈ ∆]. To do this, suppose again that
∆ = [a, b]. Then
" m
# " m  #
 X  1 X a b
(X1 , . . . , Xm )  Xi ∈ [a, b] ≡ (X1 , . . . , Xm )  √ Xi ∈ √ , √ .
 
i=1
m i=1 m m
Pm
Now we can generate [(X1 , . . . , Xm ) | i=1 Xi ∈ ∆] in two steps:

h Pm  Pm h ii
Step 1: Generate √1 X √1 X ∈ √a , √b . This is easy to do since

m i=1 i i=1
m i
 m m

Pm  h i
1 √a , √b
√ X ∼ N(0, 1) so we just need to generate N(0, 1) N(0, 1) ∈ which we can do using

m i=1 i  m m
the method described in Example 10. Let w be the generated value.

h  Pm i
Step 2: Now generate (X1 , . . . , Xm )  √1 Xi = w which we can do by the second result above and

m i=1
the comments that follow it.

Example 12 (Pricing a Barrier Option)

Recall again the problem of pricing an option that has payoff

max(0, ST − K1 ) if ST /2 ≤ L,
h(X) =
max(0, ST − K2 ) otherwise.
where X = (ST /2 , ST ). We can write the price of the option as
h  i
C0 = E e−rT max(0, ST − K1 )I{ST /2 ≤L} + max(0, ST − K2 )I{ST /2 >L}

where we again assume that St ∼ GBM (r, σ 2 ). Using conditional Monte-Carlo, we can write (why?)
C0 = E[Y ] where
 
Y := e−rT /2 c(ST /2 , T /2, K1 , r, σ)I{ST /2 ≤L} + c(ST /2 , T /2, K2 , r, σ)I{ST /2 ≥L} (17)

and where c(x, t, k, r, σ) is the price of a European call option with strike k, interest rate r, volatility σ, time to
maturity t, and initial stock price x.

Question 1: Having conditioned on ST /2 , could we now also use stratified sampling?

Question 2: Could we use importance sampling?
Question 3: What about using importance sampling before doing the conditioning?
Monte-Carlo Methods for Risk Management 16

4 A Brief Introduction to Counterparty Risk and CVA

We now briefly discuss10 the concept of credit valuation adjustment (CVA) and how Monte-Carlo is typically
required to calculate it. Since the 2008 financial crisis CVA has played an increasingly important role in the
valuation of OTC derivatives. Whereas the pre-crisis world generally ignored counterparty credit risk in the
pricing of derivatives, in the post-crisis world OTC derivative security prices must now reflect the possibility of
counterparty default.

Measuring Counterparty Credit Exposure: PFE and EPE

Counterparty credit exposure is the exposure to loss due to the failure of a counterparty to fulfill their
obligations. Note that counterparty risk is at the root of traditional banking since the earliest form of financial
instruments were bonds and their value was driven by the perceived credit worthiness of the issuers. The key
statistical measures of counterparty credit exposure are potential future exposure(PFE) and expected
positive exposure (EPE). Given a security (or portfolio of securities with the same counterparty) with
maturity T , PFEt for any t ∈ [0, T ] is defined to be the quantile (usually 97.5% or 99%) of the time t price
distribution as determined using the market information today at time 0. Similarly, EPEt is defined to be the
mean of the positive part of the time t price distribution. In Figure 1 we can see the price distributions at
various points in time for a plain vanilla swap. We can also see PFEt and EPEt plotted as a function of t.

Figure 1: Computing exposures of a vanilla swap by simulation. (Courtesy of Giovanni Cesari and UBS)

Credit Value Adjustment (CVA)

The CVA of a derivative security is defined to be the price of the counterparty credit exposure for that security.
It is an adjustment to the price of the derivative that takes into account counterparty credit exposure. In order
to calculate the CVA consider the loss (to us) that occurs on the derivative position should the counterparty
default at time t < T . In that case our loss is given by

Losst = (1 − RV ) 1{τ =t} Vτ+ (18)

where τ denotes the random default time of the counterparty, Vt denotes the time t value (from our
perspective) of the derivative security and RV is the recovery rate. The CVA is the time t = 0 value of the sum
10 We will closely follow a tutorial presentation of Giovanni Cesari.
Monte-Carlo Methods for Risk Management 17

of the losses, Losst , where the sum is taken over all t ∈ [0, T ]. Like all security prices, it is determined via
risk-neutral pricing and so we obtain
Z T  + 
Q Vu
CVA0,T = (1 − RV ) E0 1{τ =u} du (19)
0 Nu
where Q denotes the risk-neutral probability measure corresponding to the numeraire price process Nt .
Evaluating (19) is generally computationally intensive. If Vu+ /Nu is Q-independent of τ then some simplification
occurs as we can then model the default event separately from the price process of the derivative security.
Moreover, CDS prices in the market-place can typically then be used to work with the term EQ

0 1{τ =u} which
then appears in the integrand in (19). When Vu+ /Nu and τ are Q-dependent, however, then matters are
considerably more complex and one must now account for the possibility of wrong-way risk in the modeling.
Very often analytic expressions are not available for the price Vt in which case it must be computed numerically
or via Monte-Carlo simulation. If the derivative security is Bermudan / American and therefore allows for the
possibility of early exercise then matters become even more complicated as an optimal stopping problem must
then be (approximately) solved11 at each time point t along each simulated path. In general then, evaluating
CVA0,T is demanding and often requires nested Monte-Carlo simulations.
In determining (18) we have only taken our point of view. However our counterparty will also perceive the
possibility that we will default on them and therefore they too will compute a CVA. The sum of the two CVAs is
known as bilateral CVA and it is this CVA that should be accounted for in the price of the derivative security.
Note that the value of our counterparty’s CVA is our debit value adjustment (DVA).

Exercise 2 Suppose the underlying derivative security moves in our favor. Do you think we will record a gain
or loss on our CVA? What about our DVA?

Differences Between CVA and Counterparty Exposure

Both CVA and counterparty exposure measures such as PFE and EPE both require price distributions at
different times in the future. Counterparty exposures are statistical measures, however. PFE is sometimes used
against credit limits while EPE is used for computing risk-weighted assets and corresponding capital
requirements. In contrast CVA is the cost of buying protection on the counterparty with the payout being the
portfolio value (less recovery) in the case of default. CVA is often hedged and enters into the calculation of VaR
and other risk measures.

5 Estimating the Greeks Via Monte Carlo

We now briefly consider three different methods for estimating the Greeks, i.s. parameter sensitivities, for
derivative securities. The first approach is based on finite difference approximations and therefore results in
biased estimators. The other two approaches, namely the pathwise and likelihood ratio methods result (when
applicable) in unbiased estimators. All three approaches are often used in practice and Chapter 7 of
Glasserman’s “Monte Carlo Methods in Financial Engineering” should be consulted for further details. We also
follow Glasserman in our discussion here.
To begin, let C(θ) := E0 [Y (θ)] be the price of a particular derivative security when the value of some parameter
is θ. Then
∂C
α(θ) :=
∂θ
is the derivative price’s sensitivity to changes in the parameter θ. For example, if C is the value of a standard
European call option in the Black-Scholes framework and θ = S0 , the underlying security price, then α(θ) is the
delta of the option and it can be calculated explicitly. In general, however, we will not have an explicit
expression for α(θ) and in such circumstances we might instead use Monte-Carlo methods to estimate it.
11 The cross-path regression approach is then typically employed as it can provide prices of the Bermudan / American security

simultaneously for all time periods and sample paths.

Monte-Carlo Methods for Risk Management 18

5.1 Finite Difference Approximations and Common Random Numbers

One approach to estimating α(θ) is to simply simulate the finite-difference ratio

Y (θ + ) − Y (θ − )
∆ := (20)
2
for some small given  > 0. We could therefore estimate α(θ) (with some hopefully small bias) by simulating
samples of ∆ and then averaging them. Noting that
Var (Y (θ + )) + Var (Y (θ − )) − 2 Cov (Y (θ + ), Y (θ − ))
Var (∆ ) = (21)
42
it is clear that a reduction in variance will follow if we can arrange our Monte-Carlo samples so that
Cov (Y (θ + ), Y (θ − )) > 0. Therefore Y (θ + ) and Y (θ − ) should not be estimated independently but
instead they should be estimated using common random numbers, i.e. generating them from the same12 U (0, 1)
random variable.
For small , the variance reduction from using common random numbers can be dramatic. That said, there is
clearly a tradeoff between bias and variance in our selection of . In general it can be shown that we should
choose  = O(n−1/5 ) in (20) where as usual n is the number of samples. In that case convergence of the
estimator in (20) is O(n−2/5 ) although O(n−1/2 ) convergence can be obtained if Y is continuous in θ (in which
case  should be taken as small as possible).
Similar results are available for estimating second-order derivatives such as gamma, vanna and volga but the
convergence rates are not as good and estimating these quantities is fundamentally harder than estimating first
order derivatives.

5.2 The Pathwise Method

Recalling that C(θ) := E0 [Y (θ)], the pathwise estimator is calculated by interchanging the order of
differentiation and integration to obtain
 
∂ ∂Y (θ)
α(θ) = E0 [Y (θ)] = E0 . (22)
∂θ ∂θ

Assuming the interchange of differentiation and integration in (22) can be justified, we could therefore use
∂Y (θ)
∂θ as an unbiased estimator of α(θ). In order to operationalize (22), however, one must first explicitly state
the relationship between Y and the parameter θ.

Example 13 (Estimating the Black-Scholes Delta)

Consider the case of a European call option with strike K and maturity T in the Black-Scholes framework.
While an explicit expression for the option delta is available, we can also estimate it via the pathwise method as
follows. We first write the option payoff as
+
Y = e−rT (ST − K) (23)
 2
 √
r− σ2 T +σ T Z
ST = S0 e (24)

where Z ∼ N(0, 1). It follows from (23) and (24) that13

∂Y ∂Y ∂ST
=
∂S0 ∂ST ∂S0
ST
= e−rT 1{ST >K} . (25)
S0
12 In practice we would simulate n samples of ∆ and use a single sequence U , . . . , U of IID U (0, 1) random variables to
 1 n
generate these samples.
13 While Y is not differentiable with respect to S when S = K this event happens with probability 0 and so the derivative
T T
exists almost surely.
Monte-Carlo Methods for Risk Management 19

The estimator (25) is easily calculated via a Monte-Carlo simulation. It should also be clear that the expression
(25) is valid for any model of security prices where St = S0 eXt for any (risk-neutral) stochastic process Xt that
does not depend on S0 .
While we haven’t justified interchanging the order of differentiation and integration in Example 13, a general
rule of thumb is that this interchange can be justified when the payoff Y is continuous in θ and that is clearly
the case there. In contrast, the interchange is generally invalid when Y is not continuous in θ. This means in
particular that the pathwise method does not work in general for barrier and digital options.
Exercise 3 Find an analogous expression to (25) for estimating the vega of a call option in the Black-Scholes
world.

5.3 The Likelihood Ratio Method

In contrast to the pathwise method, the likelihood ratio method differentiates a probability density (rather than
a discounted payoff Y ) with respect to the parameter of interest, θ. It provides a good potential alternative to
the pathwise method when Y is not continuous in θ. In order to develop the method we now write the payoff
Y = f (X1 , . . . , Xm ) as a function of the random vector X = (X1 , . . . , Xm ) whose components could for
example represent the price of an underlying security at different dates, or the prices of several underlying
securities at the same date. We assume that X has a density g and that θ is a parameter of this density. We
will therefore write gθ for this density and use Eθ to denote that expectations are taken with respect to gθ . We
can therefore write Z
Eθ [Y ] = f (x)gθ (x) dx. (26)
Rm
We can now differentiate across (26) to obtain
∂
α(θ) = Eθ [Y ]
∂θ
Z
∂
= f (x) gθ (x) dx (27)
R m ∂θ
where we have assumed the interchange of differentiation and integration is again justified. Writing ġθ for
∂gθ /∂θ we can multiply and divide the integrand in (27) by gθ to obtain
Z
ġθ (x)
α(θ) = f (x) gθ (x) dx
Rm gθ (x)
 
ġθ (X)
= Eθ f (X) . (28)
gθ (X)
The ratio ġθ (X)/gθ (X) is known as the score function. While the interchange of differentiation and
integration needs to be justified this is typically not a problem since (unlike option payoffs) density functions are
usually smooth functions of their parameters. It’s also worth noting that there is considerable flexibility in
whether we choose to view θ as a parameter of the payoff Y or of the density g. In (24), for example, it is clear
that S0 is a parameter of the path and not of the density (which is N(0, 1)). But we could also have written
(how?) the density as a function of S0 . In that case we can (after some algebra) obtain
 
∂C Z
= E e−rT (ST − K)+ √ (29)
∂S0 S0 σ T
where Z is a standard normal random variable. The expression inside the expectation in (29) is therefore our
likelihood ratio estimator for the option delta.
In practice, use of the likelihood ratio method tends to be limited by either not having gθ available explicitly or
because the estimator’s variance may be too high. We can see from (29), for example, that the variance of the
likelihood ratio estimator will be very high when T is close to 0. This can be a serious14 problem for the method
more generally when the payoff of the derivative security depends on the underlying price at a range of times
with small increments between them.
14 This is related to importance sampling and associated absolute continuity conditions. See Glasserman for further details.