Accuracy of discrete approximation for integral functionals of Markov processes

The article is devoted to the estimation of the rate of convergence of integral functionals of a Markov process. Under the assumption that the given Markov process admits a transition probability density which is 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. Some examples are provided.


Introduction
Let X t , t ≥ 0, be a Markov process with values in R d . Consider an integral functional of the form where h : R d → R is a given measurable function. In this paper we investigate the accuracy of the approximation of I T (h) by the Riemannian sums h(X (kT )/n ), n ≥ 1.
The function h is assumed to be bounded, only; i.e., we do not impose any regularity assumptions on h. In particular, under this assumption the class of integral functionals which we investigate contains the class of occupation time type functionals (for which h = 1 A for a fixed A ∈ B(R d )), 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 the approximation methods. As an example of such a situation, we mention the so-called occupation time option [14], 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 when the underlying process is a Lévy process which is "spectrally negative" (i.e. does not have positive jumps, see [5]), and practically more realistic cases of general Lévy processes, solutions to Lévy driven SDE's, etc. can be treated only numerically. To estimate the rate of convergence of the respective Monte-Carlo approximative methods, one needs to estimate the accuracy of various approximation steps involved in the algorithm. In this paper we focus on solving of such a problem for the discrete approximation of the integral functional of type (1).
For diffusion processes, this problem was studied in [4] and recently in [11], by means of the methods involving the particular structural features of the process, e.g. the Malliavin calculus tools. On the other hand, in two recent papers [3], [13], 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.
X. The transition probability P t (x, dy) of X admits a density p t (x, y) w.r.t. the Lebesgue measure on R d . This density is differentiable w.r.t. t, and its derivative possesses the bound for some B T,X ≥ 1, β ≥ 1, and measurable function q, such that for each fixed t, x the function q t,x (·) is a distribution density. In [3], [13], the condition similar to X was formulated with β = 1. Such a condition is verified for the particularly important classes of diffusion process and symmetric α-stable processes. However, in some natural cases one can expect to get (2) only with β > 1. As the simplest and the most illustrative example one can take an α-stable process with drift: where c = 0 and Z is a (e.g. symmetric) α-stable process. Then where g (α) denotes the distribution density of Z 1 . Straightforward calculation shows that (2) holds true with β = max(1, 1/α), which is strictly greater than 1 when α < 1. Since the Lévy noises are now used extensively in various applied models, the simple calculation made above shows that it is highly desirable to extend the results of [3] and [13], 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 the tools which 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 [10], with the process X being a solution of the SDE 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 R d . Recently, apart from the stable noises, various types of "locally stable" noises are frequently used in applied models: tempered stable processes, damped stable processes, etc. Heuristically, for a "locally stable" process its "small jumps behavior" is the same as for the stable one, but the "large jumps behavior" of the former is drastically different from "tail behavior" of the Lévy measure. Since (2) is genuinely related to "local behavior" of the process, one can expect that the results of [10] should have a natural extension to the case of "locally stable" Z. However, to make such a conjecture rigorous is a sophisticated problem; 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 p t (x, y), and hence any bound for p t (x, y), which one may expect to design within this method, is at least as complicated as 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], [9], [12], [15], [16], [17], [18], [19], [21]; 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 the exact two-sided bounds for p t (x, y) one should impose quite non-trivial structural assumptions on the "tails" of the Lévy measure even in a comparatively simple "locally stable" case. Motivated by this observation on one hand, and by the initial approximation problem which suggests the condition (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 one can get the bound (2), where at the left hand side we have the transition probability density of the SDE driven by a "locally stable" process, and at the right hand side we have a (properly shifted) 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" e.g. for exponentially tempered α-stable law. The gain is, however, that under a mild set of assumptions we obtain a uniform-in-class upper bound, which is clearly easy to use 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), see below. 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 two our 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 will prove two results. The first one concerns the "strong approximation rate", i.e. the control on the L p -distance between I T (h) and its approximation I T,n (h).
The second result concerns the "weak approximation", i.e. the control on the difference between the expectations of certain terms, which involve I T (h) together with its approximation I T,n (h).
Theorem 2.2. Suppose that X holds. Then for any k ∈ N and any bounded function f we have Remark 2.2. This theorem extends [13, Theorem 1.1], where it was assumed that β = 1. In the proof below, we concentrate on the case β > 1.
Using the Taylor expansion, one can obtain directly the following corollary of Theorem 2.2.
Corollary 2.1. Suppose that X holds, and let ϕ be an analytic function defined in some neighbourhood of 0. In addition, suppose that the constants D ϕ , R ϕ > 0 are such that Then for any bounded function f and a function h such that T h < R ϕ , we have the following bound: 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 to make the exposition more transparent.
Proof. Let us introduce the notation used throughout the whole section: for t ∈ [kT /n, (k + 1)T /n), k ≥ 0, we put η n (t) = kT n , ζ n (t) = (k+1)T n ; that is, η n (t) is the point of the partition {T k/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, which is strictly larger than t.
We have for some 1 ≤ k n,β ≤ n, which will be chosen later. We estimate each term separately.
For M 1 we have Further, using (2), we get Now we finalize the argument.
1) If β = 1, put k n,β = 1, n ≥ 1. Then we get 2) If β > 1, put k n,β = [n 1−1/β ] + 1. Then Therefore, Proof of Theorem 2.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 By definition, the function t → J t,n (h) is absolutely continuous. Then using the Newton-Leibnitz formula twice we get Therefore, Let us estimate separately the expectations of H 1 T,n,p (h) and H 2 T,n,p (h). By the Hölder inequality, Further, observe that for every s the variables are F ζn(s) -measurable; here and below {F t , t ≥ 0} denotes the natural filtration for the process X. Hence, By Proposition 2.1 and the Markov property of X, we have Therefore, using the Hölder inequality, we get Note that n −1 ≤ D T,β (n), hence the above bounds for E x H 1 T,n,p (h) and E x H 2 T,n,p (h) finally yield the estimate It can be easily seen that the above inequality also holds true if J T,n (h) in the left hand side is replaced by J t,n (h). Taking the integral over t ∈ [0, T ], we get Because h is bounded, the left hand side expression in the above inequality is finite. Hence, resolving this inequality, we get which together with (8) gives the required statement.
Let us estimate the r-th term in the last line in (9). We have Since the case β = 1 was already treated in [13], for the rest of the proof we assume that β > 1.
Consider two cases: a) s r − s r−1 > k n,β T /n and b) s r − s r−1 ≤ k n,β T /n. In case a), using condition X and the Chapman-Kolmogorov equation, we derive v −β q v,y r−1 (y r )dv dy r−1 dy r .
Since k n,β ≥ 2, then in case a) we have s r − s r−1 ≥ 2T /n, and hence Therefore, using the fact that q t,y (·) is the probability density for any t > 0 and y ∈ R d , we finally get In case b) we simply apply the Chapman-Kolmogorov equation: p ηn(s r−1 )−ηn(s 0 ) (x, y r−1 ) p sr−s r−1 (y r−1 , y r ) + p ηn(sr)−ηn(s r−1 ) (y r−1 , y r ) dy r−1 dy r ≤ 2.
Therefore, summarizing the estimates obtained in cases a) and b) we get, using (7), the estimates where in the fourth and the fifth lines we used that s r−1 r−1 ≤ s r−1 r . Taking into account that in (9) we have k terms, and the common multiplier k!, we finally arrive at (5).

Condition X for solutions to Lévy driven SDEs
Consider the SDE where Z is a real-valued Lévy process. In [10] it was shown that if Z t is a symmetric α-stable process and b(·) is bounded and Lipschitz continuous, then the solution to equation (10) satisfies condition X with β = max(1, 1/α) (in fact, in [10] more general multidimensional SDEs of the form (3) are considered). In this section we outline the argument which makes it possible to extend the class of "Lévy noises". Namely, we will omit the requirement on Z to be symmetric, and relax the stability assumption, demanding Z to be "locally α-stable" in the sense we specify below.
Recall that for a real-valued Lévy process the characteristic function is of the form where the characteristic exponent ψ admits the Lévy-Khinchin representation In what follows, we assume that σ 2 = 0 and the Lévy measure µ is of the form with some C ± ≥ 0 and m(u) ≥ 0 such that m(u) = 0 for |u| ≤ 1, and 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" the α-stable Lévy measure. Let us impose three minor conventions, which will simplify the technicalities below. First, since we are mostly interested in the case β > 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 to this integral; that is, ψ has the form Clearly, this does not restrict the generality because one can change the constant a by changing respectively the drift coefficient b(·) in (10). Finally, in order to avoid the usage of the Rademacher theorem (see [10,Lemma 7.4] for the case when b is just Lipschitz continuous), let us assume that b ∈ C 1 (R).
In what follows we show how the parametrix construction developed in [10] can be modified to provide the representation and the bonds for the transition probability density p t (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 Generically, the parametrix construction provides the representation of the required transition probability density in the form (14) p Here p 0 t (x, y) is a "zero approximation term" for the unknown p t (x, y), the function Ψ t (x, y) is given by the "convolution series" the function Φ t (x, y) depends on the particular choice of p 0 t (x, y), and equals (16) Φ t (x, y) is the formal generator of the process X. The subscript x in above expressions means that the operator is applied with respect to the variable x. Note that to make the above construction feasible, one should properly choose the "zero approximation term" p 0 t (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 p 0 t (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, and the characteristic exponent Note that since the process Z (α,C ± ) possesses the scaling property . By the scaling property we have Denote also by Z (α) the symmetric α-stable process; that is, the process of the form introduced above with C + = C − = 1. Let g (α) t be the respective distribution density and g (α) := g (α) 1 . Finally, denote by χ t (x) and θ t (y), respectively, the solutions to the ODEs Note that these solutions exist, because b(·) is Lipschitz continuous. Now we are ready to formulate the main statement of this section.
Then the convolution series (15) is well defined, and the formula (14) gives the representation of the transition probability density p t (x, y) of the process X. This density and its time derivative have the following upper bounds: Consequently, the process X satisfies condition X with β = 1/α and q t,x (y) = g (α) Proof. First we evaluate Φ t (x, y). If it is not stated otherwise, we assume in all estimates obtained below that t ∈ (0, T ] for some T > 0, and x, y ∈ R. Observe that g (α,C ± ) ∈ C 2 b (R). Indeed, this property easily follows from the Fourier inversion formula and the expression for the characteristic Let us approximate the price C(T ) of our option by C n (T ) = exp(−rT )E exp −ρT /n n−1 k=0 I {S kT /n ≤L} dt (S T − 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.
Proposition 4.1. Suppose that the process X satisfies condition X, and assume that there exists with constants C T,λ/(λ−1) and D T,β (n) are given by (4).
Proof. The proof is a simple corollary of Theorem 2.1. Denote h(x) = ρI x≤log L . Then keeping the notation of Section 2 we get We also can control the accuracy of the approximation using the weak rate bound from Theorem 2.2. Observe that the bound given below is sharper than those obtained in the previous proposition precisely when λ > 2. To complete the proof, put N = n 1/(βλ) .