Approximation of the distribution of a stationary Markov process with
  application to option pricing
Approximation of the distribution of a stationary Markov process with application to option pricing
Bernoulli 15(1), 2009, 146–177
DOI: 10.3150/08-BEJ142
Approximation of the distribution of a
stationary Markov process with application to
option pricing
GILLES PAGÈS1 and FABIEN PANLOUP2
Laboratoire de Probabilités et Modèles Aléatoires, UMR 7599, Université Paris 6, Case 188, 4
pl. Jussieu, F-75252 Paris Cedex 5. E-mail: gpa@ccr.jussieu.fr
Laboratoire de Statistiques et Probabilités, Université Paul Sabatier & INSA Toulouse, 135,
Avenue de Rangueil, 31077 Toulouse Cedex 4. E-mail: fpanloup@insa-toulouse.fr
We build a sequence of empirical measures on the space D(R+,R
d) of Rd-valued cadlag functions
on R+ in order to approximate the law of a stationary R
d-valued Markov and Feller process (Xt).
We obtain some general results on the convergence of this sequence. We then apply them to
Brownian diffusions and solutions to Lévy-driven SDE’s under some Lyapunov-type stability
assumptions. As a numerical application of this work, we show that this procedure provides an
efficient means of option pricing in stochastic volatility models.
Keywords: Euler scheme; Lévy process; numerical approximation; option pricing; stationary
process; stochastic volatility model; tempered stable process
1. Introduction
1.1. Objectives and motivations
In this paper, we deal with an Rd-valued Feller Markov process (Xt) with semigroup
(Pt)t≥0 and assume that (Xt) admits an invariant distribution ν0. The aim of this work is
to propose a way to approximate the whole stationary distribution Pν0 of (Xt). More pre-
cisely, we want to construct a sequence of weighted occupation measures (ν(n)(ω,dα))n≥1
on the Skorokhod space D(R+,R
d) such that ν(n)(ω,F )
n→+∞−→
F (α)Pν0(dα) a.s. for a
class of functionals F :D(R+,R
d) which includes bounded continuous functionals for the
Skorokhod topology.
One of our motivations is to develop a new numerical method for option pricing in sta-
tionary stochastic volatility models which are slight modifications of the classical stochas-
tic volatility models, where we suppose that the volatility evolves under its stationary
regime.
This is an electronic reprint of the original article published by the ISI/BS in Bernoulli,
2009, Vol. 15, No. 1, 146–177. This reprint differs from the original in pagination and
typographic detail.
1350-7265 c© 2009 ISI/BS
http://arxiv.org/abs/0704.0335v3
http://isi.cbs.nl/bernoulli/
http://dx.doi.org/10.3150/08-BEJ142
mailto:gpa@ccr.jussieu.fr
mailto:fpanloup@insa-toulouse.fr
http://isi.cbs.nl/BS/bshome.htm
http://isi.cbs.nl/bernoulli/
http://dx.doi.org/10.3150/08-BEJ142
Approximation of the distribution of a stationary Markov process 147
1.2. Background and construction of the procedure
This work follows on from a series of recent papers due to Lamberton and Pagès ([12, 13]),
Lemaire ([14, 15]) and Panloup ([18, 19, 20]), where the problem of the approximation
of the invariant distribution is investigated for Brownian diffusions and for Lévy-driven
SDE’s.1 In these papers, the algorithm is based on an adapted Euler scheme with de-
creasing step (γk)k≥1. To be precise, let (Γn) be the sequence of discretization times:
Γ0 = 0, Γn =
k=1 γk for every n≥ 1, and assume that Γn → +∞ when n→+∞. Let
(X̄Γn)n≥0 be the Euler scheme obtained by “freezing” the coefficients between the Γn’s
and let (ηn)n≥1 be a sequence of positive weights such that Hn :=
k=1 ηk →+∞ when
k→+∞. Then, under some Lyapunov-type stability assumptions adapted to the stochas-
tic processes of interest, one shows that for a large class of steps and weights (ηn, γn)n≥1,
ν̄n(ω, f) :=
ηkf(X̄Γk−1)
n→+∞−→
f(x)ν0(dx) a.s., (1)
(at least)2 for every bounded continuous function f .
Since the problem of the approximation of the invariant distribution has been deeply
studied for a wide class of Markov processes (Brownian diffusions and Lévy-driven SDE’s)
and since the proof of (1) can be adapted to other classes of Markov processes under some
specific Lyapunov assumptions, we choose in this paper to consider a general Markov pro-
cess and to assume the existence of a time discretization scheme (X̄Γk)k≥0 such that (1)
holds for the class of bounded continuous functions. The aim of this paper is then to inves-
tigate the convergence properties of a functional version of the sequence (ν̄n(ω,dα))n≥1.
Let (Xt) be a Markov and Feller process and let (X̄t)t≥0 be a stepwise constant time
discretization scheme of (Xt) with non-increasing step sequence (γn)n≥1 satisfying
γn = 0, Γn :=
n→+∞−→ +∞. (2)
Letting Γ0 := 0 and X̄0 = x0 ∈Rd, we assume that
X̄t = X̄Γn ∀t ∈ [Γn,Γn+1[ (3)
and that (X̄Γn)n≥0 can be simulated recursively.
We denote by (Ft)t≥0 and (F̄t)t≥0 the usual augmentations of the natural filtrations
(σ(Xs,0≤ s≤ t))t≥0 and (σ(X̄s,0≤ s≤ t))t≥0, respectively.
1Note that computing the invariant distribution is equivalent to computing the marginal laws of the
stationary process (Xt) since ν0Pt = ν0 for every t≥ 0.
2The class of functions for which (1) holds depends on the stability of the dynamical system. In
particular, in the Brownian diffusion case, the convergence may hold for continuous functions with
subexponential growth, whereas the class of functions strongly depends on the moments of the Lévy
process when the stochastic process is a Lévy-driven SDE.
148 G. Pagès and F. Panloup
For k ≥ 0, we denote by (X̄(k)t )t≥0 the shifted process defined by
t := X̄Γk+t.
In particular, X̄
t = X̄t. We define a sequence of random probabilities (ν
(n)(ω,dα))n≥1
on D(R+,R
d) by
ν(n)(ω,dα) =
ηk1{X̄(k−1)(ω)∈dα},
where (ηk)k≥1 is a sequence of weights. For t ≥ 0, (ν(n)t (ω,dx))n≥1 will denote the se-
quence of “marginal” empirical measures on Rd defined by
t (ω,dx) =
ηk1{X̄(k−1)
(ω)∈dx}
1.3. Simulation of (ν(n)(ω,F ))
For every functional F :D(R+,R
d)→R, the following recurrence relation holds for every
n≥ 1:
ν(n+1)(ω,F ) = ν(n)(ω,F ) +
(F (X(n)(ω))− ν(n)(ω,F )). (4)
Then, if T is a positive number and F :D(R+,R
d) → R is a functional depending only
on the trajectory between 0 and T , (ν(n)(ω,F ))n≥1 can be simulated by the following
procedure.
Step 0. (i) Simulate (X̄
t )t≥0 on [0, T ], that is, simulate (X̄Γk)k≥0 for k =
0, . . . ,N(0, T ), where
N(n,T ) := inf{k ≥ n,Γk+1 − Γn > T }
= max{k ≥ 0,Γk − Γn ≤ T }, n≥ 0, T > 0.
Note that n 7→N(n, t) is an increasing sequence since (γn) is non-increasing, and that
ΓN(n,T ) − Γn ≤ T < ΓN(n,T )+1 − Γn.
(ii) Compute F ((X̄
t )t≥0) and ν
(1)(ω,F ). Store the values of (X̄Γk) for k =
1, . . . ,N(0, T ).
Step n (n≥ 1). (i) Since the values (X̄Γk)k≥0 are stored for k = n, . . . ,N(n− 1, T ),
simulate (X̄Γk)k≥0 for k =N(n−1, T )+1, . . . ,N(n,T ) in order to obtain a path of (X̄
on [0, T ].
(ii) Compute F ((X̄
t )t≥0) and use (4) to compute ν
(n+1)(ω,F ). Store the values of
(X̄Γk) for k = n+ 1, . . . ,N(n,T ).
Approximation of the distribution of a stationary Markov process 149
Remark 1. As shown in the description of the procedure, one generally has to store
the vector [X̄Γn , . . . , X̄ΓN(n,T) ] at time n. Since (γn) is a sequence with infinite sum that
decreases to 0, it follows that the size of this vector increases “slowly” to +∞. For
instance, if γn = Cn
−ρ with ρ ∈ (0,1), its size is of order nρ. However, it is important
to remark that even though the number of values to be stored tends to +∞, that is
not always the case for the number of operations at each step. Indeed, since X̄(n+1)
is obtained by shifting X̄(n), it is usually possible to use, at step n+ 1, the preceding
computations and to simulate the sequence (F (X̄(n)))n≥0 in a “quasi-recursive” way.
For instance, such remark holds for Asian options because the associated pay-off can be
expressed as a function of an additive functional (see Section 5 for simulations).
Before outlining the sequel of the paper, we list some notation linked to the spaces
D(R+,R
d) and D([0, T ],Rd) of cadlag Rd-valued functions on R+ and [0, T ], respectively,
endowed with the Skorokhod topology. First, we denote by d1 the Skorokhod distance
on D([0,1],Rd) defined for every α, β ∈D([0,1],Rd) by
d1(α,β) = inf
t∈[0,1]
|α(t)− β(λ(t))|, sup
0≤s<t≤1
λ(t)− λ(s)
where Λ1 denotes the set of increasing homeomorphisms of [0,1]. Second, for T > 0,
φT :D(R+,R
d) 7→D([0,1],Rd) is the function defined by (φT (α))(s) = α(sT ) for every s ∈
[0,1]. We then denote by d the distance on D(R+,R
d) defined for every α,β ∈D(R+,Rd)
d(α,β) =
e−t(1∧ d1(φt(α), φt(β))) dt. (6)
We recall that (D(R+,R
d), d) is a Polish space and that the induced topology is the usual
Skorokhod topology on D(R+,R
d) (see, e.g., Pagès [16]). For every T > 0, we set
σ(πu,0≤ u≤ s),
where πs :D(R+,R
d)→Rd is defined by πs(α) = α(s). For a functional F :D(R+,Rd)→
R, FT denotes the functional defined for every α ∈D(R+,Rd) by
FT (α) = F (α
T ) with αT (t) = α(t ∧ T ) ∀t≥ 0. (7)
Finally, we will say that a functional F :D(R+,R
d)→R is Sk-continuous if F is contin-
uous for the Skorokhod topology on D(R+,R
d) and the notation “
=⇒” will denote the
weak convergence on D(R+,R
In Section 2, we state our main results for a general Rd-valued Feller Markov process.
Then, in Section 3, we apply them to Brownian diffusions and Lévy-driven SDE’s. Section
4 is devoted to the proofs of the main general results. Finally, in Section 5, we complete
this paper with an application to option pricing in stationary stochastic volatility models.
150 G. Pagès and F. Panloup
2. General results
In this section, we state the results on convergence of the sequence (ν(n)(ω,dα))n≥1 when
(Xt) is a general Feller Markov process.
2.1. Weak convergence to the stationary regime
As explained in the Introduction, since the a.s. convergence of (ν
0 (ω,dx))n≥1 to the
invariant distribution ν0 has already been deeply studied for a large class of Markov
processes (Brownian diffusions and Lévy driven SDE’s), our approach will be to derive
the convergence of (ν(n)(ω,dα))n≥1 toward Pν0 from that of (ν
0 (ω,dx))n≥1 to the
invariant distribution ν0. More precisely, we will assume in Theorem 1 that
(C0,1): (Xt) admits a unique invariant distribution ν0 and
0 (ω,dx)
=⇒ ν0(dx) a.s.,
whereas in Theorem 2, we will only assume that
(C0,2): (ν
0 (ω,dx))n≥1 is a.s. tight on R
We also introduce three other assumptions, (C1), (C2) and (C3,ε), regarding the conti-
nuity in probability of the flow x 7→ (Xxt ), the asymptotic convergence of the shifted time
discretization scheme to the true process (Xt) and the steps and weights, respectively.
(C1): For every x0 ∈Rd, ǫ > 0 and T > 0,
limsup
0≤t≤T
|Xxt −Xx0t | ≥ ǫ
= 0. (8)
(C2): (X̄t) is a non-homogeneous Markov process and for every n≥ 0, it is possible to
construct a family of stochastic processes (Y
(n,x)
t )x∈Rd such that
(i) L(Y (n,x)) D(R+,R
= L(X̄(n)|X̄(n)0 = x);
(ii) for every compact set K of Rd, for every T ≥ 0,
0≤t≤T
|Y (n,x)t −Xxt |
n→+∞−→ 0 in probability. (9)
(C3,ε): For every n≥ 1, ηn ≤CγnHεn.
Remark 2. Assumption (C2) implies, in particular, that asymptotically and uniformly
on compact sets of Rd, the law of the approximate process (X̄(n)), given its initial value,
is close to that of the true process.
If there exists a unique invariant distribution ν0, the second part of (C2) can be relaxed
to the following, less stringent, assertion: for all ǫ > 0, there exists a compact set Aǫ ⊂Rd
such that ν0(A
ǫ)≤ ǫ and such that
0≤t≤T
|Y (n,x)t −Xxt |
n→+∞−→ 0 in probability. (10)
Approximation of the distribution of a stationary Markov process 151
This weaker assumption can some times be needed in stochastic volatility models like
the Heston model (see Section 5 for details).
The preceding assumptions are all that we require for the convergence of (ν(n)(ω,dα))n≥1
to Pν0 along the bounded Sk-continuous functionals, that is, for the a.s. weak conver-
gence on D(R+,R
d). However, the integration of non-bounded continuous functionals
F :D([0, T ],Rd)→ R will need some additional assumptions, depending on the stability
of the time discretization scheme and on the steps and weights sequences. We will sup-
pose that F is dominated (in a sense to be specified later) by a function V : Rd → R+
that satisfies the following assumptions for some s≥ 2 and ε < 1.
H(s, ε): For every T > 0,
(i) sup
0≤t≤T
Vs(Y (n,x)t )
≤CTVs(x),
(ii) sup
0 (V)<+∞,
(iii)
E[V2(X̄Γk−1)]<+∞,
∆N(k,T )
E[Vs(1−ε)(X̄Γk−1)]<+∞,
where T 7→CT is locally bounded on R+ and ∆N(k,T ) =N(k,T )−N(k− 1, T ).
For every ε < 1, we then set
K(ε) = {V ∈ C(Rd,R+),H(s, ε) holds for some s≥ 2}.
Remark 3. Apart from assumption (i), which is a classical condition on the finite time
horizon control, the assumptions in H(s, ε) strongly rely on the stability of the time
discretization scheme (and then, to that of the true process). More precisely, we will see
when we apply our general results to SDE’s that these properties are some consequences
of the Lyapunov assumptions needed for the tightness of (ν
0 (ω,dx))n≥1.
We can now state our first main result.
Theorem 1. Assume (C0,1), (C1), (C2) and (C3,ε) with ε ∈ (−∞,1). Then, a.s., for
every bounded Sk-continuous functional F :D(R+,R
d)→R,
ν(n)(ω,F )
n→+∞−→
F (α)Pν0 (dα), (11)
where Pν0 denotes the stationary distribution of (Xt) (with initial law ν0).
Furthermore, for every T > 0, for every non-bounded Sk-continuous functional
F :D(R+,R
d)→ R, (11) holds a.s. for FT (defined by (7)) if there exists V ∈ K(ε) and
152 G. Pagès and F. Panloup
ρ ∈ [0,1) such that
|FT (α)| ≤C sup
0≤t≤T
Vρ(αt) ∀α ∈D(R+,Rd). (12)
In the second result, the uniqueness of the invariant distribution is not required and
the sequence (ν
0 (ω,dx))n≥1 is only supposed to be tight.
Theorem 2. Assume (C0,2), (C1), (C2) and (C3,ε) with ε ∈ (−∞,1). Assume that
0 (ω,dx))n≥1 is a.s. tight on R
d. We then have the following.
(i) The sequence (ν(n)(ω,dα))n≥1 is a.s. tight on D(R+,R
d) and a.s., for ev-
ery convergent subsequence (nk(ω))n≥1, for every bounded Sk-continuous functional
F :D(R+,R
d)→R,
ν(nk(ω))(ω,F )
n→+∞−→
F (α)Pν∞(dα), (13)
where Pν∞ is the law of (Xt) with initial law ν∞ being a weak limits for (ν
0 (ω,dx))n≥1.
Furthermore, for every T > 0, for every non-bounded Sk-continuous functional
F :D(R+,R
d)→R, (13) holds a.s. for FT if (12) is satisfied with V ∈K(ε) and ρ ∈ [0,1).
(ii) If, moreover,
l≥k+1
|∆ηℓ|
n→+∞−→ 0, (14)
then ν∞ is necessarily an invariant distribution for the Markov process (Xt).
Remark 4. Condition (14) holds for a large class of steps and weights. For instance,
if ηn = C1n
−ρ1 and γn = C2n
−ρ2 with ρ1 ∈ [0,1] and ρ2 ∈ (0,1], then (14) is satisfied if
ρ1 = 0 or if ρ1 ∈ (max(0,2ρ2 − 1),1).
2.2. Extension to the non-stationary case
Even though the main interest of this algorithm is the weak approximation of the pro-
cess when stationary, we observe that when ν0 is known, the algorithm can be used to
approximate Pµ0 if µ0 is a probability on R
d that is absolutely continuous with respect
to ν0.
Indeed, assume that µ0(dx) = φ(x)ν0(dx), where φ :R
d → R is a continuous non-
negative function. For a functional F :D(R+,R
d)→ R, denote by Fφ the functional de-
fined on D(R+,R
d) by Fφ(α) = F (α)φ(α(0)).
Then, if ν(n)(ω,dα)
(Sk)⇒ Pν0(dα) a.s., we also have the following convergence: a.s., for
every bounded Sk-continuous functional F :D(R+,R
d)→R,
ν(n)(ω,Fφ)
n→+∞−→
Fφ(α)Pν0 (dα) =
F (α)Pµ0 (dα).
Approximation of the distribution of a stationary Markov process 153
3. Application to Brownian diffusions and
Lévy-driven SDE’s
Let (Xt)t≥0 be a cadlag stochastic process solution to the SDE
dXt = b(Xt−) dt+ σ(Xt−) dWt + κ(Xt−) dZt, (15)
where b :Rd → Rd, σ :Rd 7→Md,ℓ (set of d× ℓ real matrices) and κ :Rd 7→Md,ℓ are con-
tinuous functions with sublinear growth, (Wt)t≥0 is an ℓ-dimensional Brownian motion
and (Zt)t≥0 is an integrable purely discontinuous R
ℓ-valued Lévy process independent of
(Wt)t≥0 with Lévy measure π and characteristic function given for every t≥ 0 by
E[ei〈u,Zt〉] = exp
ei〈u,y〉 − 1− i〈u, y〉π(dy)
Let (γn)n≥1 be a non-increasing step sequence satisfying (2). Let (Un)n≥1 be a sequence
of i.i.d. random variables such that U1
=N (0, Iℓ) and let ξ := (ξn)n≥1 be a sequence of
independent Rℓ-valued random variables, independent of (Un)n≥1. We then denote by
(X̄t)t≥0 the stepwise constant Euler scheme of (Xt) for which (X̄Γn)n≥0 is recursively
defined by X̄0 = x ∈Rd and
X̄Γn+1 = X̄Γn + γn+1b(X̄Γn) +
γn+1σ(X̄Γn)Un+1 + κ(X̄Γn)ξn+1. (16)
We recall that the increments of (Zt) cannot be simulated in general. That is why we
generally need to construct the sequence (ξn) with some approximations of the true
increments. We will come back to this construction in Section 3.2.
As in the general case, we denote by (X̄(k))k≥0 and (ν
(n)(ω,dα))n≥1 the sequences of
associated shifted Euler schemes and empirical measures, respectively.
Let us now introduce some Lyapunov assumptions for the SDE. Let EQ(Rd) denote
the set of essentially quadratic C2-functions V :Rd → R∗+ such that limV (x) = +∞ as
|x| →+∞, |∇V | ≤C
V and D2V is bounded. Let a ∈ (0,1] denote the mean reversion
intensity. The Lyapunov (or mean reversion) assumption is the following.
(Sa): There exists a function V ∈ EQ(Rd) such that:
(i) |b|2 ≤CV a, Tr(σσ∗(x)) + ‖κ(x)‖2 |x|→+∞= o(V a(x));
(ii) there exist β ∈R and ρ > 0 such that 〈∇V, b〉 ≤ β − ρV a.
From now on, we separate the Brownian diffusions and Lévy-driven SDE cases.
3.1. Application to Brownian diffusions
In this part, we assume that κ= 0. We recall a result by Lamberton and Pagès [13].
Proposition 1. Let a ∈ (0,1] such that (Sa) holds. Assume that the sequence (ηn/γn)n≥1
is non-increasing.
154 G. Pagès and F. Panloup
(a) Let (θn)n≥1 be a sequence of positive numbers such that
n≥1 θnγn < +∞ and
that there exists n0 ∈N such that (θn)n≥n0 is non-increasing. Then, for every positive r,
θnγnE[V
r(X̄Γn−1)]<+∞.
(b) For every r > 0,
0 (ω,V
r)<+∞ a.s. (17)
Hence, the sequence (ν
0 (ω,dx))n≥1 is a.s. tight.
(c) Moreover, every weak limit of this sequence is an invariant probability for the SDE
(15). In particular, if (Xt)t≥0 admits a unique invariant probability ν0, then for every
continuous function f such that f ≤CV r with r > 0, limn→∞ ν(n)0 (ω, f) = ν0(f) a.s.
Remark 5. For instance, if V (x) = 1 + |x|2, then the preceding convergence holds for
every continuous function with polynomial growth. According to Theorem 3.2 in Lemaire
[14], it is possible to extend these results to continuous functions with exponential growth,
but it then strongly depends on σ. Further the conditions on steps and weights can be
less restrictive and may contain the case ηn = 1, for instance (see Remark 4 of Lamberton
and Pagès [13] and Lemaire [14]).
We then derive the following result from the preceding proposition and from Theorems
1 and 2.
Theorem 3. Assume that b and σ are locally Lipschitz functions and that κ = 0. Let
a ∈ (0,1] such that (Sa) holds and assume that (ηn/γn) is non-increasing.
(a) The sequence (ν(n)(ω,dα))n≥1 is a.s. tight on C(R+,Rd)3 and every weak limit
of (ν(n)(ω,dα))n≥1 is the distribution of a stationary process solution to (15). In par-
ticular, when uniqueness holds for the invariant distribution ν0, a.s., for every bounded
continuous functional F :C(R+,Rd)→R,
ν(n)(ω,F )
n→+∞−→
F (x)Pν0 (dx). (18)
(b) Furthermore, if there exists s ∈ (2,+∞) and n0 ∈N such that
∆N(k,T )
is non-increasing and
∆N(k,T )
<+∞, (19)
3C(R+,R
d) denotes the space of continuous functions on R+ with values in R
d endowed with the
topology of uniform convergence on compact sets.
Approximation of the distribution of a stationary Markov process 155
then, for every T > 0, for every non-bounded continuous functional F :C(R+,Rd)→ R,
(18) holds for FT if the following condition is satisfied:
∃r > 0 such that |FT (α)| ≤C sup
0≤t≤T
V r(αt) ∀α ∈ C(R+,Rd).
Remark 6. If ηn =C1n
−ρ1 and γn =C2n
−ρ2 with 0< ρ2 ≤ ρ1 ≤ 1, then for s ∈ (1,+∞),
(19) is fulfilled if and only if s > 1/(1− ρ1). It follows that there exists s ∈ (2,+∞) such
that (19) holds as soon as ρ1 < 1.
Proof of Theorem 3. We want to apply Theorem 2. First, by Proposition 1, assumption
(C0,2) is fulfilled and every weak limit of (ν
0 (ω,dx)) is an invariant distribution. Second,
it is well known that (C1) and (C2) are fulfilled when b and σ are locally Lispchitz
sublinear functions. Then, since (C3,ε) holds with ε = 0, (18) holds for every bounded
continuous functional F . Finally, one checks that H(s,0) holds with V := V r (r > 0).
It is classical that assumption (a) is true when b and σ are sublinear. Assumption (b)
follows from Proposition 1(b). Let θn,1 = ηn/(γnH
n) and θn,2 =∆N(n,T )/(γnH
n). Using
(19) and the fact that (ηn/γn) is non-increasing yields that (θn,1) and (θn,2) satisfy the
conditions of Proposition 1 (see (35) for details). Then, (iii) and (iv) of H(s,0) are
consequences of Proposition 1(a). This completes the proof. �
3.2. Application to Lévy-driven SDE’s
When we want to extend the results obtained for Brownian SDE’s to Lévy-driven SDE’s,
one of the main difficulties comes from the moments of the jump component (see Panloup
[18] for details). For simplification, we assume here that (Zt) has a moment of order
2p≥ 2, that is, that its Lévy measure π satisfies the following assumption with p≥ 1:
(H1p) :
|y|>1
π(dy)|y|2p <+∞.
We also introduce an assumption about the behavior of the moments of the Lévy measure
at 0:
(H2q) :
|y|≤1
π(dy)|y|2q <+∞, q ∈ [0,1].
This assumption ensures that (Zt) has finite 2q-variations. Since
|y|≤1
|y|2π(dy) is finite,
this is always satisfied for q = 1.
Let us now specify the law of (ξn) introduced in (16). When the increments of (Zt) can
be exactly simulated, we denote by (E) the Euler scheme and by (ξn,E) the associated
sequence
= Zγn ∀n≥ 1.
156 G. Pagès and F. Panloup
When the increments of (Zt) cannot be simulated, we introduce some approximated Euler
schemes (P) and (W) built with some sequences (ξn,P ) and (ξn,W ) of approximations of
the true increment (see Panloup [19] for more detailed presentations of these schemes).
In scheme (P),
=Zγn,n,
where (Z·,n)n≥1 a sequence of compensated compound Poisson processes obtained by
truncating the small jumps of (Zt)t≥0:
Zt,n :=
0<s≤t
∆Zs1{|∆Zs|>un} − t
|y|>un
yπ(dy) ∀t≥ 0, (20)
where (un)n≥1 is a sequence of positive numbers such that un → 0. We recall that
n→+∞−→ Z locally uniformly in L2 (see, e.g., Protter [21]).
As shown in Panloup [19], the error induced by this approximation is very large when
the local behavior of the small jumps component is irregular. However, it is possible to
refine this approximation by a Wienerization of the small jumps, that is, by replacing
the small jumps by a linear transform of a Brownian motion instead of discarding them
(see Asmussen and Rosinski [2]). The corresponding scheme is denoted by (W) with ξn,W
satisfying
= ξn,P +
γnQnΛn ∀n≥ 1,
where (Λn)n≥1 is a sequence of i.i.d. random variables, independent of (ξn,P )n≥1 and
(Un)n≥1, such that Λ1
=N (0, Iℓ) and (Qn) is a sequence of ℓ× ℓ matrices such that
n)i,j =
|y|≤uk
yiyjπ(dy).
We recall the following result obtained in Panloup [18] in our slightly simplified frame-
work.
Proposition 2. Let a ∈ (0,1], p≥ 1 and q ∈ [0,1] such that (H1p), (H2q) and (Sa) hold.
Assume that the sequence (ηn/γn)n≥1 is non-increasing. Then, the following assertions
hold for schemes (E), (P) and (W).
(a) Let (θn) satisfy the conditions of Proposition 1. Then,
n≥1 θnγnE[V
p+a−1(X̄Γn−1)]<
(b) We have
0 (ω,V
p/2+a−1)<+∞ a.s. (21)
Hence, the sequence (ν
0 (ω,dx))n≥1 is a.s. tight as soon as p/2+ a− 1> 0.
Approximation of the distribution of a stationary Markov process 157
(c) Moreover, if Tr(σσ∗)+ ‖κ‖2q ≤CV p/2+a−1, then every weak limit of this sequence
is an invariant probability for the SDE (15). In particular, if (Xt)t≥0 admits a unique
invariant probability ν0, for every continuous function f such that f = o(V
p/2+a−1),
limn→∞ ν
0 (ω, f) = ν0(f) a.s.
Remark 7. For schemes (E) and (P), the above proposition is a direct consequence of
Theorem 2 and Proposition 2 of Panloup [18]. As concerns scheme (W), a straightforward
adaptation of the proof yields the result.
Our main functional result for Lévy-driven SDE’s is then the following.
Theorem 4. Let a ∈ (0,1] and p≥ 1 such that p/2+ a− 1> 0 and let q ∈ [0,1]. Assume
(H1p), (H
q) and (Sa). Assume that b, σ and κ are locally Lipschitz functions. If, more-
over, (ηn/γn)n≥1 is non-increasing, then the following result holds for schemes (E), (P)
and (W).
(a) The sequence (ν(n)(ω,dα))n≥1 is a.s. tight on D(R+,R
d). Moreover, if
Tr(σσ∗) + ‖κ‖2q ≤CV p/2+a−1 or 1
l≥k+1
|∆ηℓ|
n→+∞−→ 0, (22)
then every weak limit of (ν(n)(ω,dα))n≥1 is the distribution of a stationary process solu-
tion to (15).
(b) Assume that the invariant distribution is unique. Let ε≤ 0 such that (C3,ε) holds.
Then, a.s., for every T > 0, for every Sk-continuous functional F :D(R+,R
d)→R, (18)
holds for FT if there exist ρ ∈ [0,1) and s≥ 2, such that
|FT (α)| ≤C sup
0≤t≤T
V (ρ(p+a−1))/s(αt) ∀α ∈D(R+,Rd)
and if
∆N(k,T )
s(1−ε)
is non-increasing and
∆N(k,T )
s(1−ε)
<+∞. (23)
Remark 8. In (22), both assumptions imply the invariance of every weak limit of
0 (ω,dx)). These two assumptions are very different. The first is needed in Proposition
2 for using the Echeverria–Weiss invariance criteria (see Ethier and Kurtz [7], page 238,
Lamberton and Pagès [12] and Lemaire [14]), whereas the second appears in Theorem
2, where our functional approach shows that under some mild additional conditions on
steps and weights, every weak limit is always invariant.
For (23), we refer to Remark 6 for simple sufficient conditions when (γn) and (ηn) are
some polynomial steps and weights.
158 G. Pagès and F. Panloup
4. Proofs of Theorems 1 and 2
We begin the proof with some technical lemmas. In Lemma 1, we show that the a.s
weak convergence of the random measures (ν(n)(ω,dα))n≥1 can be characterized by the
convergence (11) along the set of bounded Lipschitz functionals F for the distance d.
Then, in Lemma 2, we show with some martingale arguments that if the functional
F depends only on the restriction of the trajectory to [0, T ], then the convergence of
(ν(n)(ω,F ))n≥1 is equivalent to that of a more regular sequence. This step is fundamental
for the sequel of the proof.
Finally, Lemma 4 is needed for the proof of Theorem 2. We show that under some mild
conditions on the step and weight sequences, any Markovian weak limit of the sequence
(ν(n)(ω,dα))n≥1 is stationary.
4.1. Preliminary lemmas
Lemma 1. Let (E,d) be a Polish space and let P(E) denote the set of probability
measures on the Borel σ-field B(E), endowed with the weak convergence topology. Let
(µ(n)(ω,dα))n≥1 be a sequence of random probabilities defined on Ω×B(E).
(a) Assume that there exists µ(∞) ∈ P(E) such that for every bounded Lipschitz func-
tion F :E→R,
µ(n)(ω,F )
n→+∞−→ µ(∞)(F ) a.s. (24)
Then, a.s., (µ(n)(ω,dα))n≥1 converges weakly to µ
(∞) on P(E).
(b) Let U be a subset of P(E). Assume that for every sequence (Fk)k≥1 of Lipschitz
and bounded functions, a.s., for every subsequence (µ(φω(n))(ω,dα)), there exists a sub-
sequence (µ(φω◦ψω(n))(ω,dα)) and a U -valued random probability µ(∞)(ω,dα) such that
for every k ≥ 1,
µ(ψω◦φω(n))(ω,Fk)
n→+∞−→ µ(∞)(ω,Fk) a.s. (25)
Then, (µ(n)(ω,dα))n≥1 is a.s. tight with weak limits in U .
Proof. We do not give a detailed proof of the next lemma, which is essentially based
on the fact that in a separable metric space (E,d), one can build a sequence of bounded
Lipschitz functions (gk)k≥1 such that for any sequence (µn)n≥1 of probability measures
on B(E), (µn)n≥1 weakly converges to a probability µ if and only if the convergence
holds along the functions gk, k ≥ 1 (see Parthasarathy [22], Theorem 6.6, page 47 for a
very similar result). �
For every n≥ 0, for every T > 0, we introduce τ(n,T ) defined by
τ(n,T ) := min{k ≥ 0,N(k,T )≥ n}=min{k ≤ n,Γk + T ≥ Γn}. (26)
Approximation of the distribution of a stationary Markov process 159
Note that for k ∈ {0, . . . , τ(n,T )− 1}, {X̄(k)t ,0≤ t≤ T } is �FΓn -measurable and
T − γτ(n,T )−1 ≤ Γn − Γτ(n,T ) ≤ T.
Lemma 2. Assume (C3,ε) with ε < 1. Let F :D(R+,R
d)→R be a Sk-continuous func-
tional. Let (Gk) be a filtration such that F̄Γk ⊂ Gk for every k ≥ 1. Then, for any T > 0:
(a) if FT (defined by (7)) is bounded,
ηk(FT (X̄
(k−1))−E[FT (X̄(k−1))/Gk−1])
n→+∞−→ 0 a.s.; (27)
(b) if FT is not bounded, (27) holds if there exists V :Rd→R+, satisfying H(s, ε) for
some s≥ 2, such that |FT (α)| ≤C sup0≤t≤T V(αt) for every α ∈D(R+,Rd); furthermore,
ν(n)(ω,FT )<+∞ a.s. (28)
Proof. We prove (a) and (b) simultaneously. Let Υ(k) be defined by Υ(k) = FT (X̄
(k)).
We have
(k−1) −E[Υ(k−1)/Gk−1])
(k−1) −E[Υ(k−1)/Gn]) (29)
ηk(E[Υ
(k−1)/Gn]−E[Υ(k−1)/Gk−1]). (30)
We have to prove that the right-hand side of (29) and (30) tend to 0 a.s. when n→+∞.
We first focus on the right-hand side of (29). From the very definition of τ(n,T ), we
have that {X̄(k)t ,0≤ t≤ T } is F̄Γn -measurable for k ∈ {0, . . . , τ(n,T )− 1}. Hence, since
FT is σ(πs,0≤ s≤ T )-measurable and F̄Γn ⊂ Gn, it follows that Υ(k) is Gn-measurable
and that Υ(k) = E[Υ(k)/Gn] for every k ≤ τ(n,T )− 1. Then, if FT is bounded, we derive
from (C3,ε) that
(k−1) −E[Υ(k−1)/Gn])
≤ 2‖FT ‖sup
k=τ(n,T )+1
k=τ(n,T )+1
H1−εn
(Γn − Γτ(n,T ))
160 G. Pagès and F. Panloup
≤ C(T )
H1−εn
n→+∞−→ 0 a.s.,
where we used the fact that (Hn)n≥1 and (γn)n≥1 are non-decreasing and non-increasing
sequences, respectively.
Assume, now, that the assumptions of (b) are fulfilled with V satisfying H(s, ε) for
some s≥ 2 and ε < 1. By the Borel–Cantelli-like argument, it suffices to show that
k=τ(n,T )+1
(k−1) −E[Υ(k−1)/Gn])
<+∞. (31)
Let us prove (31). Let ak := η
(s−1)/s
k and bk(ω) := η
(k−1) − E[Υ(k−1)/Gn]). The
Hölder inequality applied with p̄= s/(s− 1) and q̄ = s yields
k=τ(n,T )+1
akbk(ω)
k=τ(n,T )+1
)s−1( n
k=τ(n,T )+1
ηk|Υ(k−1) −E[Υ(k−1)/Gn]|s
Now, since FT (α) ≤ sup0≤t≤T V(α), it follows from the Markov property and from
H(s, ε)(i) that
E[|FT (X̄(k))|s/F̄Γk ]≤CE
0≤t≤T
Vs(X̄(k)t )/F̄Γk
≤CTVs(X̄Γk).
Then, using the two preceding inequalities and (C3,ε) yields
k=τ(n,T )+1
(k−1) −E[Υ(k−1)/Gn])
k=τ(n,T )+1
)s−1( n
k=τ(n,T )+1
ηkE[Vs(X̄Γk−1)]
k=τ(n,T )+1
k=τ(n,T )+1
Vs(X̄Γk−1)
k=τ(n,T )+1
t∈[0,S(n,T )]
Vs(X̄τ(n,T )t )
where S(n,T ) = Γn−1 − Γτ(n,T ) and C does not depend n. By the definition of τ(n,T ),
S(n,T )≤ T . Then, again using H(s, ε)(i) yields
k=τ(n,T )
(k−1) −E[Υ(k−1)/Gn])
s(1−ε)
E[Vs(X̄(τ(n,T )))].
Approximation of the distribution of a stationary Markov process 161
Since n 7→ N(n,T ) is an increasing function, n 7→ τ(n,T ) is a non-decreasing function
and Card{n, τ(n,T ) = k}=∆N(k+1, T ) :=N(k+1, T )−N(k,T ). Then, since n 7→Hn
increases, a change of variable yields
k=τ(n,T )+1
(k−1) −E[Υ(k−1)/Gn])
∆N(k,T )
s(1−ε)
E[Vs(X̄Γk−1)]<+∞,
by H(s, ε)(iv).
Second, we prove that (30) tends to 0. For every n≥ 1, we let
(E[Υ(k−1)/Gn]−E[Υ(k−1)/Gk−1]). (32)
The process (Mn)n≥1 is a (Gn)-martingale and we want to prove that this process is
L2-bounded. Set Φ(k,n) = E[FT (X̄
(k))/Gn]− E[FT (X̄(k))/Gk]. Since FT is σ(πs,0 ≤ s ≤
T )-measurable, the random variable Φ(k,n) is F̄ΓN(k,T) -measurable. Then, for every i ∈
{N(k,T ), . . . , n}, Φ(k,n) is Gi-measurable so that
E[Φ(i,n)Φ(k,n)] =E[Φ(k,n)E[Φ(i,n)/Gi]] = 0.
It follows that
E[M2n] =
E[(Φ(k−1,n))
] + 2
N(k−1,T )∧n
i=k+1
E[Φ(i−1,n)Φ(k−1,n)]. (33)
Then,
E[M2n] ≤
E[(Φ(k−1,n))
] + 2
N(k−1,T )
i=k+1
E[Φ(i−1,n)Φ(k−1,n)]
H2−εk
E[(Φ(k−1,n))
] (34)
H2−εk
N(k−1,T )
i=k+1
γi sup
E[Φ(i−1,n)Φ(k−1,n)]
162 G. Pagès and F. Panloup
where, in the second inequality, we used assumption (C3,ε) and the decrease of i 7→
1/H1−εi . Hence, if FT is bounded, using the fact that
∑N(k−1,T )
i=k+1 γi ≤ T yields
E[M2n]≤C
H2−εk
H2−ε1
<+∞ (35)
since ε < 1. Assume, now, that the assumptions of (b) hold and let FT be dominated
by a function V satisfying H(s, ε). By the Markov property, the Jensen inequality and
H(s, ε)(i),
E[(Φ(k,n))
0≤t≤T
V2(X̄(k)t )/F̄Γk
≤CTE[V2(X̄Γk)].
We then derive from the Cauchy–Schwarz inequality that for every n, k ≥ 1, for every
i ∈ {k, . . . ,N(k,T )},
|E[Φ(i,n)Φ(k,n)]| ≤C
E[V2(X̄Γi)]
E[V2(X̄Γk)]≤C sup
t∈[0,T ]
E[V2(X̄(k)t )]≤CE[V2(X̄Γk)],
where, in the last inequality, we once again used H(s, ε)(i). It follows that
E[M2n]≤C
H2−εk
E[V2(X̄Γk−1)]<+∞,
by H(s, ε)(iii). Therefore, (34) is finite and (Mn) is bounded in L
2. Finally, we derive
from the Kronecker lemma that
ηk(E[FT (X̄
(k−1))/Gn]−E[FT (X̄(k−1))/Gk−1])
n→+∞−→ 0 a.s.
As a consequence, supn≥1 ν
(n)(ω,FT )<+∞ a.s. if and only if
E[FT (X̄
(k−1))/Fk−1]<+∞ a.s.
This last property is easily derived from H(s, ε)(i) and (ii). This completes the proof. �
Lemma 3. (a) Assume (C1) and let x0 ∈Rd. We then have limx→x0 E[d(Xx,Xx0)] = 0.
In particular, for every bounded Lispchitz (w.r.t. the distance d) functional F :D(R+,R
R, the function ΦF defined by ΦF (x) = E[F (Xx)] is a (bounded) continuous function on
(b) Assume (C2). For every compact set K ⊂Rd,
E[d(Y n,x,Xx)]
n→+∞−→ 0. (36)
Approximation of the distribution of a stationary Markov process 163
Set ΦFn (x) = E[F (Y
n,x)]. Then, for every bounded Lispchitz functional F :D(R+,R
d)→R,
|ΦF (x)−ΦFn (x)|
n→+∞−→ 0 for every compact set K ⊂Rd. (37)
Proof. (a) By the definition of d, for every α, β ∈D(R+,Rd) and for every T > 0,
d(α,β)≤
1∧ sup
0≤t≤T
|α(t)− β(t)|
+ e−T . (38)
It easily follows from assumption (C1) and from the dominated convergence theorem
limsup
E[d(Xx,Xx0)]≤ e−T for every T > 0.
Letting T →+∞ implies that limx→x0 E[d(Xx,Xx0)] = 0.
(b) We deduce from (38) and from assumption (C2) that for every compact setK ⊂Rd,
for every T > 0,
limsup
E[d(Y n,x,Xx)]≤ e−T .
Letting T →+∞ yields (36). �
Lemma 4. Assume that (ηn)n≥1 and (γn) satisfy (C3,ε) with ε < 1 and (14). Then:
(i) for every t≥ 0, for every bounded continuous function f :Rd→R,
t (ω, f)− ν
0 (ω, f)
n→+∞−→ 0 a.s.;
(ii) if, moreover, a.s., every weak limit ν(∞)(ω,dα) of (ν(n)(ω,dα))n≥1 is the dis-
tribution of a Markov process with semigroup (Qωt )t≥0, then, a.s., ν
(∞)(ω,dα) is the
distribution of a stationary process.
Proof. (i) Let f :Rd →R be a bounded continuous function. Since X̄(k)t = X̄ΓN(k,t) , we
t (ω, f)− ν
0 (ω, f) =
ηk(f(X̄ΓN(k−1,t))− f(X̄Γk−1)).
From the very definition of N(n,T ) and τ(n,T ), one checks that N(k − 1, T )≤ n− 1 if
and only if τ(n,T )≥ k. Then,
ηkf(X̄Γk−1) =
τ(n,t)
ηN(k−1,t)+1f(X̄ΓN(k−1,t))
ηkf(X̄Γk−1)1{k−1/∈N({0,...,n},t)}.
164 G. Pagès and F. Panloup
It follows that
t (ω, f)− ν
0 (ω, f) =
τ(n,t)
(ηk − ηN(k−1,t)+1)f(X̄ΓN(k−1,t))
τ(n,t)+1
ηkf(X̄ΓN(k−1,t))
ηkf(X̄Γk−1)1{k−1/∈N({0,...,n},t)}.
Then, since f is bounded and since
ηk1{k−1/∈N({0,...,n},t)} =
τ(n,t)
ηN(k−1,t)+1
τ(n,t)
|ηk − ηN(k−1,t)+1|+
k=τ(n,t)+1
we deduce that
|ν(n)t (ω, f)− ν
0 (ω, f)| ≤ 2‖f‖∞
τ(n,t)
|ηk − ηN(k−1,t)+1|+
k=τ(n,t)+1
Hence, we have to show that the sequences of the right-hand side of the preceding in-
equality tend to 0. On the one hand, we observe that
|ηk − ηN(k−1,t)+1| ≤
N(k−1,T )+1
ℓ=k+1
|ηℓ − ηℓ−1| ≤ max
ℓ≥k+1
|∆ηℓ|
N(k−1,T )+1
Using the fact that
∑N(k−1,T )+1
ℓ=k γℓ ≤ T + γ1 and condition (14) yields
τ(n,t)
|ηk − ηN(k−1,t)+1|
n→+∞−→ 0.
On the other hand, by (C3,ε), we have
k=τ(n,T )+1
H1−εn
k=τ(n,T )+1
H1−εn
n→+∞−→ 0 a.s.,
which completes the proof of (i).
Approximation of the distribution of a stationary Markov process 165
(ii) Let Q+ denote the set of non-negative rational numbers. Let (fℓ)ℓ≥1 be an every-
where dense sequence in CK(Rd) endowed with the topology of uniform convergence on
compact sets. Since Q+ and (fℓ)ℓ≥1 are countable, we derive from (i) that there exists
Ω̃⊂Ω such that P(Ω̃) = 1 and such that for every ω ∈ Ω̃, every t ∈Q+ and every ℓ≥ 1,
t (ω, fℓ)− ν
0 (ω, fℓ)
n→+∞−→ 0.
Let ω ∈ Ω̃ and let ν(∞)(ω,dα) denote a weak limit of (ν(n)(ω,dα))n≥1. We have
t (ω, fℓ) = ν
0 (ω, fℓ) ∀t ∈Q+ ∀ℓ≥ 1
and we easily deduce that
t (ω, f) = ν
0 (ω, f) ∀t ∈R+ ∀f ∈ CK(Rd).
Hence, if ν(∞)(ω,dα) is the distribution of a Markov process (Yt) with semigroup (Q
t )t≥0,
we have, for all f ∈ CK(Rd),
Qωt f(x)ν
0 (ω,dx) =
f(x)ν
0 (ω,dx) ∀t≥ 0.
0 (ω,dx) is then an invariant distribution for (Yt). This completes the proof. �
4.2. Proof of Theorem 1
Thanks to Lemma 1(a) applied with E =D(R+,R
d) and d defined by (6),
ν(n)(ω,dα)
=⇒ Pν0(dα) a.s.⇐⇒ ν(n)(ω,F )
n→+∞−→
F (x)Pν0 (dx) a.s. (39)
for every bounded Lipschitz functional F :D(R+,R
d)→ R. Now, consider such a func-
tional. By the assumptions of Theorem 1, we know that a.s., (ν
0 (ω,dx))n≥1 converges
weakly to ν0. Set Φ
F (x) := E[F (Xx)], x ∈Rd. By Lemma 3(a), ΦF is a bounded contin-
uous function on Rd. It then follows from (C0,1) that
F (X̄
(k−1)
n→+∞−→
ΦF (x)ν0(dx) =
F (x)Pν0 (dx) a.s.
Hence, the right-hand side of (39) holds for F as soon as
ηk(F (X̄
(k−1))−ΦF (X̄(k−1)0 ))
n→+∞−→ 0 a.s. (40)
166 G. Pagès and F. Panloup
Let us prove (40). First, let T > 0 and let FT be defined by (7). By Lemma 2,
ηkFT (X̄
(k−1))− 1
ηkE[FT (X̄
(k−1))/F̄Γk−1 ]
n→+∞−→ 0 a.s. (41)
With the notation of Lemma 3(b), we derive from assumption (C2)(i) that
E[FT (X̄
(k−1))/F̄Γk−1 ] = Φ
k (X̄
(k−1)
Let N ∈N. On one hand, by Lemma 3(b),
k (X̄
(k−1)
0 )−ΦFT (X̄
(k−1)
0 ))1{|X̄(k−1)
n→+∞−→ 0 a.s. (42)
On the other hand, the tightness of (ν
0 (ω,dx))n≥1 on R
d yields
ψ(ω,N) := sup
0 (ω, (B(0,N)
N→+∞−→ 0 a.s.
It follows that, a.s.,
ηk|ΦFTk (X̄
(k−1)
0 )−ΦFT (X̄
(k−1)
0 )|1{|X̄(k−1)
≤ 2‖F‖∞ψ(ω,N)
N→+∞−→ 0.
Hence, a combination of (42) and (43) yields
∀T > 0 1
k (X̄
(k−1)
0 )−ΦFT (X̄
(k−1)
n→+∞−→ 0 a.s. (44)
Finally, let (Tℓ)ℓ≥1 be a sequence of positive numbers such that, Tℓ→+∞ when ℓ→+∞.
Combining (44) and (41), we obtain that, a.s., for every ℓ≥ 1,
limsup
ηk(F (X̄
(k−1))−ΦF (X̄(k−1)))
≤ lim sup
ηk(F (X̄
(k−1))−FTℓ(X̄(k−1)))
+ limsup
FTℓ (X̄
(k−1)
0 )−ΦF (X̄
(k−1)
Approximation of the distribution of a stationary Markov process 167
By the definition of d, |F − FTℓ | ≤ e−Tℓ . Then, a.s.,
limsup
ηk(F (X̄
(k−1))−ΦF (X̄(k−1)0 ))
≤ 2e−Tℓ ∀ℓ≥ 1.
Letting ℓ→+∞ implies (40).
The generalization to non-bounded functionals in Theorem 1 is then derived from (28)
and from a uniform integrability argument.
4.3. Proof of Theorem 2
(i) We want to prove that the conditions of Lemma 1(b) are fulfilled. Since (ν
0 (ω,dx))n≥1
is supposed to be a.s. tight, one can check that for every bounded Lipschitz functional
F :D(R+,R
d)→R, (40) is still valid. Then, let (Fℓ)ℓ≥1 be a sequence of bounded Lipschitz
functionals. There exists Ω̃⊂Ω with P(Ω̃) = 1 such that for every ω ∈ Ω̃, (ν(n)0 (ω,dx))n≥1
is tight and
ηk(Fℓ(X̄
(k−1)(ω))−ΦFℓ(X̄(k−1)0 (ω)))
n→+∞−→ 0 ∀ℓ≥ 1. (45)
Let ω ∈ Ω̃ and let φω :N 7→N be an increasing function. As (ν(φω(n))0 (ω,dx))n≥1 is tight,
there exists a convergent subsequence (ν
(φω◦ψω(n))
0 (ω,dx))n≥1. We denote its weak limit
by ν∞. Since Φ
Fℓ is continuous for every ℓ≥ 1 (see Lemma 3(a)),
(φω◦ψω(n))
0 (ω,Φ
n→+∞−→ ν∞(ΦFℓ) =
Fℓ(α)Pν∞(dα) ∀ℓ≥ 1.
We then derive from (45) that for every ℓ≥ 1
ν(φω◦ψω(n))(ω,Fℓ)
n→+∞−→
Fℓ(α)Pν∞(dα).
It follows that the conditions of Lemma 1(b) are fulfilled with U = {Pµ, µ ∈ I}, where
µ ∈P(Rd),∃ω ∈ Ω̃ and an increasing function φ :N 7→N, µ= lim
ν(φ(n))(ω,dα)
Hence, by Lemma 1(b), we deduce that (ν(n)(ω,dα))n≥1 is a.s. tight with U -valued limits.
Finally, Theorem 2(ii) is a consequence of condition (14) and Lemma 4(ii).
168 G. Pagès and F. Panloup
5. Path-dependent option pricing in stationary
stochastic volatility models
In this section, we propose a simple and efficient method to price options in stationary
stochastic volatility (SSV) models. In most stochastic volatility (SV) models, the volatil-
ity is a mean reverting process. These processes are generally ergodic with a unique
invariant distribution (the Heston model or the BNS model for instance (see below) but
also the SABR model (see Hagan et al. [8]), . . .). However, they are usually considered
in SV models under a non-stationary regime, starting from a deterministic value (which
usually turns out to be the mean of their invariant distribution). However, the instanta-
neous volatility is not easy to observe on the market since it is not a traded asset. Hence,
it seems to be more natural to assume that it evolves under its stationary regime than
to give it a deterministic value at time 0.4
From a purely calibration viewpoint, considering an SV model in its SSV regime will
not modify the set of parameters used to generate the implied volatility surface, although
it will modify its shape, mainly for short maturities. This effect can in fact be an asset
of the SSV approach since it may correct some observed drawbacks of some models (see,
e.g., the Heston model below).
From a numerical point of view, considering SSV models is no longer an obstacle, es-
pecially when considering multi-asset models (in the unidimensional case, the stationary
distribution can be made more or less explicit like in the Heston model; see below) since
our algorithm is precisely devised to compute by simulation some expectations of func-
tionals of processes under their stationary regime, even if this stationary regime cannot
be directly simulated.
As a first illustration (and a benchmark) of the method, we will describe in detail
the algorithm for the pricing of Asian options in a Heston model. We will then show
in our numerical results to what extent it differs, in terms of smile and skew, from the
usual SV Heston model for short maturities. Finally, we will complete this section with
a numerical test on Asian options in the BNS model where the volatility is driven by a
tempered stable subordinator. Let us also mention that this method can be applied to
other fields of finance like interest rates, and commodities and energy derivatives where
mean-reverting processes play an important role.
4When one has sufficiently close observations of the stock price, it is in fact possible to derive a rough
idea of the size of the volatility from the variations of the stock price (see, e.g., Jacod [10]). Then, using
this information, a good compromise between a deterministic initial value and the stationary case may
be to assume that the distribution µ0 of the volatility at time 0 is concentrated around the estimated
value (see Section 2.2 for application of our algorithm in this case).
Approximation of the distribution of a stationary Markov process 169
5.1. Option pricing in the Heston SSV model
We consider a Heston stochastic volatility model. The dynamic of the asset price process
(St)t≥0 is given by S0 = s0 and
dSt = St(rdt+
(1− ρ2)vt dW 1t + ρ
vt dW
dvt = k(θ− vt) dt+ ς
vt dW
where r denotes the interest rate, (W 1,W 2) is a standard two-dimensional Brownian
motion, ρ ∈ [−1,1] and k, θ and ς are some non-negative numbers. This model was
introduced by Heston in 1993 (see Heston [9]). The equation for (vt) has a unique (strong)
pathwise continuous solution living in R+. If, moreover, 2kθ > ς
2, then (vt) is a positive
process (see Lamberton and Lapeyre [11]). In this case, (vt) has a unique invariant
probability ν0. Moreover, ν0 = γ(a, b) with a= (2k)/ς
2 and b= (2kθ)/ς2. In the following,
we will assume that (vt) is in its stationary regime, that is, that
L(v0) = ν0.
5.1.1. Option price and stationary processes
Using our procedure to price options in this model naturally needs to express the option
price as the expectation of a functional of a stationary stochastic process.
Näıve method. (may work) Since (vt)t≥0 is stationary, the first idea is to express the
option price as the expectation of a functional of (vt)t≥0: by Itô calculus, we have
St = s0 exp
rt− 1
vs ds
vs dW
1− ρ2
vs dW
. (46)
Since
vs dW
s =Λ(t, (vt)) :=
vt − v0 − kθt+ k
vs ds
it follows by setting Mt =
vs dW
s that
St =Ψ(t, (vs), (Ms)), (47)
where Ψ is given for every t≥ 0, u and w ∈ C(R+,R) by
Ψ(t, u,w) = s0 exp
rt− 1
u(s) ds
+ ρΛ(t, u) +
1− ρ2w(t)
Then, let F :C(R+,R) → R be a non-negative measurable functional. Conditioning by
FW 2T yields
E[FT ((St)t≥0)] = E[F̃T ((vt)t≥0)],
170 G. Pagès and F. Panloup
where, for every u ∈ C(R+,R),
F̃T (u) = E
t, u,
u(s) dW 1s
For some particular options such as the European call or put (thanks to the Black–
Scholes formula), the functional F̃ is explicit. In those cases, this method seems to be
very efficient (see Panloup [20] for numerical results). However, in the general case, the
computation of F̃ will need some Monte Carlo methods at each step. This approach is
then very time-consuming in general – that is why we are going to introduce another
representation of the option as a functional of a stationary process.
General method. (always works) We express the option premium as the expectation
of a functional of a two-dimensional stationary stochastic process. This method is based
on the following idea. Even though (vt,Mt) is not stationary, (St) can be expressed as a
functional of a stationary process (vt, yt). Indeed, consider the following SDE given by
dyt =−yt dt+
vt dW
dvt = k(θ− vt) dt+ ς
vt dW
First, one checks that the SDE has a unique strong solution and that assumption (S1) is
fulfilled with V (x1, x2) = 1+ x
2. This ensures the existence of an invariant distribu-
tion ν̃0 for the SDE (see, e.g., Pagès [17]). Then, since (vt) is positive and has a unique
invariant distribution, the uniqueness of the invariant distribution follows. Then, assume
that L(y0, v0) = ν̃0. Since (vt,Mt) = (vt, yt − y0 +
ys ds), we have, for every positive
measurable functional F :C(R+,R)→R,
E[FT ((St)t≥0)] = E[FT ((ψ(t, vt,Mt))t≥0)]
= Eν̃0
t, vt, yt − y0 +
ys ds
where Pν̃0 is the stationary distribution of the process (vt, yt). Every option price can
then be expressed as the expectation of an explicit functional of a stationary process. We
will develop this second general approach in the numerical tests below.
Remark 9. The idea of the second method holds for every stochastic volatility model
for which (St) can be written as follows:
St =Φ
t, vt,
hi(|vs|) dY is
, (50)
where, for every i ∈ {1, . . . , p}, hi :R+ →R is a positive function such that hi(x) = o(|x|)
as |x| → +∞, (Y it ) is a square-integrable centered Lévy process and (vt) is a mean
reverting stochastic process solution to a Lévy driven SDE.
Approximation of the distribution of a stationary Markov process 171
In some complex models, showing the uniqueness of the invariant distribution may be
difficult. In fact, it is important to note at this stage that the uniqueness of the invariant
distribution for the couple (vt, yt) is not required. Indeed, by construction, the local
martingale (Mt) does not depend on the choice of y0. It follows that if L(y0, v0) = µ̃,
with µ̃ constructed such that L(v0) = ν0, (49) still holds. This implies that it is only
necessary that uniqueness holds for the invariant distribution of the stochastic volatility
process.
5.1.2. Numerical tests on Asian options
We recall that (vt) is a Cox–Ingersoll–Ross process. For this type of processes, it is well
known that the genuine Euler scheme cannot be implemented since it does not preserve
the non-negativity of the (vt). That is why some specific discretization schemes have
been studied by several authors (Alfonsi [1], Deelstra and Delbaen [5] and Berkaoui et al.
[4, 6]). In this paper, we consider the scheme studied by the last authors in a decreasing
step framework. We denote it by (v̄t). We set v̄0 = x > 0 and
v̄Γn+1 = |v̄Γn + kγn+1(θ− v̄Γn) + ς
v̄Γn(W
−W 2Γn)|.
We also introduce the stepwise constant Euler scheme (ȳt) of (yt)t≥0 defined by
ȳΓn+1 = ȳΓn − γn+1ȳΓn +
v̄Γn(W̃
− W̃ 1Γn), ȳ0 = y ∈R
Denote by (v̄
t ) and (ȳ
t ) the shifted processes defined by v̄
t := v̄Γk+t and ȳ
ȳΓk+t, and let (ν
(n)(ω,dα))n≥1 be the sequence of empirical measures defined by
ν(n)(ω,dα) =
ηk1{(v̄(k−1),ȳ(k−1))∈dα}.
The specificity of both the model and the Euler scheme implies that Theorems 1 and 2
cannot be directly applied here. However, a specific study using the fact that (9) holds
for every compact set of R∗+ ×R when 2kθ/ς2 > 1+ 2
6/ς (see Theorem 2.2 of Berkaoui
et al. [4] and Remark 9) shows that
ν(n)(ω,dα)
=⇒ Pν̃0(dα) a.s.
when 2kθ/ς2 > 1+ 2
6/ς . Details are left to the reader.
Let us now state our numerical results obtained for the pricing of Asian options with
this discretization. We denote by Cas(ν0,K,T ) and Pas(ν0,K,T ) the Asian call and put
prices in the SSV Heston model. We have
Cas(ν0,K,T ) = e
Ss ds−K
172 G. Pagès and F. Panloup
Pas(ν0,K,T ) = e
K − 1
Ss ds
With the notation of (49), approximating Cas(ν0,K,T ) and Pas(ν0,K,T ) by our proce-
dure needs to simulate the sequences (Cnas)n≥1 and (P
as)n≥1 defined by
Cnas =
Ψ(s, v̄(k−1), M̄ (k−1)) ds−K
Pnas =
K − 1
Ψ(s, v̄(k−1), M̄ (k−1)) ds
These sequences can be computed by the method developed in Section 1.3. Note that
the specific properties of the exponential function and the linearity of the integral imply
that (
Ψ(t, v̄(n−1), M̄ (n−1)) ds) can be computed quasi-recursively.
Let us state our numerical results for the Asian call with parameters
s0 = 50, r = 0.05, T = 1, ρ= 0.5,
θ = 0.01, ς = 0.1, k = 2.
We also assume that K ∈ {44, . . . ,56} and choose the following steps and weights: γn =
ηn = n
−1/3. In Table 1, we first state the reference value for the Asian call price obtained
for N = 108 iterations. In the two following lines, we state our results for N = 5.104 and
N = 5.105 iterations. Then, in the last lines, we present the numerical results obtained
Table 1. Approximation of the Asian call price
K 44 45 46 47 48 49 50
Asian call (ref.) 6.92 5.97 5.04 4.12 3.25 2.46 1.78
N = 5 · 104 6.89 6.07 5.07 4.13 3.18 2.49 1.77
N = 5 · 105 6.90 6.02 5.00 4.11 3.24 2.46 1.79
N = 5 · 104 (CP parity) 6.92 5.96 5.04 4.13 3.26 2.46 1.78
N = 5 · 105 (CP parity) 6.92 5.97 5.04 4.12 3.25 2.47 1.78
K 51 52 53 54 55 56
Asian call (ref.) 1.23 0.82 0.53 0.33 0.21 0.12
N = 5 · 104 1.21 0.81 0.51 0.34 0.22 0.11
N = 5 · 105 1.23 0.82 0.53 0.33 0.21 0.13
N = 5 · 104 (CP parity) 1.23 0.82 0.53 0.31 0.21 0.12
N = 5 · 105(CP parity) 1.23 0.82 0.53 0.33 0.21 0.13
Approximation of the distribution of a stationary Markov process 173
using the call-put parity
Cas(ν0,K,T )− Pas(ν0, S0,K,T ) =
(1− e−rT )−Ke−rT (52)
as a means of variance reduction. The computation times for N = 5.104 and N = 5.105
(using MATLAB with a Xeon 2.4 GHz processor) are about 5 s and 51 s, respectively. In
particular, the complexity is quasi-linear and the additional computations needed when
we use the call-put parity are negligible.
5.2. Implied volatility surfaces of Heston SSV and SV models
Given a particular pricing model (with initial value s0 and interest rate r) and its asso-
ciated European call prices denoted by Ceur(K,T ), we recall that the implied volatility
surface is the graph of the function (K,T ) 7→ σimp(K,T ), where σimp(K,T ) is defined for
every maturity T > 0 and strike K as the unique solution of
CBS(s0,K,T, r, σimp(K,T )) =Ceur(K,T ),
where CBS(s0,K,T, r, σ) is the price of the European call in the Black–Scholes model
with parameters s0, r and σ. When Ceur(K,T ) is known, the value of σimp(K,T ) can be
numerically computed using the Newton method or by dichotomy if the first method is
not convergent.
In this last part, we compare the implied volatility surfaces induced by the SSV and SV
Heston models where we suppose that the initial value of (vt) in the SV Heston model is
the mean of the invariant distribution, that is, we suppose that v0 = θ.
5 We also assume
that the parameters are those of (51), except the correlation coefficient ρ.
In Figures 1 and 2, the volatility curves obtained when T = 1 are depicted, whereas in
Figures 3 and 4, we set the strikeK atK = 50 and let the time vary. These representations
show that when the maturity is long, the differences between the SSV and SV Heston
models vanish. This is a consequence of the convergence of the stochastic volatility to its
stationary regime when T →+∞.
The main differences between these models then appear for short maturities. That is
why we complete this part by a representation of the volatility curve when T = 0.1 for
ρ= 0 and ρ= 0.5 in Figures 5 and 6, respectively. We observe that for short maturities,
the volatility smile is more curved and the skew is steeper. These phenomena seem
interesting for calibration since one well-known drawback of the standard Heston model
is that it can have overly flat volatility curves for short maturities.
5.3. Numerical tests on Asian options in the BNS SSV model
The BNS model introduced in Barndorff-Nielsen and Shephard [3] is a stochastic volatility
model where the volatility process is a Lévy-driven positive Ornstein–Uhlenbeck process.
5This choice is the most usual in practice.
174 G. Pagès and F. Panloup
Figure 1. ρ= 0, K 7→ σimp(K,1).
The dynamic of the asset price (St) is given by St = S0 exp(Xt),
dXt = (r− 12vt) dt+
vt dWt + ρdZt, ρ≤ 0,
dvt = −µvt dt+dZt, µ > 0,
Figure 2. ρ= 0.5, K 7→ σimp(K,1).
Approximation of the distribution of a stationary Markov process 175
Figure 3. ρ= 0, T 7→ σimp(50, T ).
where (Zt) is a subordinator without drift term and Lévy measure π. In the following,
we assume that (Zt) is a tempered stable subordinator, that is, that
π(dy) = 1{y>0}
c exp(−λy)
dy, c > 0, λ > 0, α∈ (0,1).
As in the Heston model, we want to use our algorithm as a way of option pricing when
the stochastic volatility evolves under its stationary regime and test it on Asian options
using the method described in detail in Section 5.1. This model does not require a specific
Figure 4. ρ= 0.5, T 7→ σimp(50, T ).
176 G. Pagès and F. Panloup
Figure 5. ρ= 0, T 7→ σimp(50, T ).
discretization and the approximate Euler scheme (P) (see Section 3.2) relative to (vt)
can be implemented using the rejection method. In Table 2, we present our numerical
results obtained for the following choices of parameters, steps and weights:
ρ=−1, λ= µ= 1, c= 0.01, α= 1
, γn = ηn = n
−1/3.
The computation times forN = 5.104 andN = 5.105 are about 8.5 s and 93 s, respectively.
Note that for this model, the convergence seems to be slower because of the approximation
of the jump component.
Figure 6. ρ= 0.5, T 7→ σimp(50, T ).
Approximation of the distribution of a stationary Markov process 177
Table 2. Approximation of the Asian call price in the BNS model
K 44 45 46 47 48 49 50
Asian call (ref.) 6.75 5.83 4.93 4.05 3.18 2.35 1.57
N = 5 · 104 6.83 5.91 5.01 4.10 3.22 2.35 1.51
N = 5 · 105 6.78 5.86 4.96 4.06 3.19 2.34 1.52
N = 5 · 104 (CP parity) 6.76 5.85 4.94 4.07 3.20 2.29 1.51
N = 5 · 105 (CP parity) 6.75 5.83 4.93 4.04 3.17 2.32 1.54
K 51 52 53 54 55 56
Asian call (ref.) 0.91 0.55 0.39 0.29 0.23 0.18
N = 5 · 104 0.77 0.46 0.33 0.27 0.22 0.19
N = 5 · 105 0.79 0.48 0.34 0.27 0.21 0.17
N = 5 · 104 (CP parity) 0.79 0.47 0.37 0.27 0.23 0.19
N = 5 · 105(CP parity) 0.83 0.50 0.36 0.28 0.22 0.17
Acknowledgement
The authors would like to thank Vlad Bally for interesting comments on the paper.
References
[1] Alfonsi, A. (2005). On the discretization schemes for the CIR (and Bessel squared) pro-
cesses. Monte Carlo Methods Appl. 11 355–384. MR2186814
[2] Asmussen, S. and Rosinski, J. (2001). Approximations of small jumps of Lévy processes
with a view towards simulation. J. Appl. Probab. 38 482–493. MR1834755
[3] Barndorff-Nielsen, O.E. and Shephard, N. (2001). Modelling by Lévy processes for financial
economics. In Lévy Processes 283–318. Boston: Birkhäuser. MR1833702
[4] Berkaoui, A., Bossy, M. and Diop, A. (2008). Euler scheme for SDE’s with non-Lipschitz
diffusion coefficient: Strong convergence. ESAIM Probab. Statist. 12 1–11. MR2367990
[5] Deelstra, G. and Delbaen, F. (1998). Convergence of discretized stochastic (interest rate)
processes with stochastic drift term. Appl. Stochastic Models Data Anal. 14 77–84.
MR1641781
[6] Diop, A. (2003). Sur la discrétisation et le comportement à petit bruit d’EDS unidimension-
nelles dont les coefficients sont à dérivées singulières. Ph.D. thesis, Univ. Nice Sophia
Antipolis.
[7] Ethier, S. and Kurtz, T. (1986). Markov Processes, Characterization and Convergence.
Wiley Series in Probability and Mathematical Statistics: Probability and Mathematical
Statistics. New York: Wiley. MR0838085
[8] Hagan, D., Kumar, D., Lesniewsky, A. and Woodward, D. (2002). Managing smile risk.
Wilmott Magazine 9 84–108.
[9] Heston, S. (1993). A closed-form solution for options with stochastic volatility with appli-
cations to bond and currency options. Review of Financial Studies 6 327–343.
[10] Jacod, J. (2008). Asymptotic properties of realized power variations and related functionals
of semimartingales. Stochastic Process. Appl. 118 517–559. MR2394762
http://www.ams.org/mathscinet-getitem?mr=2186814
http://www.ams.org/mathscinet-getitem?mr=1834755
http://www.ams.org/mathscinet-getitem?mr=1833702
http://www.ams.org/mathscinet-getitem?mr=2367990
http://www.ams.org/mathscinet-getitem?mr=1641781
http://www.ams.org/mathscinet-getitem?mr=0838085
http://www.ams.org/mathscinet-getitem?mr=2394762
178 G. Pagès and F. Panloup
[11] Lamberton, D. and Lapeyre, B. (1996). Introduction to Stochastic Calculus Applied to Fi-
nance. London: Chapman and Hall/CRC. MR1422250
[12] Lamberton, D. and Pagès, G. (2002). Recursive computation of the invariant distribution
of a diffusion. Bernoulli 8 367–405. MR1913112
[13] Lamberton, D. and Pagès, G. (2003). Recursive computation of the invariant distribution
of a diffusion: The case of a weakly mean reverting drift. Stoch. Dynamics 4 435–451.
MR2030742
[14] Lemaire, V. (2007). An adaptive scheme for the approximation of dissipative systems.
Stochastic Process. Appl. 117 1491–1518. MR2353037
[15] Lemaire, V. (2005). Estimation numérique de la mesure invariante d’un processus de diffu-
sion. Ph.D. thesis, Univ. Marne-La Vallée.
[16] Pagès, G. (1985). Théorèmes limites pour les semi-martingales. Ph.D. thesis, Univ. Paris
[17] Pagès, G. (2001). Sur quelques algorithmes récursifs pour les probabilités numériques.
ESAIM Probab. Statist. 5 141–170. MR1875668
[18] Panloup, F. (2008). Recursive computation of the invariant measure of a SDE driven by a
Lévy process. Ann. Appl. Probab. 18 379–426. MR2398761
[19] Panloup, F. (2008). Computation of the invariant measure of a Lévy driven SDE: Rate of
convergence. Stochastic Process. Appl. 118 1351–1384.
[20] Panloup, F. (2006). Approximation du régime stationnaire d’une EDS avec sauts. Ph.D.
thesis, Univ. Paris VI.
[21] Protter, P. (1990). Stochastic Integration and Differential Equations. Berlin: Springer.
MR1037262
[22] Parthasarathy, K.R. (1967). Probability Measures on Metric Spaces. New York: Academic
Press. MR0226684
Received April 2007 and revised March 2008
http://www.ams.org/mathscinet-getitem?mr=1422250
http://www.ams.org/mathscinet-getitem?mr=1913112
http://www.ams.org/mathscinet-getitem?mr=2030742
http://www.ams.org/mathscinet-getitem?mr=2353037
http://www.ams.org/mathscinet-getitem?mr=1875668
http://www.ams.org/mathscinet-getitem?mr=2398761
http://www.ams.org/mathscinet-getitem?mr=1037262
http://www.ams.org/mathscinet-getitem?mr=0226684
	Introduction
	Objectives and motivations
	Background and construction of the procedure
	Simulation of ((n)(,F))n1
	General results
	Weak convergence to the stationary regime
	Extension to the non-stationary case
	Application to Brownian diffusions and Lévy-driven SDE's
	Application to Brownian diffusions
	Application to Lévy-driven SDE's
	Proofs of Theorems 1 and 2
	Preliminary lemmas
	Proof of Theorem 1
	Proof of Theorem 2
	Path-dependent option pricing in stationary stochastic volatility models
	Option pricing in the Heston SSV model
	Option price and stationary processes
	Numerical tests on Asian options
	Implied volatility surfaces of Heston SSV and SV models
	Numerical tests on Asian options in the BNS SSV model
	Acknowledgement
	References
