VMSTA Modern Stochastics: Theory and Applications 2351-6054 2351-6046 2351-6046 VTeXMokslininkų g. 2A, 08412 Vilnius, Lithuania VMSTA46 10.15559/15-VMSTA46 Research Article Accuracy of discrete approximation for integral functionals of Markov processes GanychenkoIuriiiurii_ganychenko@ukr.neta KnopovaVictoriavicknopova@googlemail.comb KulikAlexeikulik.alex.m@gmail.comc Taras Shevchenko National University of Kyiv, Kyiv, Ukraine V.M. Glushkov Institute of Cybernetics NAS of Ukraine, Kyiv, Ukraine Institute of Mathematics, National Academy of Sciences of Ukraine, Kyiv, Ukraine Corresponding author. 2015 3012201524401420 4122015 15122015 15122015 © 2015 The Author(s). Published by VTeX2015 Open access article under the CC BY license.

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.

Markov process integral functional approximation rate 60H07 60H35
Introduction

Let Xt , t0 , be a Markov process with values in Rd . Consider an integral functional of the form IT(h)=0Th(Xt)dt, where h:RdR 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=0n1h(X(kT)/n),n1. 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 AB(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 , 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 ), 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  and recently in  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 , 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),tT, for some BT,X1 , β1 , and a measurable function q such that for any fixed t,x , the function qt,x(·) is a distribution density.

In , 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 β>1 1$]]>. As the simplest and the most illustrative example, we can take an α-stable process with drift: Xt=ct+Zt, where c0 , and Z is an (e.g., symmetric) α-stable process. Then pt(x,y)=td/αg(α)(yxctt1/α), 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  and , 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 , 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  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., ; 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) . Suppose that X holds. Then, for any p>0 0$]]>, (Ex|IT(h)IT,n(h)|p)1/pCT,ph(DT,β(n))1/2, where h=supx|h(x)| , DT,β(n)=n1logn,β=1,max(1,T1ββ1)n1/β,β>1,CT,p=(14p(p1)BT,X)1/2T,p2,CT,2=(28BT,X)1/2T,p(0,2). 1,\end{array}\right.\\{} & \displaystyle C_{T,p}=\left\{\begin{array}{l@{\hskip10.0pt}l}{(14p(p-1)B_{T,X})}^{1/2}T,\hspace{1em}& p\ge 2,\\{} C_{T,2}={(28B_{T,X})}^{1/2}T,\hspace{1em}& p\in (0,2).\end{array}\right.\end{array}\]]]>

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 that X holds. Then, for any kN and any bounded function f, we have |Ex(IT(h))kf(XT)Ex(IT,n(h))kf(XT)|2β2k2BT,XTk+1hkfDT,β(n).

This theorem extends [3, Theorem 1.1], where it was assumed that β=1 . In the proof, we will concentrate ourselves on the case β>1 1$]]>. Using the Taylor expansion, we can directly obtain the following corollary of Theorem 2.2. Suppose that X holds, and let φ be an analytic function defined in a neighborhood of 0. In addition, suppose that the constants Dφ,Rφ>0 0$]]> are such that |φ(m)(0)m!|Dφ(1Rφ)m,m0. Then, for any bounded function f and a function h such that Th<Rφ , we have the following bound: |Exφ(IT(h))f(XT)Exφ(IT,n(h))f(XT)|CT,X,h,φfDT,β(n), where CT,X,h,φ=2β2DφBT,XT2hRφ(1+ThRφ)(1ThRφ)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 f1 , but we give it separately in order to make the exposition more transparent.

Suppose that X holds. Then |ExIT(h)ExIT,n(h)|5BT,XThDT,β(n).

Let us introduce the notation used throughout the whole section: for t[kT/n,(k+1)T/n),k0 , we put ηn(t)=kTn,ζn(t)=(k+1)Tn ; that is, ηn(t) is the point of the partition {Tk/n,k0} 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=0TRdh(y)[ps(x,y)pηn(s)(x,y)]dyds=M1+M2, where M1=0kn,βT/nRdh(y)[ps(x,y)pηn(s)(x,y)]dydsandM2=kn,βT/nTRdh(y)[ps(x,y)pηn(s)(x,y)]dyds for some 1kn,βn , which will be chosen later. We estimate each term separately.

For M1 , we have |M1|h0kn,βT/nRd[ps(x,y)+pηn(s)(x,y)]dyds=2hTkn,βn. Further, using (2), we get |M2|hkn,βT/nTRd|ps(x,y)pηn(s)(x,y)|dydshkn,βT/nTηn(s)sRd|upu(x,y)|dydudsBT,Xhkn,βT/nTηn(s)sRduβqu,x(y)dyduds=BT,Xhkn,βT/nTηn(s)suβduds=BT,Xh i=kn,βn1iT/n(i+1)T/niT/nsuβduds=BT,Xh i=kn,βn1iT/n(i+1)T/nu(i+1)T/nuβdsduTnBT,Xh i=kn,βn1iT/n(i+1)T/nuβdu=TnBT,Xhkn,βT/nTuβdu.

Now we finalize the argument.

1) If β=1 , put kn,β=1,n1 . Then we get |M1|2hTn1,|M2|TnBT,XhT/nTu1du=BT,XThn1logn.

2) If β>1 1$]]>, put kn,β=[n11/β]+1 . Then |M1|2hT[n11/β]+1n2hTn11/β+1n4hTn1/β. To estimate M2 , observe that kn,βT/nTuβduT1ββ1(kn,βn)1βT1ββ1(n11/βn)1βT1ββ1n11/β. Therefore, |M2|TnBT,Xhkn,βT/nTuβduBT,Xβ1T2βhn1/β. □ 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 p2 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 tJt,n(h) is absolutely continuous. Then using the Newton–Leibniz formula twice, we get |JT,n(h)|p=p(p1)0T|Js,n(h)|p2Δn(s)(sTΔn(t)dt)ds. Therefore, |JT,n(h)|pp(p1)(HT,n,p1(h)+HT,n,p2(h)), where HT,n,p1(h)=0T|Js,n(h)|p2|Δn(s)||sζn(s)Δn(t)dt|ds,HT,n,p2(h)=0T|Js,n(h)|p2Δ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)(Ex0T|Js,n(h)|pds)12/p×(Ex0T|Δn(s)|p/2|sζn(s)Δn(t)dt|p/2ds)2/p(Ex0T|Js,n(h)|pds)12/p((2h)p/2T(2h)p/2(Tn)p/2)2/p=4T1+2/pn1h2(Ex0T|Js,n(h)|pds)12/p. Further, observe that for every s the variables Δn(s),|Js,n(h)|p2Δn(s) are Fζn(s) -measurable, where {Ft,t0} is the natural filtration for the process X. Hence, ExHT,n,p2(h)=Ex(0T|Js,n(h)|p2Δn(s)Ex(ζn(s)TΔn(t)dt|Fζn(s))ds)Ex(0T|Js,n(h)|p2|Δ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×(Ex0T|Js,n(h)|pds)12/p(Ex0T|Δn(s)|p/2ds)2/p10BT,XT1+2/pDT,β(n)h2(Ex0T|Js,n(h)|pds)12/p. Note that n1DT,β(n) ; hence, the bounds for ExHT,n,p1(h) and ExHT,n,p2(h) finally yield the estimate Ex|JT,n(h)|p14p(p1)BT,XT1+2/pDT,β(n)h2×(Ex0T|Js,n(h)|pds)12/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 Ex0T|Jt,n(h)|pdt14p(p1)BT,XT2+2/pDT,β(n)h2×(Ex0T|Js,n(h)|pds)12/p. Because h is bounded, the left-hand side expression in this inequality is finite. Hence, resolving this inequality, we get Ex0T|Js,n(h)|pds(14p(p1)BT,X)p/2Tp+1(DT,β(n))p/2hp, which together with (8) gives the required statement. □ Denote Sk,a,b:={(s1,s2,,sk)Rk:as1<s2<<skb},kN,a,bR. We have Ex[(IT(h))k(IT,n(h))k]f(XT)=k!ExSk,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=1kpsisi1(yi1,yi))×pTsk(yk,z)dz j=1kdyj i=1kdsik!Sk,0,T(Rd)k+1( i=1kh(yi))f(z)×( i=1kpηn(si)ηn(si1)(yi1,yi))pTηn(sk)(yk,z)dz j=1kdyj i=1kdsi=k! r=1kSk,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=1r1pηn(si)ηn(si1)(yi1,yi))×(psrsr1(yr1,yr)pηn(sr)ηn(sr1)(yr1,yr))( i=rkpsi+1si(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=1kdsihkfSk,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 , for the rest of the proof, we assume that β>1 1$]]>.

Consider two cases: a) srsr1>kn,βT/n k_{n,\beta }T/n$]]> and b) srsr1kn,βT/n . Recall that kn,β=[n11/β]+1,β>1 1$]]>, 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=1kdyjBT,X(Rd)2pηn(sr1)ηn(s0)(x,yr1)×|ηn(sr)ηn(sr1)srsr1vβqv,yr1(yr)dv|dyr1dyr. Since kn,β2 , in the case a) we have srsr12T/n , and hence ηn(sr)ηn(sr1)srTnsr1srsr12. Therefore, using the fact that qt,y(·) is the probability density for any t>0 0$]]> and yRd , we finally get (Rd)k+1|Js1,,sk,T(r)(x,y1,,yk,z)|dz j=1kdyjBT,Xsrsr1T/nsrsr1vβdvBT,XTn(srsr12)β. 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(sr1)ηn(s0)(x,yr1)×(psrsr1(yr1,yr)+pηn(sr)ηn(sr1)(yr1,yr))dyr1dyr2. 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=1kdsiBT,XTn2β0T0srkn,βT/nsr1r1(r1)!(srsr1)β(Tsr)kr(kr)!dsr1dsr+20Tsrkn,βT/nsrsr1r1(r1)!(Tsr)kr(kr)!dsr1dsrBT,XTn2β0Tsrr1(r1)!(Tsr)kr(kr)!dsr(kn,βT/nTuβdu)+2Tkn,βn0Tsrr1(r1)!(Tsr)kr(kr)!dsr2βBT,XTk+2β(β1)(k1)!n1/β+4Tk+1(k1)!n1/β2β2Tk+1BT,XDT,β(n)(k1)!, where in the fourth and fifth lines we used that sr1r1srr1 . 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 , 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 , 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=etψ(ξ),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(1eiξ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+u1α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 β>1 1$]]>, 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(1eiξ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 bC1(R) .

In what follows, we show how the parametrix construction developed in  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 (fg)(x,y):=Rdf(x,z)g(z,y)dz,(fg)t(x,y):=0t(ftsgs)(x,y)ds=0tRdfts(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)+0tRpts0(x,z)Ψs(z,y)dzds,t>0,x,yR. 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Φtk(x,y),t>0,x,yR, 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):=(Lxt)pt0(x,y), where Lf(x):=b(x)f(x)+R(f(x+u)f(x))μ(du),fCb2(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+u1α1u>0du+C(u)1α1u<0, 0}du+C_{-}{(-u)}^{-1-\alpha }\mathbb{1}_{u<0},\]]]> and the characteristic exponent ψ(α,C±)(ξ)=R(1eiξ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)=t1/αg(α,C±)(xt1/α),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.

Let pt0(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 density pt(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,yR, tpt(x,y)BT,Xt1/α(gt+1(α)+gt(α))(yχt(x)),t(0,T],x,yR. Consequently, the process X satisfies condition X with β=1/α and qt,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>0 0$]]> and that x,yR . 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±)(yx) is the fundamental solution to the Cauchy problem for tL(α,C±) , where L(α,C±) denotes the generator of the process Z(α,C±) : L(α,C±)f(x)=R(f(x+u)f(x))μ(α,C±)(du),fCb2(R). Since (tLx(α,C±))gt(α,C±)(yx)=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±)(wut1/α)1t1/αg(α,C±)(wt1/α))×(m(u)dum(α,C±)(du))]|w=θt(y)x. Therefore, using the relation tθt(y)=b(θt(y)) , we get Φt(x,y)=(Lxt)pt0(x,y)=b(x)b(θt(y))t2/α(g(α,C±))(θt(y)xt1/α)+1t1/α|u|1(g(α,C±)(θt(y)xut1/α)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),xR, |(g(α,C±))(x)|D(1+|x|)1g(α)(x),xR, |(g(α,C±))(x)|D(1+|x|)2g(α)(x),xR. 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|1Dg(α)(u), which implies |Φt2(x,y)|Dt1/α|u|1g(α,C±)(θt(y)xut1/α)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: Rgts(α)(θts(z)x)gs(α)(θs(y)z)dzDgt(α)(θt(y)x). By this property we get Rgt+1s(α)(θts(z)x)gs(α)(θs(y)z)dzDgt+1(α)(θt(y)x),Rgts+1(α)(θts(z)x)gs+1(α)(θs(y)z)dzDgt+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 |Φtk(x,y)|C0(Dt)k1k!(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>0 0$]]>, 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)|Dt1/αgt(α)(θt(y)x), |tΦt(x,y)|Dt1/α(gt+1(α)+gt(α))(θt(y)x).

To show that the convolution powers Φtk(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)=0tRΦtsk(x,z)Φs(z,y)dzds=0t/2RΦtsk(x,z)Φs(z,y)dzds+0t/2RΦsk(x,z)Φts(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Φts(x,z) has a nonintegrable singularity (ts)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 k1 , the function Φtk(x,y) is continuously differentiable in t, satisfies tΦt(k+1)(x,y)=0t/2Rd(tΦk)ts(x,z)Φs(z,y)dzds+0t/2RdΦsk(x,z)(tΦ)ts(z,y)dzds+RdΦt/2k(x,z)Φt/2(z,y)dz, and possesses the bound |tΦtk(x,y)|C0(Dt)k1t1/αk!(gt+1(α)+gt(α))(θt(y)x),k1. 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Ψtk(x,y)|Dt1/α(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/2Rdpts0(x,z)Ψs(z,y)dzds+0t/2Rdps0(x,z)Ψts(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)|Dt1/α(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 ), 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{StL}dt)(STK)+, and the price C(T) of the option is defined as C(T)=exp(rT)E[exp(ρ0TI{StL}dt)(STK)+] where r is the risk-free interest rate (see ).

Assume that the price of an asset S={St,t0} 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 JR 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=0n1I{SkT/nL}dt)(STK)+]. 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 condition X and assume that there exists λ>1 1$]]> such that G(λ):=Eexp(λXT)=E(ST)λ<+ . Then |Cn(T)C(T)|exp(rT)ρG(λ)1/λCT,λ/(λ1)(DT,β(n))1/2 with constants CT,λ/(λ1) and DT,β(n) given by (4). The proof is a simple corollary of Theorem 2.1. Denote h(x)=ρIxlogL . Then, keeping the notation of Section 2, we get C(T)=erTEeIT(h)(STK)+,Cn(T)=erTEeIT,n(h)(STK)+. By the Hölder inequality with p=λ and q=λ/(λ1) , |Cn(T)C(T)|erT(E(STK)+λ)1/λ×(E|eIT(h)eIT,n(h)|λ/(λ1))(λ1)/λ. Since |eaeb||ab| for positive a and b, we have |Cn(T)C(T)|erTG(λ)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 λ>2 2$]]>.

Under the assumptions of Proposition 4.1, we have |Cn(T)C(T)|2β2+1max{BT,XρT2(1+ρT)eρT,G(λ)}erTD˜T,β(n), where D˜T,β(n)=n(11/λ)logn,β=1,max(1,T1ββ1)n1/β(11/λ),β>1. 1.\end{array}\right.\]]]>

For N>0 0\$]]>, denote CN(T)=exp(rT)E[exp(ρ0TI{StL}dt)((STK)+N)],CnN(T)=exp(rT)E[exp(ρT/n k=0n1I{SkT/nL}dt)((STK)+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(ρTrT)DT,β(n). For the last two terms, we get |C(T)CN(T)|+|Cn(T)CnN(T)|2exp(rT)E[(STK)+(STK)+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.

References Bogdan, K., Grzywny, T., Ryznar, M.: Density and tails of unimodal convolution semigroups. J. Funct. Anal. 266(6), 35433571 (2014). MR3165234. doi:10.1016/j.jfa.2014.01.007 Ganychenko, I., Kulik, A.: Rates of approximation of nonsmooth integral-type functionals of Markov processes. Mod. Stoch., Theory Appl. 2, 117126 (2014). MR3316480. doi:10.15559/vmsta-2014.12 Ganychenko, I., Kulik, A.: Weak approximation rates for integral functionals of Markov processes. Mod. Stoch., Theory Appl. 3, 251266 (2015) Gobet, E., Labart, C.: Sharp estimates for the convergence of the density of the Euler scheme in small time. Electron. Commun. Probab. 13, 352363 (2008). MR2415143. doi:10.1214/ECP.v13-1393 Guerin, 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.3130 Knopova, V.: Compound kernel estimates for the transition probability density of a Lévy process in Rn . Theory Probab. Math. Stat. 89, 5770 (2014). MR3235175. doi:10.1090/s0094-9000-2015-00935-2 Knopova, V., Kulik, A.: Parametrix construction of the transition probability density of the solution to an SDE driven by α-stable noise. arXiv:1412.8732 Knopova, V., Kulik, A.: Intrinsic small time estimates for distribution densities of Lévy processes. Random Oper. Stoch. Equ. 21(4), 321344 (2013). MR3139314. doi:10.1515/rose-2013-0015 Kohatsu-Higa, A., Makhlouf, A., Ngo, H.L.: Approximations of non-smooth integral type functionals of one dimensional diffusion precesses. Stoch. Process. Appl. 124, 18811909 (2014). MR3170228. doi:10.1016/j.spa.2014.01.003 Kulczycki, T., Ryznar, M.: Gradient estimates of harmonic functions and transition densities for Lévy processes. arXiv:1307.7158. MR3413864. doi:10.1090/tran/6333 Linetsky, V.: Step options. Math. Finance 9, 5596 (1999). MR1849356. doi:10.1111/1467-9965.00063 Mimica, A.: Heat kernel upper estimates for symmetric jump processes with small jumps of high intensity. Potential Anal. 36(2), 203222 (2012). MR2886459. doi:10.1007/s11118-011-9225-1 Pruitt, W.E., Taylor, S.J.: The potential kernel and hitting probabilities for the general stable process in Rn . Trans. Am. Math. Soc. 146, 299321 (1969). MR0250372 Rosinski, J., Singlair, J.: Generalized tempered stable processes. Stab. Probab. 90, 153170 (2010). MR2798857. doi:10.4064/bc90-0-10 Sztonyk, P.: Estimates of tempered stable densities. J. Theor. Probab. 23(1), 127147 (2010). MR2591907. doi:10.1007/s10959-009-0208-8 Sztonyk, P.: Transition density estimates for jump Lévy processes. Stoch. Process. Appl. 121, 12451265 (2011). MR2794975. doi:10.1016/j.spa.2011.03.002 Watanabe, T.: Asymptotic estimates of multi-dimensional stable densities and their applications. Trans. Am. Math. Soc. 359(6), 28512879 (2007). MR2286060. doi:10.1090/S0002-9947-07-04152-9