MARKOVIAN PROJECTION METHOD FOR VOLATILITY CALIBRATION

VLADIMIR PITERBARG

Abstract: We present the “Markovian projection method”, a method to obtain closed-

form approximations to European option prices on various underlyings that, in principle,
is applicable to any (diffusive) model. Successful applications of the method have already
appeared in the literature, in particular for interest rate models (short rate and forward
Libor models with stochastic volatility), and interest rate/FX hybrid models with FX
skew. The purpose of this note is thus not to present other instances where the Markovian
projection method is applicable (even though more examples are indeed given) but to
distill the essence of the method into a conceptually simple “plan of attack”, a plan that
anyone who wants to obtain European option approximations can follow.

First version: 25 May 2006

This version: 25 May 2006

1. I NTRODUCTION
European-style options are usually the most liquid options available in any market.
More often than not, they are the only options that are liquid enough to be used for
model calibration. Thus, efficient methods for valuing European-style options are a criti-
cal requirement for any model. In this paper we develop the Markovian projection (MP)
method, a very general and powerful approach to deriving accurate approximations to
European option prices in a wide range of models.
Various ideas related to the Markovian projection method have already been used, and
very successfully so, in a wide range of contexts. In Piterbarg (2004) and Piterbarg
(2005b), closed-form approximations to European swaptions prices in a forward Libor
model with stochastic volatility have been obtained. In Andreasen (2005), that same has
been accomplished for a short rate model with stochastic volatility. Piterbarg (2006a)
and Piterbarg (2005a) have introduced an FX skew in the interest rate/FX hybrid model
and demonstrated how closed-form approximations for FX options can be derived. Same
techniques can easily be adapted to interest rate/equity hybrids, for example. In this note,
we formalize the ideas from these papers as the MP method.
The Markovian projection method is based on combining pioneering ideas of Dupire
(see Dupire (1997)) on local and stochastic volatility (and recently rediscovered, in the
quantitative finance context, results of Gyöngy, see Gyöngy (1986)) with parameter-
averaging techniques that have been developed for stochastic volatility models relatively
Key words and phrases. Local volatility, stochastic volatility, Markovian projection, parameter averag-
ing, Dupire’s local volatility, index options, basket options, spread options.
I would like to thanks my colleagues at Barclays Capital, as well as Leif Andersen and Martin Forde, for
valuable comments and discussions.
1
2 VLADIMIR PITERBARG

recently. Both Dupire and Gyöngy show that if two underlyings, each governed by its
own SDEs, are given, and the expected values of the SDE coefficients conditioned on the
appropriate underlying (the so-called “Dupire’s local volatilities”) are equal, then Euro-
pean options prices on (or, equivalently, one-dimensional marginal distributions of) the
two underlyings are the same. This observation is the key tool to replace a model with
complicated dynamics for the underlying with a much simpler one. Once that is accom-
plished, the (local or stochastic volatility) model that is obtained is further simplified, by
applying the parameter averaging techniques, to a form that admits closed form solutions.
The MP method proceeds in four steps.
Step 1: For the underlying of interest, its SDE, driven by a single Brownian motion,
is written down by computing its quadratic variance (and combining all dt terms
if exist).
Step 2: In this SDE, the diffusion and drift (if exists) coefficients are replaced with
their expected values conditional on the underlying. As shown in Section 2, this
does not affect the values of European options.
Step 3: The conditional expected values from Step 2 are computed or, more com-
monly, approximated. Methods based on, or related to, Gaussian approximations
are often used for this step, see Section 4 below.
Step 4: Parameter averaging techniques (see Piterbarg (2005c), Piterbarg (2005b))
are used to relate the time-dependent coefficients of the SDE obtained in Step 3 to
time-independent ones. This, typically, allows for a quick and direct computation
of European option values.

2. O NE - DIMENSIONAL DISTRIBUTIONS OF I TO PROCESSES

Suppose a model is given. Let S (t) be the quantity of interest, European options on
which are to be valued. For example, it could be a forward price of an equity share, a
forward FX rate, a Libor or a swap rate, and so on. The market observable may or may
not be a “primitive” variable in the model. For example, equity models are most often
expressed in terms of the share price, and then the dynamics are imposed on S (·) directly.
On the other hand, in the interest rate modeling, the dynamics are often imposed on short
rates or on Libor rates. If S (·) is the swap rate, then it is not a “primitive” variable in such
a model; the dynamics of S (·) are implicitly defined by the model expressed in terms of
other variables. In either case, an Ito process for S (·) , given by the model, can always be
written down. It is of the general form
(1) dS (t) = Σ (t, . . . ) dW (t) .
Here Σ (t, . . . ) is, in general, a stochastic process that may depend on S (·) and other
sources of randomness. Note that we do not include the dt term. It will be clear from
what follows that this is not a limitation, as the drifts can be handled as well. We write
the process without the dt term largely by convention, as most often there is measure un-
der which S (·) is a martingale, and it is precisely the measure that we are interested in
studying the dynamics of S (·) under (eg the swap measure for a swap rate, or the forward
measure for a forward FX rate).
 MARKOVIAN PROJECTION 3

Another feature of the way (1) is written is the fact that the Brownian motion dW (·)
is one-dimensional. Again, this is not a limitation of the method. Even if S (·) is driven
by a multi-dimensional Brownian motion, a one-dimensional Brownian motion, and the
diffusion coefficient, can always be found such that the law of the original process and
the new one are exactly the same. Let us elaborate. Suppose S follows
d
(2) dS (t) = ∑ Σi (t, . . . ) dWi (t) ,
i=1
where (Σ1 (t, . . . ) , . . . , Σd (t, . . . )) is a d-dimensional process and (W1 (t) , . . . ,Wd (t)) is
 d
a d-dimensional Brownian motion with the correlation matrix ρi j i, j=1 . Then, define
Σ (t, . . . ) by
d
Σ2 (t, . . . ) = ∑ Σi (t, . . . ) Σ j (t, . . . ) ρi j
i, j=1
and dW (t) by !
d
dW (t) = Σ−1 (t, . . . ) ∑ Σi (t, . . . ) dWi (t) .
i=1
Simple quadratic variation calculations confirm that dW is indeed a (one-dimensional)
Brownian motion and the law of the process S (·) given by (2) is the same as the one given
by (1).
To price a European option on S (·) with expiry T and strike K, we need to know
the one-dimensional distribution of S (·) (at time T ). On the other hand, the knowledge
of prices of all European options for all strikes K for a given expiry T is equivalent to
knowing the (one-dimensional) distribution of S (·) at time T. In view of these remarks,
the importance of the following result is evident ((Gyöngy, 1986, Theorem 4.6)).
Theorem 2.1 (Gyongy 1986). Let X (t) be given by
(3) dX (t) = α (t) dt + β (t) dW (t) ,
where α (·) , β (·) are adapted bounded stochastic processes such that (3) admits a unique
solution. Define a (t, x) , b (t, x) by
a (t, x) = E (α (t)| X (t) = x) ,
b2 (t, x) = E β 2 (t) X (t) = x .

Then the SDE

(4) dY (t) = a (t,Y (t)) dt + b (t,Y (t)) dW (t) ,
Y (0) = X (0) ,
admits a weak solution Y (t) that has the same one-dimensional distributions as X (t) .
Remark 2.1. The process Y (·) follows what is known in financial mathematics as a “local
volatility” process. The function b (t, x) is often called
“Dupire’s local volatility” for
X (·), after Dupire who directly linked E β 2 (t) X (t) = x to prices of European options
on X (·) when α (·) ≡ 0 (see Dupire (1997)), a link that can be seen to lead to the same
result.
4 VLADIMIR PITERBARG

Remark 2.2. Since X (·) and Y (·) have the same one-dimensional distributions, the prices
of European options on X (·) and Y (·) for all strikes K and expiries T will be the same.
Thus, for the purposes of European option valuation and/or calibration to European op-
tions, a potentially very complicated process X (·) can be replaced with a much simpler
Markov process Y (·) , the Markovian projection1 of X (·) .
Remark 2.3. Efficient valuation and approximation methods for European options in lo-
cal volatility models, ie models of the type (4), are available, see for example
Andersen and Brotherton-Ratcliffe (2001).
Proof. The original proof in Gyöngy (1986) is fairly involved and we do not reproduce it
here. Instead, we present an outline of the proof inspired by Savine (2000), a much more
financially-motivated approach. As explained above, the case α (t) = 0 is most common
in financial applications, so that is what we consider.
Denote
c (t, K) = E0 (X (t) − K)+ .

Here {c (t, K)}t,K are the values of European call options on X (·) for expiries t and strikes
K. Dupire in Dupire (1994) was the first to find the expression for the local volatility in
a local volatility model that reproduces all European options {c (t, K)}t,K for all expiries
t and strikes K. The key insight of this proof is to note that the same technique can be
applied to European option prices obtained from another model, instead of the market
option prices as originally done by Dupire.
As explained above, matching all European option prices (for all expiries and strikes) is
equivalent to matching all one-dimensional distributions. “Dupire local volatility” b (t, x)
in the model
dY (t) = b (t,Y (t)) dW (t)
for the set of option prices {c (t, K)} is given by
∂
2 ∂t c (t, K)
(5) b (t, K) = 2 ∂ 2 ,
∂K 2 c (t, K)
see Dupire (1994). To compute the right-had side, we first write (the use of delta-functions
in the integrands can be justified by Tanaka’s formula, see Karatzas and Shreve (1997))
1
d (X (t) − K)+ = 1{X(t)>K} dX (t) + δ{X(t)=K} d hX (t)i ,
2
and, since X (t) is a martingale under the measure considered,
1 t
Z
E (X (t) − K)+ − (X (0) − K)+ =


E δ{X(s)=K} d hX (s)i .
2 0
Clearly


E δ{X(t)=K} d hX (t)i = E δ{X(t)=K} × E (d hX (t)i| X (t) = K)
and

∂2 + ∂2
E δ{X(t)=K} = E (X (t) − K) = c (t, K) .
∂ K2 ∂ K2
1Often we will use the same letter to denote both, the original process and its Markovian projection.
 MARKOVIAN PROJECTION 5

From (3)
d hX (t)i = β 2 (t) dt,
so that
∂2
c (t, K)2 × E β 2 (t) X (t) = K dt.



E δ{X(t)=K} d hX (t)i = 2
∂K
In particular,
∂ ∂
E (X (t) − K)+ − (X (0) − K)+


c (t, K) =
∂t ∂t
1 ∂2 2


= × c (t, K) × E β (t) X (t) = K .
2 ∂ K2
Substituting this equality into (5) we obtain
b2 (t, K) = E β 2 (t) X (t) = K ,

and the theorem is proved.

For the variable S (·) considered above, the theorem immediately implies that we can
replace (1) with the Markovian projection
(6) dS (t) = Σ̃ (t, S (t)) dW (t)
to price European options, where Σ̃ (t, x) , the deterministic function of time and state, is
obtained by
Σ̃2 (t, x) = E Σ2 (t, . . . ) S (t) = x .


(7)
We emphasize that the Markovian projection is exact for European options but, of course,
does not preserve the dependence structure of S (·) at different times or, equivalently,
the marginal distributions of the process of orders higher than one. Thus, the prices of
securities dependent on sampling of S (·) at multiple times, such as barriers, American
options and so on, are different between the original model (1) and the projected model
(6).
In practice, the conditional expected value in (7) needs to be computed/approximated.
This is often impossible to do exactly, and various approximations need to be employed.
We discuss common approaches later in the paper. Next, however, we consider another
ways that Theorem 2.1 can be used.

3. A PPLICATIONS TO STOCHASTIC VOLATILITY MODELS

Suppose a stochastic volatility model is given. Let S (·) be the variable of interest as
before, governed in the model by the equation (1). The results of Theorem 2.1 allow us to
replace the dynamics with local volatility dynamics for the purposes of pricing European
options. However, substituting a local volatility model for a stochastic volatility model is
not always ideal. The two types of models are very different, and that makes computing
the conditional expected value in (7) particularly difficult. This renders various possible
approximations rather inaccurate. As a general rule, when one model is approximated
with another, it is always best to keep as many features as possible the same in the two
6 VLADIMIR PITERBARG

models. For example, a stochastic volatility model should be approximated with a (poten-
tially simpler) stochastic volatility model, a model with jumps by a (simpler) model with
jumps, and so on.
Upon some reflection, it should be clear that Theorem 2.1 provides us with a tool that
is more general than just a way of approximating any model with a local volatility one.
The following trivial Corollary to Theorem 2.1 provides the required framework.
Corollary 3.1. If two processes have the same Dupire’s local volatility, the European
option prices on both are the same for all strikes and expiries.
Consider a stochastic volatility model. Let X1 (t) follow
p
dX1 (t) = b1 (t, X1 (t)) z1 (t) dW (t) ,
where z1 (t) is some variance process. Suppose we would like to match the European
option prices on X1 (·) (for all expiries and strikes) in a model of the form
p
dX2 (t) = b2 (t, X2 (t)) z2 (t) dW (t) ,
where z2 (t) is a different, and potentially simpler, variance process. Then Corollary 3.1
and Theorem 2.1 imply that to achieve that, we need to choose b2 (t, x) such that
E b22 (t, x) z2 (t) X2 (t) = x = E b21 (t, x) z1 (t) X1 (t) = x ,

or, rearranging the terms,

E (z1 (t)| X1 (t) = x)
(8) b22 (t, x) = b21 (t, x) .
E (z2 (t)| X2 (t) = x)
To apply this formula, the conditional expected values need to be computed. Note that
it is beneficial to choose z2 as similar to z1 as possible. Then, various approximations
employed to compute both conditional expected values will tend to cancel each other,
thus increasing the accuracy of calculating b2 (t, x).
It is generally the case that a realistic stochastic volatility model with a local volatility
component should have a local volatility that does not have a lot of curvature (such as
a linear or a power function of the underlying). Moreover, such a choice makes many
calculations and approximations much easier, as then local volatility functions can just be
described by two numbers, a value and a slope (at the forward value of the underlying,
for example). In practice it means that when a stochastic volatility model is constructed
to fit implied volatility smiles (be that market volatilities or implied volatilities from a
different model), it is better to choose the stochastic volatility that explains as much of the
curvature of the implied volatility smile as possible. Then the local volatility function will
be “almost linear”, and the benefits mentioned above will apply. For the concrete example
of equation (8), one should typically try to find z2 (·) that has the same term variance as
z1 (·) , as it is the volatility of stochastic variance that defines the curvature of the smile.
Another typical application of Corollary 3.1 proceeds as follows. Suppose
p
dX1 (t) = β (t) z (t) dW (t) ,
where z (t) is the stochastic variance process, but β (t) is not deterministic, but a stochastic
process in its own right (eg a complicated function of the state variables in a term structure
 MARKOVIAN PROJECTION 7

model of interest rates, with S (·) being a swap rate). We would like to replace the SDE
with a local-stochastic volatility model,
p
dX2 (t) = b (t, X2 (t)) z (t) dW (t) .
Then, according to Corollary 3.1, we need to set
E β 2 (t) z (t) X1 (t) = x


2
(9) b (t, x) = .
E (z (t)| X2 (t) = x)
This formula can be simplified when, for example, β (·) and z (·) are (approximately)
conditionally independent given X1 (·). Then

E (z (t)| X1 (t) = x)
(10) b2 (t, x) ≈ E β 2 (t) X1 (t) = x .
E (z (t)| X2 (t) = x)
In many real-world situations, it can be safely assumed that E (z (t)| X1 (t) = x) ≈ E (z (t)| X2 (t) = x),
in which case the formula simplifies further,
b2 (t, x) ≈ E β 2 (t) X1 (t) = x .


(11)

4. G AUSSIAN APPROXIMATIONS FOR CONDITIONAL EXPECTED VALUE

CALCULATIONS

Recall the set up of Section 2, namely the general equation of the type (1) for an under-
lying of interest. Calculations of conditional expected values, as in (7) (or (8)) are often
the most difficult part of the MP method. Tractable formulas are hard to find for
any type
of process, with one notable exception. If in (7) the process S (t) , Σ (t, . . . ) is jointly
2

Gaussian, then the calculations are straightforward. While this is the case in only the
trivial situation of Σ2 (t, . . . ) being deterministic, the availability of closed-form solutions
in the Gaussian case gives rise to the approximation in which the actual dynamics of S (t)
and Σ2 (t)(= Σ2 (t, . . . )) are replaced by the Gaussian ones. While the details vary from
case to case, the general approach is given by the following proposition.
Proposition 4.1. Let the dynamics of S (t) , Σ2 (t) be written in the following form,

(12) dS (t) = Σ (t) dW (t) ,

dΣ2 (t) = η (t) dt + ε (t) dB (t) ,
where Σ (t) , η (t) and ε (t) are adapted stochastic processes (that may depend on S (·) ,
Σ2 (·) , etc) and W (t) , B (·) are both Brownian motions. Then the conditional expected
value in (7) can be approximated by
Σ̃2 (t, x) , E Σ2 (t) S (t) = x


(13)
= Σ2 (t) + r (t) (x − S0 ) ,
0 ε̄ (s) Σ̄ (s) ρ̄ (s) ds
Rt
r (t) = ,
0 Σ̄ (s) ds
Rt
2
8 VLADIMIR PITERBARG

where ε̄ (t) , Σ2 (t) , ρ̄ (t) , Σ̄ (t) are, respectively, deterministic approximations to ε (t) ,
Σ2 (t), ρ (t) , hdW (t) , dB (t)i /dt, and Σ (t) . In particular, one can take
ε̄ (t) = Eε (t) ,
Z t
Σ2 (t) = (Eη (s)) ds,
0
ρ̄ (t) = E (hdW (t) , dB (t)i /dt) ,
q
Σ̄ (t) = Σ2 (t).

Proof. If (X,Y ) is a Gaussian vector, then

covar (X,Y )
(14) E (X|Y ) = EX + (Y − EY ) .
Var (Y )
The result follows by computing the required moments in a Gaussian approximation to
the dynamics of (12).
Remark 4.1. If more is known about the form of the dynamics (12), other deterministic
approximations to ε (t) (and other parameters) can be more useful.
Remark 4.2. The approximation (13) is linear in x. The Gaussian approximation, hence,
works best when the local volatility functions that are used are either linear or close to
linear. This provides another angle to the comment at the end of Section 3, namely that
the most practical technique for computing conditional expected values, the Gaussian
approximation, is not going to work well when the curvature in the volatility smile is ex-
plained by the local volatility function, rather than by the stochastic volatility component.
Sometimes it is beneficial to apply approximations that are not Gaussian, but linked to
them. For example, a lognormal approximation can be used in the place of a Gaussian
one. Then, the approximation is based on the known expression for E (X|Y ) when (X,Y )
are jointly lognormal. Other twists involving “mapping functions” Y → f (Y ) are also
possible.

5. PARAMETER AVERAGING
Parameter averaging methods are extensively described elsewhere (see Piterbarg (2005c),
Piterbarg (2005b)). Here, we present a short summary of the results.
Consider a stochastic volatility model of the form
p
(15) dz (t) = θ (1 − z (t)) dt + γ (t) z (t) dV (t) ,
p
dS (t) = (β (t) S (t) + (1 − β (t)) S (0)) σ (t) z (t) dW (t) ,
z (0) = 1, S (0) = S0 .

with time-dependent volatility σ (t) , skew β (t) , and volatility of variance γ (t). Fix a
particular T > 0.
 MARKOVIAN PROJECTION 9

Theorem 5.1. Prices of European options in the model (15) with expiry T > 0 are well-
approximated by the European option prices in the following model with constant coeffi-
cients,
p
(16) dz (t) = θ (1 − z (t)) dt + η z (t) dV (t) ,
p
dS (t) = (bS (t) + (1 − b) S (0)) λ z (t) dW (t) ,
where
RT 2
2γ (t) ρ (t) dt
η = 0 RT ,
0 ρ (t) dt
RT
0 β (t) v2 (t) σ 2 (t) dt
b= RT ,
0 v2 (t) σ 2 (t) dt
and λ is a solution to
 00 00
g (ζ ) T 2
  Z T 
2 g (ζ )
Z
E exp σ (t) z (t) dt = E exp λ 0 z (t) dt .
g0 (ζ ) 0 g (ζ ) 0
Here
Z T Z T
ρ (r) = ds dt σ 2 (t) σ 2 (s) e−θ (t−s) e−2θ (s−r) ,
r s
eθ s − e−θ s
Z t Z t
2 2 2 −θt
v (t) = σ (s) ds + η e σ 2 (s) ds,
0 0 2θ
S0 √


g (x) = 2Φ b x/2 − 1 ,
b
Z T
ζ= σ 2 (t) dt,
0

and Φ (·) is the standard Gaussian cumulative distribution function.

Proof. See Piterbarg (2005c), Piterbarg (2005b).
Remark 5.1. Efficient numerical methods to value European options in the model (16)
are available, see Andersen and Andreasen (2002).
Remark 5.2. Since the parameter averaging methods are especially well developed for
the model of the type (15), in the applications of the Markovian projection method it is
beneficial to try to reduce the dynamics of the quantity of interest, S (·), to this particu-
lar form. One should realize, however, that the only restriction imposed by (15) is the
particular form of the dynamics of the stochastic variance z (·) (a mean reverting square
root process). The form of the local volatility component of S (·) can be seen to be just a
linearization of whatever the actual local volatility function is used.
Remark 5.3. The averaging results have been derived in the case of zero asset-volatility
correlation, hdW, dV i = 0. Experience shows, however, that the approximations work
very well even if the correlation is not zero.
10 VLADIMIR PITERBARG

6. E XAMPLES
In this section we consider a couple of examples of varying degrees of complexity,
to demonstrate typical applications of the method we developed. The reader will note
that in the examples, the use of the parameter averaging methods is not required, as the
approximate equations already have constant coefficients. This is done on purpose, as the
applications of the PA methods have been extensively exposed elsewhere, and we try to
simplify the presentation on the account of clarity.

6.1. Example 1: Options on an index in local volatility models. Consider a collection

of N assets S (t) = (S1 (t) , . . . , SN (t))> each driven by its own local volatility model
dSn (t) = ϕn (Sn (t)) dWn (t) ,
n = 1, . . . , N,
with a vector of Brownian motions (W1 (t) , . . . ,WN (t)) with an N × N correlation matrix
dWi , dW j = ρi j , i, j = 1, . . . , N.
We restrict the set of local volatility functions ϕn (x) to the set of functions well-described
by its value and the first order derivative at the forward. In effect, we require ϕn (x) to be
linear, or close to linear (such as CEV). As always, this is not really a restriction as using
more complicated local volatility functions is not a good idea anyway.
Define an index S (t) by
N
S (t) = ∑ wnSn (t) = w>S (t) .
n=1

An index is a weighted average of the assets. We would like to derive an approximation

formula to (European) options on the index S (t). This problem was considered first in
Avellaneda et al. (2002).
Applying Ito’s lemma to S (t) we obtain that
N
dS (t) = ∑ wnϕn (Sn (t)) dWn (t) .
n=1

If we define
N
2
σ (t) = ∑ wn wm ϕn (Sn (t)) ϕm (Sm (t)) ρnm ,
n,m=1

1 N
dW (t) = ∑ wnϕn (Sn (t)) dWn (t) ,
σ (t) n=1
then
dS (t) = σ (t) dW (t) ,
and, by computing its quadratic variation, we obtain that dW is a standard Brownian
motion. The process σ (t) is, of course, a complicated function of S (t). By the results of
 MARKOVIAN PROJECTION 11

Section 2 and, in particular, Theorem 2.1, for the purposes of European option valuation,
we can replace this SDE with a simple one,
dS (t) = ϕ (t, S (t)) dW (t)
where
ϕ 2 (t, x) = E σ 2 (t) S (t) = x .

We have,
N
ϕ 2 (t, x) = ∑ wn wm ρnm E (ϕn (Sn (t)) ϕm (Sm (t))| S (t) = x) .
n,m=1

By expanding each term ϕn (Sn (t)) ϕm (Sm (t)) to the first order around the forward (a
good approximation by near-linearity assumption) we obtain,
ϕn (Sn (t)) ϕm (Sm (t)) ≈ pn pm
+ pn qm (Sm (t) − Sm (0))
+ pm qn (Sn (t) − Sn (0)) ,
so that
N
(17) ϕ 2 (t, x) ≈ ∑ wn wm ρnm pn pm (1 + χn (x) + χm (x)) ,
n,m=1

where
qn
(18) χn (x) = E (Sn (t) − Sn (0)| S (t) = x) ,
pn
and
pn = ϕn (Sn (0)) , qn = ϕn0 (Sn (0)) .
To compute the conditional expected value E (Sn (t) − Sn (0)| S (t) = x) we apply the Gauss-
ian approximation,


(19) E (Sn (t) − Sn (0)| S (t) = x) ≈ E S̄n (t) − Sn (0) S̄ (t) = x ,
where
d S̄n (t) = pn dWn (t) ,
d S̄ (t) = p dW̄ (t) ,
and
N
(20) 2
p = ∑ wn wm pn pm ρnm ,
n,m=1

1 N
dW̄ (t) = ∑ wn pn dWn (t) .
p n=1
12 VLADIMIR PITERBARG

In this model,

S̄n (t) , S̄ (t)
(21) E S̄n (t) − Sn (0) S̄ (t) = x = (x − S (0))
S̄ (t) , S̄ (t)
pn
= ρn (x − S (0)) ,
p
where
(22) ρn , hdW̄ (t) , dWn (t)i /dt
1 N
= ∑ wm pmρnm.
p m=1
Combining (18), (19), (21) we obtain
ρn
χn (x) = qn (x − S (0)) .
p
Substituting it into (17), we obtain the following proposition.
Proposition 6.1. For the purposes of European option pricing, the dynamics of the index
S (t) can be approximated with the following SDE,
dS (t) = ϕ (S (t)) dW (t) ,
where ϕ (x) is such that
ϕ (S (0)) = p, ϕ 0 (S (0)) = q,
with p given by (20),
1 N qn ρn + qm ρm
(23) q, 2 ∑
p n,m=1
wn wm pn pm ρnm
2
,

and {ρn }’s given by (22).

Remark 6.1. Note that in our approximation, ϕ (x) does not depend on t; this will not
always be the case. If the local volatility function depends on t, then parameter averaging
methods can be used to find the “‘aggregate” volatility and skew.
Remark 6.2. The parameter q is the slope of the local volatility function for the index.
Note that the expression (23) represents it as a weighted average of the slopes of the
local volatility functions of the index’s components qn , and correlations ρn between the
Brownian motions driving the components and the index. Another expression can be
easily derived from (23) by manipulating the order of summation,
1 N
(24) q, ∑ wn pnρn2qn.
p n=1
Remark 6.3. The approximations we develop work best when all qn ’s are “not very dif-
ferent”. This would be the case when, for example, a swap rate is represented as a
weighted average of contiguous Libor rates, and we can expect their skews to be similar.
In other cases, when qn ’s can be significantly different, as for example, in an FX rate
basket, non-linear approximations of Avellaneda et al. (2002) have better accuracy (but
are more complicated).
 MARKOVIAN PROJECTION 13

6.2. Example 2: Spread options with stochastic volatility. Consider the problem of
valuing spread options, ie (European) options on S (t) , S1 (t) − S2 (t) , where each of
the assets S1 , S2 follows its own stochastic volatility process (as pointed out in Piterbarg
(2006b), spread options are very sensitive to stochastic variance decorrelation, an often
overlooked fact). In particular we define
p
(25) dSi (t) = ϕi (Si (t)) zi (t) dWi (t) ,
p
dzi (t) = θ (1 − zi (t)) dt + ηi zi (t) dW2+i (t) , zi (0) = 1,
i = 1, 2,
with the correlations given by
dWi (t) , dW j (t) = ρi j , i, j = 1, . . . , 4.
The spread S (t) is the underlying of interest; the first step is to write an SDE for it. We
have,
(26) dS (t) = σ (t) dW (t) ,
where
(27) σ 2 (t) = (ϕ1 (S1 (t)) u1 (t))2 − 2 (ϕ1 (S1 (t)) u1 (t)) (ϕ2 (S2 (t)) u2 (t)) ρ12
+ (ϕ2 (S2 (t)) u2 (t))2 ,
1
dW (t) = (ϕ1 (S1 (t)) u1 (t) dW1 (t) − ϕ2 (S2 (t)) u2 (t) dW2 (t)) ,
σ (t)
p
ui (t) = zi (t), i = 1, 2.
Using the same notations as in the previous example, we set
pi = ϕi (Si (0)) , qi = ϕi0 (Si (0)) , i = 1, 2.
As explained previously, directly replacing σ 2 (t) (in (26)) with its expected value con-
ditional on S (·) is not a good idea when stochastic volatility is involved. Per Corollary
3.1 and the discussion in Section 3, we shall try to find a stochastic volatility process z (·)
such that the curvature of the smile of the spread S (·) is explained by it, and the local
volatility function is only used to induce the volatility skew.
To identify a suitable candidate for z (·), let us consider what the expression for σ 2 (t)
would be if the local volatility functions for the components, ϕi (x) , i = 1, 2, were constant
functions. We see that in this case, the expression for σ 2 (t) in (27) would not involve the
processes Si (·) , i = 1, 2 and thus would be a good candidate for the stochastic variance
process. Hence we define (the division by p2 , σ 2 (0) is to preserve the scaling z (0) = 1)
1  p 
(28) z (t) = 2 p21 z1 (t) − 2p1 p2 ρ12 z1 (t) z2 (t) + p22 z2 (t) ,
p
where

1/2
(29) p = σ (0) = p21 − 2p1 p2 ρ12 + p22 .
14 VLADIMIR PITERBARG

Proposition 6.2. For the purposes of European option pricing, the dynamics of the spread
S (t) can be approximated with the following SDE,
p
dS (t) = ϕ (S (t)) z (t) dW (t) ,
where z (·) given by (28), the function ϕ (·) satisfies
ϕ (S (0)) = p, ϕ 0 (S (0)) = q,
with
1
p1 ρ12 q1 − p2 ρ22 q2 ,


q,
p
and ρ1 , ρ2 are defined by
1
(30) ρ1 = (p1 − p2 ρ12 ) ,
p
1
ρ2 = (p1 ρ12 − p2 ) .
p
Proof. By Corollary 3.1, the right ϕ (x) id given by
E σ 2 (t) S (t) = x


2
(31) ϕ (t, x) = .
E (z (t)| S (t) = x)
The expression for E σ 2 (t) S (t) = x is a linear combinations of the conditional ex-


pected values of the terms


ϕi (Si (t)) ϕ j S j (t) ui (t) u j (t) .
Each one is approximated to the first order by
 
qi qj


pi p j 1 + (Si (t) − Si (0)) + S j (t) − S j (0) + (ui (t) − 1) + u j (t) − 1 .
pi pj
As in Example 1, the conditional expected values E (Si (t) − Si (0)| S (t) = x) can be easily
computed in the Gaussian approximation. The same idea can be applied to E (ui (t) − 1| S (t) = x) .
In particular,


E (Si (t) − Si (0)| S (t) = x) ≈ E S̄i (t) − S̄i (0) S̄ (t) = x


E (ui (t) − 1| S (t) = x) ≈ E ūi (t) − 1| S̄ (t) = x ,
where (dt terms are ignored, although they may be included for more accurate approxi-
mations)
(32) d S̄ (t) = p dW̄ (t) ,
d S̄i (t) = pi dWi (t)
ηi
d ūi (t) = dW2+i (t) ,
2
1
dW̄ (t) = (p1 dW1 (t) − p2 dW2 (t)) ,
p
 MARKOVIAN PROJECTION 15

and

pi ρ i
E S̄i (t) − S̄i (0) S̄ (t) = x = (x − S (0)) ,
p

ηi ρ2+i
E ūi (t) − 1| S̄ (t) = x = (x − S (0)) ,
2p
i = 1, 2.
Here, as in the previous example, we have denoted ρi , hdW̄ (t) , dWi (t)i /dt, i = 1, . . . , 4,
so that (see (30)),
1
(33) ρi = (p1 ρi1 − p2 ρi2 ) , i = 1, . . . , 4.
p
Under the linear approximations
E σ 2 (t) S (t) = x


 
2 x − S (0) x − S (0)
= p1 1 + 2q1 ρ1 + η1 ρ3
p p
   
x − S (0) η1 ρ3 + η2 ρ4 x − S (0)
− 2p1 p2 ρ12 1 + (q1 ρ1 + q2 ρ2 ) +
p 2 p
 
x − S (0) x − S (0)
+ p22 1 + 2q2 ρ2 + η2 ρ 4
p p
Then,
 
2

2 x − S (0)
(34) E σ (t) S (t) = x = p 1 + 2qu ,
p
where
1  2
qu = p (2q1 ρ1 + η1 ρ3 )
2p2 1
  
η1 ρ3 + η2 ρ4
+ p22 (2q2 ρ2 + η2 ρ4 ) .

−2p1 p2 ρ12 (q1 ρ1 + q2 ρ2 ) +
2
Furthermore,
 
2 x − S (0)
E (z (t)| S (t) = x) = p1 1 + η1 ρ3 /p2
p
   
η1 ρ3 + η2 ρ4 x − S (0)
− 2p1 p2 ρ12 1 + /p2
2 p
 
2 x − S (0)
+ p2 1 + η2 ρ4 /p2 ,
p
so that
x − S (0)
(35) E (z (t)| S (t) = x) = 1 + 2ql ,
p
   
1 2 η1 ρ3 + η2 ρ4 2
ql = 2 p1 η1 ρ3 − 2p1 p2 ρ12 + p2 η2 ρ4 .
2p 2
16 VLADIMIR PITERBARG

The quantities ϕ (S (0)) , ϕ 0 (S (0)) are computed from the equation (31) and the expres-
sions (34), (35) for the numerator and the denominator, as we note that
q = qu − ql .

To effectively apply the result of Proposition 6.2, we still need to derive a simple form
for the stochastic variance process z (·) , defined by (28). Ideally, we should try to write
the SDE for it in the same form as the SDEs for zi ’s in (25). This cannot be done exactly;
so we approximate.
Proposition 6.3. For the purposes of European option pricing, the dynamics of the spread
S (t) can be approximated with the following SDE,
p
dS (t) = ϕ (S (t)) z (t) dW (t) ,
 γ  p
dz (t) = θ 1 + − z (t) dt + η̄ z (t)d B̄ (t) ,
θ
where ϕ (x) is given in Proposition 6.2,
2 1  2 2

η̄ = 2 (p1 η1 ρ1 ) − 2 (p1 η1 ρ1 ) (p2 η2 ρ2 ) ρ34 + (p2 η2 ρ2 ) ,
p
1
d B̄ (t) = (p1 η1 ρ1 dW3 (t) − p2 η2 ρ2 dW4 (t)) ,
η̄ p
and γ is given by
p1 p2 ρ12 2
η1 − 2η1 η2 ρ34 + η22 .


(36) γ, 2
4p
The (approximation of the) correlation between dW and d B̄ is given by
1 2 2


hdW, d B̄i = p 1 η 1 ρ1 ρ13 − p 1 p2 (η 1 ρ1 ρ 23 + η 2 ρ 2 ρ14 ) + p2 η2 ρ2 ρ 24 .
p2 η̄
Proof. For z (·) defined by (28)
dz (t) = δ1 (t) dt + δ2 (t) dt + δ3 (t) dt + η (t) dB (t) ,
where
s !
p2 p2 z2 (t)
δ1 (t) = θ 21 1 − ρ12 (1 − z1 (t)) ,
p p1 z1 (t)
s !
p22 p1 z1 (t)
δ2 (t) = θ 1 − ρ12 (1 − z2 (t)) ,
p2 p2 z2 (t)
s s !
p1 p2 ρ12 z2 (t) 2 z1 (t) 2
δ3 (t) = η − 2η1 η2 ρ34 + η ,
4p2 z1 (t) 1 z2 (t) 2
 MARKOVIAN PROJECTION 17

and

η 2 (t) = v21 (t) + 2ρ34 v1 (t) v2 (t) + v22 (t) ,

1
dB (t) = (v1 (t) dW3 (t) + v2 (t) dW4 (t)) ,
η (t)
p21 p
 
p2 p
v1 (t) = η1 2 z1 (t) − ρ12 z2 (t) ,
p p1
p22 p
 
p1 p
v2 (t) = η2 2 z2 (t) − ρ12 z1 (t) .
p p2

The expressions are quite complicated, and also not “closed” in z (·) . Theorem 2.1 can
be applied again, now to the process for z (·) . The motivation for that is that the curvature
of the volatility smile (of options on S (·)) is driven by the variance of the stochastic vari-
ance, and that is preserved under the Markovian projection of z (·).pWere we do that, we

would need to compute conditional expected values of the type E zi (t) z j (t) z (t) = x
p

and E zi (t) /z j (t) z (t) = x , for which we would apply the Gaussian approximation
method. For the purposes of this example we, however, employ a much more ad-hoc (and
simpler) approximations, in which we
p p p
• replace qz1 (t),q z2 (t) with z (t);
• replace zz2 (t) (t) ,
z1 (t)
z (t) with 1.
1 2

Under this approximations,

δ1 (t) + δ2 (t) = θ (1 − z) ,
δ3 (t) = γ,
p
η (t) = η̄ z (t),
B (t) = B̄ (t) ,

where γ, η̄, and B̄ (·) are given in the statement of the proposition. The proposition
follows.

Test results show that the approximation is quite accurate, although it does deteriorate
for longer expiries. The accuracy can be improved by calculating various terms qmore
accurately. For example, let us demonstrate how this can be done for the terms zz2 (t) ,
q 1 (t)
z1 (t) p
z2 (t) , zi (t) z j (t). The results of Proposition 6.3 were obtained by replacing them with
1. A more accurate approximation is obtained by calculating their expected values in a
log-normal approximation.
Recall that
ηi2 ρ2+i,2+i  −2θt

(37) Ezi (t) = 1, Var (zi (t)) = 1−e , i = 1, 2.
2θ
18 VLADIMIR PITERBARG

The covariance between z1 (t) and z2 (t) can be approximated by2

η1 η2 ρ34  −2θt

(38) covar (z1 (t) , z2 (t)) = 1−e .
2θ
Let us approximate
 
1
zi (t) ≈ exp ξi − Var (ξi (t)) i = 1, 2,
2
where the Gaussian vector (ξ1 , ξ2 ) has zero mean and the variance-covariance matrix
such that the second-order moments of (z1 (t) , z2 (t)) given by (37), (38) are recovered.
the latter can be achieved by inverting the moments relationships
Var (zi (t)) = exp (Var (ξi )) − 1, i = 1, 2,
covar (z1 (t) , z2 (t)) = exp (covar (ξ1 , ξ2 )) − 1.
Then,
 
p 1
E z1 (t) z2 (t) = exp − Var (ξ1 − ξ2 )
8
 
p 1
E z1 (t) /z2 (t) = exp − (Var (ξ1 ) + 2covar (ξ1 , ξ2 ) − 3Var (ξ2 )) ,
8
and so on.
6.3. Example 3: Fitting the market with a stochastic volatility model. While the fo-
cus of this paper is on approximating more complicated models with simpler ones for
calibration purposes, direct calibration of stochastic volatility models to the market is
another possible application of the techniques we use.
Let S (t) be the value of a given market variable (for example, an equity). Let S0
be the spot value. Suppose market-implied Black volatilities σmkt (T, K) are known for
all expiries T and strikes K. These can be easily converted (see (5)) into the market-
implied Dupire’s local volatility bmkt (t, x) , ie a local volatility function such that the
market European option prices are consistent with the model
dY (t) = bmkt (t,Y (t)) dW (t) , Y (0) = S0 .
Suppose we would like to construct a stochastic volatility model
p
(39) dS (t) = b (t, S (t)) z (t) dW (t) , S (0) = S0 ,
consistent with market European option prices. As follows from Theorem 2.1, the local
volatility function b (t, x) is then given by
b2mkt (t, x)
(40) b2 (t, x) = .
E (z (t)| S (t) = x)
As mentioned before, the challenge of computing E (z (t)| S (t) = x) makes the method
difficult to implement in practice. Numerical methods are often applied to compute the
conditional expected value in a forward Kolmogorov PDE for (S (t) , z (t)) . Instead, we
would like to discuss a number of approaches to obtain (analytic) approximations.
2Usiang Gaussian approximation to the dynamics of z , z
1 2
 MARKOVIAN PROJECTION 19

A direct Gaussian-approximation-based approach (see Section 4) is possible, but is

likely to be inaccurate in this case if bmkt (t, x) has a high degree of convexity. What we
would like to explore instead are methods based on comparing (39) with another stochas-
tic volatility model, a model in which European option prices can be cheaply computed3.
Suppose we have such a “proxy” model, as defined by the (known) local volatility func-
tion b̃ (t, x), for the underlying X (t) and (the same) stochastic variance process z (t) ,
p
(41) dX (t) = b̃ (t, X (t)) z (t) dW (t) , X (0) = S0 .
Typically b̃ (t, x) is a simple function, most common being b̃ (t, x) = β (“normal” case),
b̃ (t, x) = νx (“log-normal” case), or b̃ (t, x) = νx + β (“shifted log-normal” case). Let
σproxy (T, K) be the implied Black volatilities of European options from this model. By
construction, σproxy (T, K) are cheap to compute. Then, they can be turned into the
“proxy” Dupire’s local volatility bproxy (t, x) (by the means of (5)). Then, rewriting (40),

b2proxy (t, x)
(42) E (z (t)| X (t) = x) = .
b̃2 (t, x)
Hence, having a stochastic volatility model which cheaply-computable European option
prices allows us to compute the conditional expected values E (z (t)| X (t) = x) easily. One
way to take advantage of this observation is to combine (40) and (42),
Proposition 6.4. The local volatility function b (t, x) that makes the model (39) consistent
with the market is given by

E (z (t)| X (t) = x) /12

 
bmkt (t, x)
b (t, x) = b̃ (t, x) × × ,
bproxy (t, x) E (z (t)| S (t) = x)

where X (t) follows the “proxy” model (41) with a known local volatility function b̃ (t, x) .
Remark 6.4. The local volatility function b (t, x) that we are trying to find is given by the
“proxy” local volatility b̃ (t, x) times two corrections, the ratio of Dupire’s local volatili-
mkt (t,x) E( z(t)|X(t)=x)
ties bbproxy (t,x) , and the ratio of conditional expected values E( z(t)|S(t)=x) .

bmkt (t,x)
The ratio bproxy (t,x) is, as discussed, cheap to compute. Approximating

E (z (t)| X (t) = x)
≈ 1,
E (z (t)| S (t) = x)
we obtain the following useful corollary.
Corollary 6.1. Approximately,
bmkt (t, x)
b (t, x) ≈ b̃ (t, x) × .
bproxy (t, x)

3This approach has originally been suggested in Forde (2006).

20 VLADIMIR PITERBARG

To obtain a more sophisticated approximation, we can look for an (approximate) func-

tional relationship between X (t) and S (t) . Denote
Z x
dy
(43) h (t, x) = ,
x0 b (t, y)
Z x
dy
h̃ (t, x) = ,
x0 b̃ (t, y)

H (t, x) = h̃−1 (t, h (t, x)) ,

where h̃−1 (t, x) is the inverse of h̃ (t, x) in the second (ie, x) argument. Furthermore,
denote
X̃ (t) = H (t, S (t)) .
Then
∂
d X̃ (t) = H (t, x) dS (t) + . . . dt
∂x x=S(t)

p
= b̃ t, X̃ (t) z (t) dW (t) + . . . dt.
We see that X̃ (t) and X (t) have the same diffusion coefficient. This inspires an approxi-
mation
(44) X (t) ≈ H (t, S (t)) ,
that leads to the following result.
Proposition 6.5. The local volatility function b (t, x) that makes the model (39) approxi-
mately consistent with the market is given by
bmkt (t, x)
(45) b (t, x) = b̃ (t, H (t, x)) ,
bproxy (t, H (t, x))
with H (·, ·) given by (43).
Proof. By (44),
E (z (t)| S (t) = x) = E (z (t)| H (t, S (t)) = H (t, x))
≈ E (z (t)| X (t) = H (t, x)) .
From (42),
b2proxy (t, H (t, x))
E (z (t)| S (t) = x) = ,
b̃2 (t, H (t, x))
and the result follows from (40).
Remark 6.5. Note that H (·, ·) depends on the (unknown) function b (t, x) , hence (45)
should be treated as a functional equation for b (t, x) .
In the small-time limit, the relationship between a local and implied volatilities are
much simpler than implied by (5). The following corollary is well-known (see for exam-
ple Andersen and Brotherton-Ratcliffe (2001)).
 MARKOVIAN PROJECTION 21

Corollary 6.2. Let X (·) follow the time-independent local volatility model
dX (t) = ϕ (X (t)) dW (t) , X (0) = X0 .
Then the implied Black volatility is given to the first order in t by

du −1
Z K 
σ (K) = log (K/X0 ) .
X0 ϕ (u)

By using this simple (harmonic average) relationship between the two, Propositions 6.4
and 6.5 can be restated in terms of the implied, rather than local, volatilities. We denote
b (x) , limt→0+ b (t, x) , etc. Proposition 6.5 in the short-time limit was first proved in
Forde (2006),
Corollary 6.3. In the short-time limit,
0  
log (x/S0 ) 0
 
log (x/S0 )
b (x) = b̃ (H (x)) .
σproxy (H (x)) σmkt (x)
One question remains unanswered, and that is the choice of the “proxy” model. If z (·)
follows the standard mean-reverting square root process (as in the Heston model), then
the following models admit closed-form expressions for European options,
b̃ (t, x) = νx
(the original Heston model), and
b̃ (t, x) = νx + β ,
the shifted Heston model, as previously explained. Moreover, using the parameter av-
eraging techniques, this could be extended to a shifted-lognormal Heston model with
time-dependent coefficients,
b̃ (t, x) = ν (t) x + β (t) .
The last model, is, probably, the best choice for the proxy model. The time-dependent
coefficients can be chosen to make the “proxy” model as similar as possible to the model
for X (·), thus improving the quality of various approximations made. Ideally we should
use
β (t) = b (t, S0 ) ,
∂
ν (t) = b (t, S0 ) ,
∂x
the first-order approximation to the local volatility b (t, x) “along the forward”. The value
and the derivative of the local volatility b (t, x) are of course unknown a-priori; one can
envision various approximations or, conceivably, an iterative procedure where the local
volatility approximation obtained on step n is used to define the proxy local volatility for
step n + 1.
22 VLADIMIR PITERBARG

7. C ONCLUSIONS
We have developed a generic method for obtaining closed-form approximations to Eu-
ropean options and demonstrated its usage in a number of simple examples. The accuracy
of the approximations critically depends on how accurately the conditional expected val-
ues can be estimated. In particular, any improvements over Gaussian-based methods will
have an immediate effect on the European option approximations.
R EFERENCES
Andersen, L. B. and Andreasen, J. (2002). Volatile volatilities. Risk, 15(12): 163–168.
Andersen, L. B. and Brotherton-Ratcliffe, R. (2001). Extended Libor market models with
stochastic volatility. Working paper.
Andreasen, J. (2005). Back to the future. Risk Magazine, 18(9): 104–109.
Avellaneda, M., Boyer-Olson, D., Busca, J., and Friz, P. (2002). Reconstructing volatility.
Risk, 15(10).
Dupire, B. (1994). Pricing with a smile. Risk, 7(1).
Dupire, B. (1997). A unified theory of volatility. Banque Paribas working paper.
Forde, M. (2006). Calibrating local-stochastic volatility models. University of Bristol
working paper, http://www.maths.bris.ac.uk/ mamsf/.
Gyöngy, I. (1986). Mimicking the one-dimensional distributions of processes having an
Ito differential. Prob. Th. Rel. Fields, 71: 501 – 516.
Karatzas, I. and Shreve, S. E. (1997). Brownian Motion and Stochastic Calculus. Springer.
Piterbarg, V. V. (2004). A stochastic volatility forward Libor model with a term structure
of volatility smiles. SSRN Working paper.
Piterbarg, V. V. (2005a). A multi-currency model with FX volatility skew. SSRN working
paper.
Piterbarg, V. V. (2005b). Stochastic volatility model with time-dependent skew. Applied
Mathematical Finance, 12(2): 147–185.
Piterbarg, V. V. (2005c). Time to smile. Risk Magazine, 18(5): 71–75.
Piterbarg, V. V. (2006a). Smiling hybrids. Risk Magazine, 19(5): 66–71.
Piterbarg, V. V. (2006b). Towards a multi-stochastic volatility model for CMS spread
exotics. Global Derivatives Trading and Risk Management Conference, Paris.
Savine, A. (2000). A theory of volatility. Reech working paper,
http://www.reech.com/company/research.papers/.
BARCLAYS C APITAL , 5 T HE N ORTH C OLONNADE , E14 4BB, L ONDON
E-mail address: [email protected]