Let $X_{t}$, $t\ge 0$, be a Markov process with values in ${\mathbb{R}}^{d}$. Consider an integral functional of the form where $h:{\mathbb{R}}^{d}\to \mathbb{R}$ is a measurable function. In this paper, we investigate the accuracy of the approximation of $I_{T}(h)$ by the Riemannian sums 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=1_{A}$ for a fixed $A\in \mathcal{B}({\mathbb{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 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.

In [2, 3], the condition similar to X was formulated with $\beta =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 $\beta >1$. As the simplest and the most illustrative example, we can take an α-stable process with drift: where $c\ne 0$, and Z is an (e.g., symmetric) α-stable process. Then where ${g}^{(\alpha )}$ denotes the distribution density of $Z_{1}$. Straightforward calculation shows that (2) holds with $\beta =\max (1,1/\alpha )$, which is strictly greater than 1 when $\alpha <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 $\beta =1$, to the more general case of X with arbitrary $\beta \ge 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 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 ${\mathbb{R}}^{d}$. 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 $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 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 $p_{t}(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.

In this section, we prove two results. The first one concerns the “strong approximation rate”, that is, 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,” that is, the control on the difference between the expectations of certain terms, which involve $I_{T}(h)$ together with its approximation $I_{T,n}(h)$.

Using the Taylor expansion, we can directly obtain the following corollary of Theorem 2.2.

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\equiv 1$, but we give it separately in order to make the exposition more transparent.

Consider the SDE where Z is a real-valued Lévy process. In [7], it was shown that if $Z_{t}$ is a symmetric α-stable process and $b(\cdot )$ is bounded and Lipschitz continuous, then the solution to Eq. (10) satisfies condition X with $\beta =\max (1,1/\alpha )$ (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 where the characteristic exponent ψ admits the Lévy–Khinchin representation In what follows, we assume that ${\sigma }^{2}=0$ and the Lévy measure μ is of the form with some $C_{\pm }\ge 0$ and $m(u)\ge 0$ such that $m(u)=0$ for $|u|\le 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” 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 $\beta >1$, we assume that $\alpha <1$. Second, the latter assumption assures that the integral is well defined, and we assume that the constant a in (11) equals this integral; that is, ψ has the form Clearly, this does not restrict the generality because we can change the constant a by changing appropriately the drift coefficient $b(\cdot )$ 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\in {C}^{1}(\mathbb{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 $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 a representation of the required transition probability density in the form Here ${p_{t}^{0}}(x,y)$ is a “zero approximation term” for the unknown $p_{t}(x,y)$, and the function $\varPsi _{t}(x,y)$ is given by the “convolution series” where the function $\varPhi _{t}(x,y)$ depends on the particular choice of ${p_{t}^{0}}(x,y)$ and equals where 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” ${p_{t}^{0}}(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_{t}^{0}}(x,y)$ and then to construct the bounds for the associated $\varPhi _{t}(x,y)$ and its convolution powers, we need some more notation.

Denote by ${Z}^{(\alpha ,C_{\pm })}$ the α-stable process with the Lévy measure $\mu _{\alpha ,C_{\pm }}(du)={m}^{(\alpha ,C_{\pm })}(u)\hspace{0.1667em}du$, and the characteristic exponent Note that since the process ${Z}^{(\alpha ,C_{\pm })}$ possesses the scaling property Denote by ${g_{t}^{(\alpha ,C_{\pm })}}$ the distribution density of ${Z_{t}^{(\alpha ,C_{\pm })}}$. By the scaling property we have Denote also by ${Z}^{(\alpha )}$ the symmetric α-stable process, that is, the process of the form introduced before with $C_{+}=C_{-}=1$. Let ${g_{t}^{(\alpha )}}$ be the respective distribution density, and ${g}^{(\alpha )}:={g_{1}^{(\alpha )}}$.

Finally, denote by $\chi _{t}(x)$ and $\theta _{t}(y)$, respectively, the solutions to the ODEs Note that these solutions exist because $b(\cdot )$ is Lipschitz continuous.

Now we are ready to formulate the main statement of this section.

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 and the price $\mathbf{C}(T)$ of the option is defined as where r is the risk-free interest rate (see [11]).

Assume that the price of an asset $S=\{S_{t},t\ge 0\}$ is of the form where X is the Markov process studied in the previous sections. Then the time spent by the process S in a set $J\subset \mathbb{R}$ equals the time spent by X in the set ${J^{\prime }}=\log J$.

Let us approximate the price $\textbf{C}(T)$ of our option by 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.

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 $\lambda >2$.