The article is devoted to the estimation of the convergence rate of integral functionals of a Markov process. Under the assumption that the given Markov process admits a transition probability density differentiable in t and the derivative has an integrable upper bound of a certain type, we derive the accuracy rates for strong and weak approximations of the functionals by Riemannian sums. We also develop a version of the parametrix method, which provides the required upper bound for the derivative of the transition probability density for a solution of an SDE driven by a locally stable process. As an application, we give accuracy bounds for an approximation of the price of an occupation time option.
Let Xt, t≥0, be a Markov process with values in Rd. Consider an integral functional of the form
IT(h)=∫0Th(Xt)dt,
where h:Rd→R is a measurable function. In this paper, we investigate the accuracy of the approximation of IT(h) by the Riemannian sums
IT,n(h)=Tn∑k=0n−1h(X(kT)/n),n≥1.
The function h is assumed to be bounded, only, that is, we do not impose any regularity assumptions on h. In particular, under this assumption, the class of integral functionals we investigate contains the class of occupation time type functionals (for which h=1A for a fixed A∈B(Rd)), which are of particular importance.
Integral functionals arise naturally in a wide class of stochastic representation formulae and applied stochastic models. It is very typical that exact calculation of the respective probabilities and/or expectations is hardly possible, which naturally suggests the usage of approximation methods. As an example of such a situation, we mention the so-called occupation time option [11], whose price is actually given by an expression similar to the Feynman–Kac formula. The exact calculation of the price is possible only in the particular case where the underlying process is a spectrally negative Lévy process (i.e., does not have positive jumps; see [5]), and practically more realistic cases of general Lévy processes, solutions to Lévy driven SDEs, etc. can be treated only numerically. To estimate the convergence rate of the respective Monte Carlo approximative methods, we need to estimate the accuracy of various approximation steps involved in the algorithm. In this paper, we focus ourselves on solving such a problem for discrete approximation of the integral functional of type (1).
For diffusion processes, this problem was studied in [4] and recently in [9] by means of the methods involving the particular structural features of the process, for example, the Malliavin calculus tools. On the other hand, in two recent papers [2, 3], an alternative method is developed, which exploits only the basic Markov structure of the process and the additive structure of the integral functional and its discrete approximations. One of the aims of this paper is to extend this method for a wider class of Markov processes. To explain our goal more in details, let us formulate our principal assumption on the process X.
The transition probability Pt(x,dy) of X admits a density pt(x,y) w.r.t. the Lebesgue measure on Rd. This density is differentiable w.r.t. t, and its derivative possesses the bound
|∂tpt(x,y)|≤BT,Xt−βqt,x(y),t≤T,
for some BT,X≥1, β≥1, and a measurable function q such that for any fixed t,x, the function qt,x(·) is a distribution density.
In [2, 3], the condition similar to X was formulated with β=1. Such a condition is satisfied for particularly important classes of diffusion processes and symmetric α-stable processes. However, in some natural cases, we can expect to get (2) only with β>11$]]>. As the simplest and the most illustrative example, we can take an α-stable process with drift:
Xt=ct+Zt,
where c≠0, and Z is an (e.g., symmetric) α-stable process. Then
pt(x,y)=t−d/αg(α)(y−x−ctt1/α),
where g(α) denotes the distribution density of Z1. Straightforward calculation shows that (2) holds with β=max(1,1/α), which is strictly greater than 1 when α<1. Since the Lévy noises are now extensively used in various applied models, the simple calculation made above shows that it is highly desirable to extend the results of [2] and [3], which deal with the “diffusive like” class of processes satisfying X with β=1, to the more general case of X with arbitrary β≥1.
Another aim of this paper is to develop tools that would allow us to get the bounds of the form (2) for a wider class of solutions of Lévy driven SDEs. One result of such a type is given in the recent preprint [7], with the process X being a solution of the SDE
dXt=b(Xt)dt+σ(Xt−)dZt,
where Z is a symmetric α-stable process. The method used therein is a version of the parametrix method, and it is quite sensitive to the form of the Lévy measure of the process Z on the entire Rd. Currently, apart from the stable noises, various types of “locally stable” noises are frequently used in applied models: tempered stable processes, damped stable processes, and so on. Heuristically, for a “locally stable” process, its “small jump behavior” is the same as for the stable one, whereas the “large jump behavior” may differ drastically and is determined by the “tail behavior” of the particular Lévy measure. Since (2) is genuinely related to the “local behavior” of the process, we can expect that the results of [7] should have a natural extension to the case of “locally stable” Z. However, it is a sophisticated problem to make such a conjecture rigorous; the main reason here is that the parametrix method treats the transition probability of a Lévy process as the “zero approximation” for the unknown transition probability density pt(x,y), and hence any bound for pt(x,y), which one may expect to design within this method, is at least as much complicated as the respective bound for the process Z. On the other hand, there is an extensive literature on the estimates of transition probability densities for Lévy processes (e.g., [1, 8, 6, 10, 12–17]; this list is far from being complete), which shows that these densities inherit the structure of the densities of the corresponding Lévy measures. In particular, in order to get exact two-sided bounds for pt(x,y), quite nontrivial structural assumptions on the “tails” of the Lévy measure should be imposed even in a comparatively simple “locally stable” case. Motivated by this observation on one hand and by the initial approximation problem which strongly relies on assumption (2) on the other hand, we pose the following general question: Is it possible to give a “rough” upper bound, which would be the same for a large class of transition probability densities of “locally stable processes”, without assuming complicated conditions on the “tails” of their Lévy measures? The answer is positive, and it roughly says that we can get the bound (2), where on the left-hand side, we have the transition probability density of the SDE driven by a “locally stable” process, and on the right-hand side, we have a (properly modified) transition probability density of an α-stable process. This bound is not necessarily precise: the power-type “tail” of the α-stable density might be essentially larger than the “tail,” for example, for an exponentially tempered α-stable law. The gain is, however, that under a mild set of assumptions, we obtain a uniform-in-class upper bound, useful in applications. To keep the exposition reasonably compact, we treat this problem in a comparatively simple case of a one-dimensional SDE of the form (10). The extension of these results to a more general multidimensional case is much more technical, and we postpone it to a separate publication.
The structure of the paper is the following. In Section 2, we formulate and prove our two main results concerning the accuracy of the strong and weak approximations of an integral functional by Riemannian sums, provided that condition X is satisfied. In Section 3, we outline a version of the parametrix method, which makes it possible to obtain (2) for solutions to Lévy driven SDEs without strong structural assumptions on the “tails” of the Lévy measure of the noise. In Section 4, an application for the price of an occupation time option is given.
Accuracy of discrete approximation for integral functionals
In this section, we prove two results. The first one concerns the “strong approximation rate”, that is, the control on the Lp-distance between IT(h) and its approximation IT,n(h).
This theorem extends [2, Theorem 2.1], where it was assumed that β=1.
The second result concerns the “weak approximation,” that is, the control on the difference between the expectations of certain terms, which involve IT(h) together with its approximation IT,n(h).
Suppose thatXholds. Then, for anyk∈Nand any bounded function f, we have|Ex(IT(h))kf(XT)−Ex(IT,n(h))kf(XT)|≤2β∨2k2BT,XTk+1‖h‖k‖f‖DT,β(n).
This theorem extends [3, Theorem 1.1], where it was assumed that β=1. In the proof, we will concentrate ourselves on the case β>11$]]>.
Using the Taylor expansion, we can directly obtain the following corollary of Theorem 2.2.
Suppose thatXholds, and let φ be an analytic function defined in a neighborhood of 0. In addition, suppose that the constantsDφ,Rφ>00$]]>are such that|φ(m)(0)m!|≤Dφ(1Rφ)m,m≥0.Then, for any bounded function f and a function h such thatT‖h‖<Rφ, we have the following bound:|Exφ(IT(h))f(XT)−Exφ(IT,n(h))f(XT)|≤CT,X,h,φ‖f‖DT,β(n),whereCT,X,h,φ=2β∨2DφBT,XT2‖h‖Rφ(1+T‖h‖Rφ)(1−T‖h‖Rφ)−3.
Before proceeding to the proof of Theorem 2.1, we give an auxiliary result this proof is based on. This result is, in fact, a weaker version of Theorem 2.2 with k=1 and f≡1, but we give it separately in order to make the exposition more transparent.
Let us introduce the notation used throughout the whole section: for t∈[kT/n,(k+1)T/n),k≥0, we put ηn(t)=kTn,ζn(t)=(k+1)Tn; that is, ηn(t) is the point of the partition {Tk/n,k≥0} of the time axis, closest to t from the left, and ζn(t) is the point closest to t from the right that is strictly greater than t.
We have
ExIT(h)−ExIT,n(h)=∫0TEx[h(Xs)−h(Xηn(s))]ds=∫0T∫Rdh(y)[ps(x,y)−pηn(s)(x,y)]dyds=M1+M2,
where
M1=∫0kn,βT/n∫Rdh(y)[ps(x,y)−pηn(s)(x,y)]dydsandM2=∫kn,βT/nT∫Rdh(y)[ps(x,y)−pηn(s)(x,y)]dyds
for some 1≤kn,β≤n, which will be chosen later. We estimate each term separately.
For M1, we have
|M1|≤‖h‖∫0kn,βT/n∫Rd[ps(x,y)+pηn(s)(x,y)]dyds=2‖h‖Tkn,βn.
Further, using (2), we get
|M2|≤‖h‖∫kn,βT/nT∫Rd|ps(x,y)−pηn(s)(x,y)|dyds≤‖h‖∫kn,βT/nT∫ηn(s)s∫Rd|∂upu(x,y)|dyduds≤BT,X‖h‖∫kn,βT/nT∫ηn(s)s∫Rdu−βqu,x(y)dyduds=BT,X‖h‖∫kn,βT/nT∫ηn(s)su−βduds=BT,X‖h‖∑i=kn,βn−1∫iT/n(i+1)T/n∫iT/nsu−βduds=BT,X‖h‖∑i=kn,βn−1∫iT/n(i+1)T/n∫u(i+1)T/nu−βdsdu≤TnBT,X‖h‖∑i=kn,βn−1∫iT/n(i+1)T/nu−βdu=TnBT,X‖h‖∫kn,βT/nTu−βdu.
Now we finalize the argument.
1) If β=1, put kn,β=1,n≥1. Then we get
|M1|≤2‖h‖Tn−1,|M2|≤TnBT,X‖h‖∫T/nTu−1du=BT,XT‖h‖n−1logn.
2) If β>11$]]>, put kn,β=[n1−1/β]+1. Then
|M1|≤2‖h‖T[n1−1/β]+1n≤2‖h‖Tn1−1/β+1n≤4‖h‖Tn−1/β.
To estimate M2, observe that
∫kn,βT/nTu−βdu≤T1−ββ−1(kn,βn)1−β≤T1−ββ−1(n1−1/βn)1−β≤T1−ββ−1n1−1/β.
Therefore,
|M2|≤TnBT,X‖h‖∫kn,βT/nTu−βdu≤BT,Xβ−1T2−β‖h‖n−1/β.
□
Since we can obtain the required bound for p<2 from the bound with p=2 by the Hölder inequality, we consider the case p≥2 only.
Define
Jt,n(h):=It(h)−It,n(h)=∫0tΔn(s)ds,Δn(s):=h(Xs)−h(Xηn(s)).
By definition, the function t↦Jt,n(h) is absolutely continuous. Then using the Newton–Leibniz formula twice, we get
|JT,n(h)|p=p(p−1)∫0T|Js,n(h)|p−2Δn(s)(∫sTΔn(t)dt)ds.
Therefore,
|JT,n(h)|p≤p(p−1)(HT,n,p1(h)+HT,n,p2(h)),
where
HT,n,p1(h)=∫0T|Js,n(h)|p−2|Δn(s)||∫sζn(s)Δn(t)dt|ds,HT,n,p2(h)=∫0T|Js,n(h)|p−2Δn(s)(∫ζn(s)TΔn(t)dt)ds.
Let us estimate separately the expectations of HT,n,p1(h) and HT,n,p2(h). By the Hölder inequality,
ExHT,n,p1(h)≤(Ex∫0T|Js,n(h)|pds)1−2/p×(Ex∫0T|Δn(s)|p/2|∫sζn(s)Δn(t)dt|p/2ds)2/p≤(Ex∫0T|Js,n(h)|pds)1−2/p((2‖h‖)p/2T(2‖h‖)p/2(Tn)p/2)2/p=4T1+2/pn−1‖h‖2(Ex∫0T|Js,n(h)|pds)1−2/p.
Further, observe that for every s the variables
Δn(s),|Js,n(h)|p−2Δn(s)
are Fζn(s)-measurable, where {Ft,t≥0} is the natural filtration for the process X. Hence,
ExHT,n,p2(h)=Ex(∫0T|Js,n(h)|p−2Δn(s)Ex(∫ζn(s)TΔn(t)dt|Fζn(s))ds)≤Ex(∫0T|Js,n(h)|p−2|Δn(s)||Ex(∫ζn(s)TΔn(t)dt|Fζn(s))|ds).
By Proposition 2.1 and the Markov property of X we have
|Ex(∫ζn(s)TΔn(t)dt|Fζn(s))|=|EXζn(s)∫0T−ζn(s)Δn(t)dt|≤5BT,XTDT,β(n)‖h‖.
Therefore, using the Hölder inequality, we get
ExHT,n,p2(h)≤5BT,XTDT,β(n)‖h‖×(Ex∫0T|Js,n(h)|pds)1−2/p(Ex∫0T|Δn(s)|p/2ds)2/p≤10BT,XT1+2/pDT,β(n)‖h‖2(Ex∫0T|Js,n(h)|pds)1−2/p.
Note that n−1≤DT,β(n); hence, the bounds for ExHT,n,p1(h) and ExHT,n,p2(h) finally yield the estimate
Ex|JT,n(h)|p≤14p(p−1)BT,XT1+2/pDT,β(n)‖h‖2×(Ex∫0T|Js,n(h)|pds)1−2/p.
It can be easily seen that this inequality also holds if JT,n(h) in the left-hand side is replaced by Jt,n(h). Integrating over t∈[0,T], we get
Ex∫0T|Jt,n(h)|pdt≤14p(p−1)BT,XT2+2/pDT,β(n)‖h‖2×(Ex∫0T|Js,n(h)|pds)1−2/p.
Because h is bounded, the left-hand side expression in this inequality is finite. Hence, resolving this inequality, we get
Ex∫0T|Js,n(h)|pds≤(14p(p−1)BT,X)p/2Tp+1(DT,β(n))p/2‖h‖p,
which together with (8) gives the required statement. □
We have
Ex[(IT(h))k−(IT,n(h))k]f(XT)=k!Ex∫Sk,0,T[h(Xs1)h(Xs2)⋯h(Xsk)−h(Xηn(s1))h(Xηn(s2))⋯h(Xηn(sk))]f(XT)∏i=1kdsi=k!∫Sk,0,T∫(Rd)k+1(∏i=1kh(yi))f(z)(∏i=1kpsi−si−1(yi−1,yi))×pT−sk(yk,z)dz∏j=1kdyj∏i=1kdsi−k!∫Sk,0,T∫(Rd)k+1(∏i=1kh(yi))f(z)×(∏i=1kpηn(si)−ηn(si−1)(yi−1,yi))pT−ηn(sk)(yk,z)dz∏j=1kdyj∏i=1kdsi=k!∑r=1k∫Sk,0,T∫(Rd)k+1(∏i=1kh(yi))f(z)Js1,…,sk,T(r)(x,y1,…,yk,z)×dz∏j=1kdyj∏i=1kdsi,
where the convention s0=0,sk+1=T,y0=x,yk+1=z is used, and the functions J(r),r=1,…,k are defined by the relations
Js1,…,sk,T(r)(x,y1,…,yk,z)=(∏i=1r−1pηn(si)−ηn(si−1)(yi−1,yi))×(psr−sr−1(yr−1,yr)−pηn(sr)−ηn(sr−1)(yr−1,yr))(∏i=rkpsi+1−si(yi,yi+1)).
Let us estimate the rth term in the last line in (9). We have
∫Sk,0,T∫(Rd)k+1(∏i=1kh(yi))f(z)Js1,…,sk,T(r)(x,y1,…,yk,z)dz∏j=1kdyj∏i=1kdsi≤‖h‖k‖f‖∫Sk,0,T∫(Rd)k+1|Js1,…,sk,T(r)(x,y1,…,yk,z)|dz∏j=1kdyj∏i=1kdsi.
Since the case β=1 was already treated in [3], for the rest of the proof, we assume that β>11$]]>.
Consider two cases: a) sr−sr−1>kn,βT/nk_{n,\beta }T/n$]]> and b) sr−sr−1≤kn,βT/n. Recall that kn,β=[n1−1/β]+1,β>11$]]>, is defined in the proof of Proposition 2.1.
In case a), using condition X and the Chapman–Kolmogorov equation, we derive
∫(Rd)k+1|Js1,…,sk,T(r)(x,y1,…,yk,z)|dz∏j=1kdyj≤BT,X∫(Rd)2pηn(sr−1)−ηn(s0)(x,yr−1)×|∫ηn(sr)−ηn(sr−1)sr−sr−1v−βqv,yr−1(yr)dv|dyr−1dyr.
Since kn,β≥2, in the case a) we have sr−sr−1≥2T/n, and hence
ηn(sr)−ηn(sr−1)≥sr−Tn−sr−1≥sr−sr−12.
Therefore, using the fact that qt,y(·) is the probability density for any t>00$]]> and y∈Rd, we finally get
∫(Rd)k+1|Js1,…,sk,T(r)(x,y1,…,yk,z)|dz∏j=1kdyj≤BT,X∫sr−sr−1−T/nsr−sr−1v−βdv≤BT,XTn(sr−sr−12)−β.
In case b), we simply apply the Chapman–Kolmogorov equation:
∫(Rd)k+1|Js1,…,sk,T(r)(x,y1,…,yk,z)|dz∏j=1kdyj≤∫(Rd)2pηn(sr−1)−ηn(s0)(x,yr−1)×(psr−sr−1(yr−1,yr)+pηn(sr)−ηn(sr−1)(yr−1,yr))dyr−1dyr≤2.
Therefore, summarizing the estimates obtained in cases a) and b) we get, using (7), the estimates
∫Sk,0,T∫(Rd)k+1|Js1,…,sk,T(r)(x,y1,…,yk,z)|dz∏j=1kdyj∏i=1kdsi≤BT,XTn2β∫0T∫0sr−kn,βT/nsr−1r−1(r−1)!(sr−sr−1)−β(T−sr)k−r(k−r)!dsr−1dsr+2∫0T∫sr−kn,βT/nsrsr−1r−1(r−1)!(T−sr)k−r(k−r)!dsr−1dsr≤BT,XTn2β∫0Tsrr−1(r−1)!(T−sr)k−r(k−r)!dsr(∫kn,βT/nTu−βdu)+2Tkn,βn∫0Tsrr−1(r−1)!(T−sr)k−r(k−r)!dsr≤2βBT,XTk+2−β(β−1)(k−1)!n−1/β+4Tk+1(k−1)!n−1/β≤2β∨2Tk+1BT,XDT,β(n)(k−1)!,
where in the fourth and fifth lines we used that sr−1r−1≤srr−1. Taking into account that in (9) we have k terms and the common multiplier k!, we finally arrive at (5). □
Condition <bold>X</bold> for solutions to Lévy driven SDEs
Consider the SDE
dXt=b(Xt)dt+dZt,X0=x,
where Z is a real-valued Lévy process. In [7], it was shown that if Zt is a symmetric α-stable process and b(·) is bounded and Lipschitz continuous, then the solution to Eq. (10) satisfies condition X with β=max(1,1/α) (in fact, in [7], more general multidimensional SDEs of the form (3) are considered). In this section, we outline the argument that makes it possible to extend the class of “Lévy noises”. Namely, we omit the requirement that Z is symmetric and relax the stability assumption, demanding Z to be “locally α-stable” in the sense we specify further.
Recall that the characteristic function of a Lévy process is of the form
EeiξZt=e−tψ(ξ),t>0,ξ∈R,0,\hspace{0.1667em}\xi \in \mathbb{R},\]]]>
where the characteristic exponent ψ admits the Lévy–Khinchin representation
ψ(ξ)=−iaξ+12σ2ξ2+∫R(1−eiξu+iξu1{|u|≤1})μ(du).
In what follows, we assume that σ2=0 and the Lévy measure μ is of the form
μ(du)=C+u−1−α1u∈(0,1)du+C−|u|−1−α1u∈(−1,0)du+m(u)du
with some C±≥0 and m(u)≥0 such that m(u)=0 for |u|≤1 and
m(u)≤c|u|−1−α,|u|≥1.
On the interval [−1,1], the Lévy measure μ given by (12) coincides with the Lévy measure of a (nonsymmetric) α-stable process. This is the reason for us to call Z a “locally α-stable” process: its “local behavior” near the origin is similar to those of the α-stable process. In that context, condition (13) means that the “tails” of the Lévy measure for μ are dominated by the “tails” of the α-stable Lévy measure.
Let us impose three minor conventions, which will simplify the technicalities. First, since we are mostly interested in the case β>11$]]>, we assume that α<1. Second, the latter assumption assures that the integral
∫{|u|≤1}uμ(du)
is well defined, and we assume that the constant a in (11) equals this integral; that is, ψ has the form
ψ(ξ)=∫R(1−eiξu)μ(du).
Clearly, this does not restrict the generality because we can change the constant a by changing appropriately the drift coefficient b(·) in (10). Finally, in order to avoid the usage of the Rademacher theorem (see [7, Lemma 7.4] for the case where b is just Lipschitz continuous), let us assume that b∈C1(R).
In what follows, we show how the parametrix construction developed in [7] can be modified to provide the representation and the bounds for the transition probability density pt(x,y) of the solution to (10) driven by the “locally stable” noise Z.
Let us introduce some notation and give some preliminaries. We denote the space and the space-time convolutions respectively by
(f∗g)(x,y):=∫Rdf(x,z)g(z,y)dz,(f⊛g)t(x,y):=∫0t(ft−s∗gs)(x,y)ds=∫0t∫Rdft−s(x,z)gs(z,y)dzds.
Generically, the parametrix construction provides a representation of the required transition probability density in the form
pt(x,y)=pt0(x,y)+∫0t∫Rpt−s0(x,z)Ψs(z,y)dzds,t>0,x,y∈R.0,\hspace{2.5pt}\hspace{2.5pt}x,y\in \mathbb{R}.\]]]>
Here pt0(x,y) is a “zero approximation term” for the unknown pt(x,y), and the function Ψt(x,y) is given by the “convolution series”
Ψt(x,y)=∑k=1∞Φt⊛k(x,y),t>0,x,y∈R,0,\hspace{2.5pt}\hspace{2.5pt}x,y\in \mathbb{R},\]]]>
where the function Φt(x,y) depends on the particular choice of pt0(x,y) and equals
Φt(x,y):=(Lx−∂t)pt0(x,y),
where
Lf(x):=b(x)f′(x)+∫R(f(x+u)−f(x))μ(du),f∈Cb2(R),
is the formal generator of the process X. The subscript x means that the operator L is applied with respect to the variable x. Note that to make this construction feasible, we should properly choose the “zero approximation term” pt0(x,y), so that the convolution series (15) converges and the space-time convolution in (14) is well defined. To introduce in our setting such pt0(x,y) and then to construct the bounds for the associated Φt(x,y) and its convolution powers, we need some more notation.
Denote by Z(α,C±) the α-stable process with the Lévy measure μα,C±(du)=m(α,C±)(u)du,
m(α,C±)(u):=C+u−1−α1u>0du+C−(−u)−1−α1u<0,0}du+C_{-}{(-u)}^{-1-\alpha }\mathbb{1}_{u<0},\]]]>
and the characteristic exponent
ψ(α,C±)(ξ)=∫R(1−eiξu)μ(α,C±)(du).
Note that since
ψ(α,C±)(cξ)=cαψ(α,C±)(ξ),c>0,0,\]]]>
the process Z(α,C±) possesses the scaling property
Law(Zct(α,C±))=Law(c1/αZt(α,C±)),c>0.0.\]]]>
Denote by gt(α,C±) the distribution density of Zt(α,C±). By the scaling property we have
gt(α,C±)(x)=t−1/αg(α,C±)(xt−1/α),g(α,C±):=g1(α,C±).
Denote also by Z(α) the symmetric α-stable process, that is, the process of the form introduced before with C+=C−=1. Let gt(α) be the respective distribution density, and g(α):=g1(α).
Finally, denote by χt(x) and θt(y), respectively, the solutions to the ODEs
dχt=b(χt)dt,χ0=x,dθt=−b(θt)dt,θ0=y.
Note that these solutions exist because b(·) is Lipschitz continuous.
Now we are ready to formulate the main statement of this section.
Letpt0(x,y):=gt(α,C±)(θt(y)−x).Then the convolution series (15) is well defined, and Eq. (14) gives a representation of the transition probability densitypt(x,y)of the process X. This density and its time derivative have the following upper bounds:pt(x,y)≤BT,X(gt+1(α)+gt(α))(y−χt(x)),t∈(0,T],x,y∈R,∂tpt(x,y)≤BT,Xt−1/α(gt+1(α)+gt(α))(y−χt(x)),t∈(0,T],x,y∈R.Consequently, the process X satisfies conditionXwithβ=1/αandqt,x(y)=(gt+1(α)+gt(α))(y−χt(x)).
First, we evaluate Φt(x,y). If not stated otherwise, in all our estimates, we further assume that t∈(0,T] for some T>00$]]> and that x,y∈R. Observe that g(α,C±)∈Cb2(R). Indeed, this property easily follows from the Fourier inversion formula and the expression for the characteristic function. It is known that gt(α,C±)(y−x) is the fundamental solution to the Cauchy problem for ∂t−L(α,C±), where L(α,C±) denotes the generator of the process Z(α,C±):
L(α,C±)f(x)=∫R(f(x+u)−f(x))μ(α,C±)(du),f∈Cb2(R).
Since
(∂t−Lx(α,C±))gt(α,C±)(y−x)=0,
we have
∂tpt0(x,y)=[∂tgt(α,C±)(w)+∂tθt(y)t2/α(g(α,C±))′(wt1/α)]|w=θt(y)−x=[1t1/α(L(α,C±)g(α,C±))(wt1/α)+∂tθt(y)t2/α(g(α,C±))′(wt1/α)]|w=θt(y)−x,
where in the last identity we used the scaling property of gt(α,C±) and the fact that L(α,C±) is a homogeneous operator of order 1/α. Next, by the definition of L and L(α,C±) we get
Lxpt0(x,y)=[1t1/α(L(α,C±)g(α))(wt1/α)−b(x)t2/α(g(α,C±))′(wt1/α)]|w=θt(y)−x+[∫|u|≥1(1t1/αg(α,C±)(w−ut1/α)−1t1/αg(α,C±)(wt1/α))×(m(u)du−m(α,C±)(du))]|w=θt(y)−x.
Therefore, using the relation ∂tθt(y)=−b(θt(y)), we get
Φt(x,y)=(Lx−∂t)pt0(x,y)=b(x)−b(θt(y))t2/α(g(α,C±))′(θt(y)−xt1/α)+1t1/α∫|u|≥1(g(α,C±)(θt(y)−x−ut1/α)−g(α,C±)(θt(y)−xt1/α))×(m(u)−m(α,C±)(u))du=:Φt1(x,y)+Φt2(x,y).
Further, we find the bounds for the absolute values of Φt1(x,y), Φt2(x,y), and Φt(x,y). In what follows, D denotes a generic constant whose value may be different from place to place. We have g(α,C±)(x)≤Dg(α)(x),x∈R,|(g(α,C±))′(x)|≤D(1+|x|)−1g(α)(x),x∈R,|(g(α,C±))″(x)|≤D(1+|x|)−2g(α)(x),x∈R. Since the argument used in the proof of (23)–(25) is quite standard (see, e.g., [7, Appendix A]), we omit the details.
By (24) and the Lipschitz continuity of b(·) we have
|Φt1(x,y)|≤D|x−θt(y)|t2/α|(g(α,C±))′(θt(y)−xt1/α)|≤Dt1/αg(α)(θt(y)−xt1/α)=Dgt(α)(θt(y)−x).
To get the estimate for |Φt2(x,y)|, we first observe that
|m(u)−m(α,C±)(u)|I|u|≥1≤Dg(α)(u),
which implies
|Φt2(x,y)|≤Dt1/α∫|u|≥1g(α,C±)(θt(y)−x−ut1/α)g(α)(u)du+Dt1/αg(α,C±)(θt(y)−xt1/α).
Taking into account (23), we deduce that
|Φt2(x,y)|≤D(gt(α)∗g1(α)+gt(α))(θt(y)−x)=D(gt+1(α)(x)+gt(α)(x))(θt(y)−x).
Combining the estimates for Φt1(x,y) and Φt2(x,y), we get
|Φt(x,y)|≤D(gt+1(α)+gt(α))(θt(y)−x).
Our next step is to estimate the convolution powers of Φ. It is shown in [7, Appendix B] that the kernel gt(α)(θt(y)−x) possess the following sub-convolution property:
∫Rgt−s(α)(θt−s(z)−x)gs(α)(θs(y)−z)dz≤Dgt(α)(θt(y)−x).
By this property we get
∫Rgt+1−s(α)(θt−s(z)−x)gs(α)(θs(y)−z)dz≤Dgt+1(α)(θt(y)−x),∫Rgt−s+1(α)(θt−s(z)−x)gs+1(α)(θs(y)−z)dz≤Dgt+2(α)(θt(y)−x)≤Dgt+1(α)(θt(y)−x),
where in the last line we used that g2(α)≤Dg1, and therefore gt+2(α)=gt(α)∗g2(α)≤Dgt+2(α). Having these estimates, we deduce in the same way as in [7, Section 3] that
|Φt⊛k(x,y)|≤C0(Dt)k−1k!(gt+1(α)+gt(α))(θt(y)−x).
Therefore, the series (15) converges absolutely for (t,x,y)∈(0,∞)×R×R, and
|Ψt(x,y)|≤D(gt+1(α)+gt(α))(θt(y)−x).
Applying once again the sub-convolution property (27), we see that the convolution p0⊛Ψ is well defined and
|(p0⊛Ψ)t(x,y)|≤D(gt+1(α)+gt(α))(θt(y)−x).
Thus, expression (14) for pt(x,y) is well defined for any (t,x,y)∈(0,∞)×R×R, and
|pt(x,y)|≤D(gt+1(α)+gt(α))(θt(y)−x).
Finally, to get (19), we use the following inequalities, which were proved in [7, Appendix B]:
c|θt(y)−x|≤|χt(x)−y|≤C|θt(y)−x|.
Since for any constant c>00$]]>, we have gt(α)(x)≍gt(α)(cx), this completes the proof of (19).
Our final step is to use representation (14) in order to find bounds for ∂tpt(x,y). Since pt0(x,y) and Φt(x,y) are given explicitly, it is straightforward to show that these functions are differentiable with respect to t and to check using (23)–(25) that |∂tpt0(x,y)|≤Dt−1/αgt(α)(θt(y)−x),|∂tΦt(x,y)|≤Dt−1/α(gt+1(α)+gt(α))(θt(y)−x).
To show that the convolution powers Φt⊛k(x,y) are differentiable in t and to get the upper bounds, we use the following trick. The expression for Φt⊛(k+1)(x,y) can be reorganized as follows:
Φt⊛(k+1)(x,y)=∫0t∫RΦt−s⊛k(x,z)Φs(z,y)dzds=∫0t/2∫RΦt−s⊛k(x,z)Φs(z,y)dzds+∫0t/2∫RΦs⊛k(x,z)Φt−s(z,y)dzds.
If k=1, then the first line in (33) does not allow us to differentiate ∂tΦt⊛(2)(x,y) because the upper bound for ∂tΦt−s(x,z) has a nonintegrable singularity (t−s)−1/α at the vicinity of the point s=t (recall that α<1). However, the identity given by the second line in (33) does not contain such singularities, and we can show using induction that for any k≥1, the function Φt⊛k(x,y) is continuously differentiable in t, satisfies
∂tΦt⊛(k+1)(x,y)=∫0t/2∫Rd(∂tΦ⊛k)t−s(x,z)Φs(z,y)dzds+∫0t/2∫RdΦs⊛k(x,z)(∂tΦ)t−s(z,y)dzds+∫RdΦt/2⊛k(x,z)Φt/2(z,y)dz,
and possesses the bound
|∂tΦt⊛k(x,y)|≤C0(Dt)k−1t−1/αk!(gt+1(α)+gt(α))(θt(y)−x),k≥1.
Since the proof is completely analogous to the proof of [7, Lemma 7.3], we omit the details.
From (34) we derive the following bound for the derivative of Ψt(x,y):
|∂tΨt⊛k(x,y)|≤Dt−1/α(gt+1(α)+gt(α))(θt(y)−x).
Reorganizing representation (14) in the same way as in (33), we get
pt(x,y)=pt0(x,y)+∫0t/2∫Rdpt−s0(x,z)Ψs(z,y)dzds+∫0t/2∫Rdps0(x,z)Ψt−s(z,y)dzds.
Using this representation of pt(x,y) together with (31) and (35), we derive the existence of the continuous derivative ∂tpt(x,y), which satisfies the inequality
|∂tpt(x,y)|≤Dt−1/α(gt+1(α)+gt(α))(θt(y)−x).
Using estimates (30) in the same way as in the proof of (19), we can change the argument θt(y)−x in the right-hand side of this estimate to y−χt(x), which completes the proof of (20). □
Application: the price of an occupation time option
In this section, we consider an occupation time option (see [11]), with the price of the option depending on the time spent by an asset price process in a given set. Comparing to the standard barrier options, which are activated or cancelled when the asset price process hits some definite level (barrier), the payoff of the occupation time option depends on the time during which the asset price process stays below or above such a barrier.
For instance, for the strike price K, the barrier level L and the knock-out rate ρ, the payoff of a down-and-out call occupation time option equals
exp(−ρ∫0TI{St≤L}dt)(ST−K)+,
and the price C(T) of the option is defined as
C(T)=exp(−rT)E[exp(−ρ∫0TI{St≤L}dt)(ST−K)+]
where r is the risk-free interest rate (see [11]).
Assume that the price of an asset S={St,t≥0} is of the form
St=S0exp(Xt),
where X is the Markov process studied in the previous sections. Then the time spent by the process S in a set J⊂R equals the time spent by X in the set J′=logJ.
Let us approximate the price C(T) of our option by
Cn(T)=exp(−rT)E[exp(−ρT/n∑k=0n−1I{SkT/n≤L}dt)(ST−K)+].
Then, using the results from the previous sections, we can get the control on the accuracy of such an approximation.
First, we apply Theorem 2.1 and derive the strong approximation rate.
Suppose that the process X satisfies conditionXand assume that there existsλ>11$]]>such thatG(λ):=Eexp(λXT)=E(ST)λ<+∞.
Then|Cn(T)−C(T)|≤exp(−rT)ρG(λ)1/λCT,λ/(λ−1)(DT,β(n))1/2with constantsCT,λ/(λ−1)andDT,β(n)given by (4).
The proof is a simple corollary of Theorem 2.1. Denote h(x)=ρIx≤logL. Then, keeping the notation of Section 2, we get
C(T)=e−rTEe−IT(h)(ST−K)+,Cn(T)=e−rTEe−IT,n(h)(ST−K)+.
By the Hölder inequality with p=λ and q=λ/(λ−1),
|Cn(T)−C(T)|≤e−rT(E(ST−K)+λ)1/λ×(E|e−IT(h)−e−IT,n(h)|λ/(λ−1))(λ−1)/λ.
Since |e−a−e−b|≤|a−b| for positive a and b, we have
|Cn(T)−C(T)|≤e−rTG(λ)1/λ(E|IT(h)−IT,n(h)|λ/(λ−1))(λ−1)/λ,
and thus the required statement follows from Theorem 2.1 with p=λλ−1. □
We also can control the accuracy of the approximation using the weak rate bound from Theorem 2.2. Note that the bound given in the next proposition is sharper than those obtained in the previous proposition precisely when λ>22$]]>.
Under the assumptions of Proposition4.1, we have|Cn(T)−C(T)|≤2β∨2+1max{BT,XρT2(1+ρT)eρT,G(λ)}e−rTD˜T,β(n),whereD˜T,β(n)=n−(1−1/λ)logn,β=1,max(1,T1−ββ−1)n−1/β(1−1/λ),β>1.1.\end{array}\right.\]]]>
For N>00$]]>, denote
CN(T)=exp(−rT)E[exp(−ρ∫0TI{St≤L}dt)((ST−K)+∧N)],CnN(T)=exp(−rT)E[exp(−ρT/n∑k=0n−1I{SkT/n≤L}dt)((ST−K)+∧N)].
Then
|Cn(T)−C(T)|≤|CnN(T)−CN(T)|+|C(T)−CN(T)|+|Cn(T)−CnN(T)|.
We estimate each term separately.
Using Theorem 2.2 and the Taylor expansion of the exponent, we derive
|CnN(T)−CN(T)|≤2β∨2BT,XNTexp(−rT)∑k=1∞ρkk!k2TkDT,β(n)=2β∨2BT,XρNT2(1+ρT)exp(ρT−rT)DT,β(n).
For the last two terms, we get
|C(T)−CN(T)|+|Cn(T)−CnN(T)|≤2exp(−rT)E[(ST−K)+−(ST−K)+∧N]≤2exp(−rT)E[STI{ST>N}]=2exp(−rT)E[STNλ−1I{ST>N}Nλ−1]≤2G(λ)Nλ−1exp(−rT).N\}}]=2\exp (-rT)E\bigg[\frac{S_{T}{N}^{\lambda -1}\mathbb{I}_{\{S_{T}>N\}}}{{N}^{\lambda -1}}\bigg]\\{} & \displaystyle \hspace{1em}\le \frac{2G(\lambda )}{{N}^{\lambda -1}}\exp (-rT).\end{array}\]]]>
To complete the proof, put N=n1/(βλ). □
Acknowledgments
The second- and third-named authors gratefully acknowledge the DFG Grant Schi 419/8-1; the third-named author acknowledges the joint DFFD-RFFI project No. 09-01-14.
ReferencesBogdan, K., Grzywny, T., Ryznar, M.: Density and tails of unimodal convolution semigroups. J. Funct. Anal.266(6), 3543–3571 (2014). MR3165234. doi:10.1016/j.jfa.2014.01.007Ganychenko, I., Kulik, A.: Rates of approximation of nonsmooth integral-type functionals of Markov processes. Mod. Stoch., Theory Appl.2, 117–126 (2014). MR3316480. doi:10.15559/vmsta-2014.12Ganychenko, I., Kulik, A.: Weak approximation rates for integral functionals of Markov processes. Mod. Stoch., Theory Appl.3, 251–266 (2015) Gobet, E., Labart, C.: Sharp estimates for the convergence of the density of the Euler scheme in small time. Electron. Commun. Probab.13, 352–363 (2008). MR2415143. doi:10.1214/ECP.v13-1393Guerin, H., Renaud, J.-F.: Joint distribution of a spectrally negative Lévy process and its occupation time, with step option pricing in view. arXiv:1406.3130Knopova, V.: Compound kernel estimates for the transition probability density of a Lévy process in Rn. Theory Probab. Math. Stat.89, 57–70 (2014). MR3235175. doi:10.1090/s0094-9000-2015-00935-2Knopova, V., Kulik, A.: Parametrix construction of the transition probability density of the solution to an SDE driven by α-stable noise. arXiv:1412.8732Knopova, V., Kulik, A.: Intrinsic small time estimates for distribution densities of Lévy processes. Random Oper. Stoch. Equ.21(4), 321–344 (2013). MR3139314. doi:10.1515/rose-2013-0015Kohatsu-Higa, A., Makhlouf, A., Ngo, H.L.: Approximations of non-smooth integral type functionals of one dimensional diffusion precesses. Stoch. Process. Appl.124, 1881–1909 (2014). MR3170228. doi:10.1016/j.spa.2014.01.003Kulczycki, T., Ryznar, M.: Gradient estimates of harmonic functions and transition densities for Lévy processes. arXiv:1307.7158. MR3413864. doi:10.1090/tran/6333Linetsky, V.: Step options. Math. Finance9, 55–96 (1999). MR1849356. doi:10.1111/1467-9965.00063Mimica, A.: Heat kernel upper estimates for symmetric jump processes with small jumps of high intensity. Potential Anal.36(2), 203–222 (2012). MR2886459. doi:10.1007/s11118-011-9225-1Pruitt, W.E., Taylor, S.J.: The potential kernel and hitting probabilities for the general stable process in Rn. Trans. Am. Math. Soc.146, 299–321 (1969). MR0250372Rosinski, J., Singlair, J.: Generalized tempered stable processes. Stab. Probab.90, 153–170 (2010). MR2798857. doi:10.4064/bc90-0-10Sztonyk, P.: Estimates of tempered stable densities. J. Theor. Probab.23(1), 127–147 (2010). MR2591907. doi:10.1007/s10959-009-0208-8Sztonyk, P.: Transition density estimates for jump Lévy processes. Stoch. Process. Appl.121, 1245–1265 (2011). MR2794975. doi:10.1016/j.spa.2011.03.002Watanabe, T.: Asymptotic estimates of multi-dimensional stable densities and their applications. Trans. Am. Math. Soc.359(6), 2851–2879 (2007). MR2286060. doi:10.1090/S0002-9947-07-04152-9