Exact Simulation for Multivariate It\^o Diffusions

TL;DR

This paper introduces the first generic exact simulation algorithm for multivariate Itô diffusions, overcoming the limitations of Lamperti transformations which are only applicable in one dimension.

Contribution

It develops a novel exact sampling method for multivariate diffusions by combining rough path theory and multilevel Monte Carlo techniques.

Findings

First generic exact simulation algorithm for multivariate diffusions

Successfully overcomes Lamperti transformation limitations in higher dimensions

Demonstrates practical applicability of combined rough path and Monte Carlo methods

Abstract

We provide the first generic exact simulation algorithm for multivariate diffusions. Current exact sampling algorithms for diffusions require the existence of a transformation which can be used to reduce the sampling problem to the case of a constant diffusion matrix and a drift which is the gradient of some function. Such transformation, called Lamperti transformation, can be applied in general only in one dimension. So, completely different ideas are required for exact sampling of generic multivariate diffusions. The development of these ideas is the main contribution of this paper. Our strategy combines techniques borrowed from the theory of rough paths, on one hand, and multilevel Monte Carlo on the other.

Equations

$$dX(t)=\mu (X(t))dt+\sigma (X(t))dW(t),X(0)=x_0,$$

$$dX(t)=\mu (X(t))dt+dW(t),X(0)=x_0,$$

$$X(t)=(X_1(t),\cdots ,X_d(t));0\le t\le 1.$$

$$L(t)=\exp \left({\int }_0^t{\mu }^T(X(t))dX(t)-\frac{1}{2}{\int }_0^t\Vert \mu (X(t)){\Vert }_2^{2}dt\right),$$

$$\mathbb{P}(A)={\mathbb{E}}^{\tilde{\mathbb{P}}}\left[I(A)L(1)\right];\forall A\in \mathcal{F},$$

$$W(t)=(W_1(t),\cdots ,W_d(t));0\le t\le 1$$

$$W(t)=X(t)-{\int }_0^t\mu (X(s))ds;0\le t\le 1.$$

$$dY(t)=\alpha (Y(t))dt+\nu (Y(t))dW(t),Y(t)=y_0$$

$$\underset{t\in [0,1]}{\sup }\Vert Y_{\epsilon }(t)-Y(t){\Vert }_2\le \epsilon$$

$$\mathbb{P}\left(\{\omega \in \Omega :Y(1)\in {\mathring{G}}_i\}\right)=\mathbb{P}\left(\bigcup _{\epsilon >0}\{\omega \in \Omega :d(Y_{\epsilon }(1),G_i^c)>ϵ\}\right).$$

$$d\Big{(}Y_{\varepsilon}(1),G_{\Xi_{\mathcal{G}}(Y_{\varepsilon}(1))}^{c}\Big{)}>\varepsilon,$$

$$\mathbb{P}\left(\bigcup _{i\in \mathbb{N}}\bigcup _{\epsilon >0}\{\omega \in \Omega :d(Y_{\epsilon }(1),G_i^c)>ϵ\}\right)=1.$$

$$\{\begin{array}{l}dL(t)=L(t)\Vert \mu (X(t)){\Vert }_2^{2}dt+L(t){\mu }^T(X(t))dW(t),\\ dX(t)=\mu (X(t))dt+dW(t),\end{array}$$

$$\mathbb{P}\Big{(}X(1)\in dx\Big{|}\Xi_{\mathcal{G}}\big{(}(L(1),X(1))\big{)}=i\Big{)}=\tilde{\mathbb{P}}\Big{(}X(1)\in dx\Big{|}\max\big{(}i,(i+1)U\big{)}<L(1)<i+1\Big{)}$$

$$\mathbb{P}\Big{(}X(1)\in dx\Big{|}\Xi_{\mathcal{G}}\big{(}(L(1),X(1))\big{)}=i\Big{)}=\frac{\mathbb{P}\Big{(}X(1)\in dx;\Xi_{\mathcal{G}}\big{(}(L(1),X(1))\big{)}=i\Big{)}}{\mathbb{P}\Big{(}\Xi_{\mathcal{G}}\big{(}(L(1),X(1))\big{)}=i\Big{)}}.$$

$$\mathbb{P}\Big{(}X(1)\in dx\Big{|}\Xi_{\mathcal{G}}\big{(}(L(1),X(1))\big{)}=i\Big{)}=\frac{\mathbb{E}^{\tilde{\mathbb{P}}}\Big{[}L(1)I(X(1)\in dx;\Xi_{\mathcal{G}}\big{(}(L(1),X(1)\big{)}=i\Big{]}}{\mathbb{P}\Big{(}\Xi_{\mathcal{G}}\big{(}(L(1),X(1))\big{)}=i\Big{)}}.$$

$$i\le L(1)\le i+1,$$

$$\displaystyle\mathbb{E}^{\tilde{\mathbb{P}}}\Big{[}L(1)I(X(1)\in dx;\Xi_{\mathcal{G}}\big{(}(L(1),X(1))\big{)}=i)\Big{]}$$

$$=$$

$$\displaystyle\mathbb{P}\Big{(}X(1)\in dx\Big{|}\Xi_{\mathcal{G}}$$

$$\displaystyle\times\tilde{\mathbb{P}}\Big{(}X(1)\in dx\Big{|}\max\big{(}i,(i+1)U\big{)}<L(1)<i+1\Big{)}.$$

$$(i+1)\tilde{\mathbb{P}}\Big{(}\max\big{(}i,(i+1)U\big{)}<L(1)<i+1\Big{)}=\mathbb{P}\Big{(}\Xi_{\mathcal{G}}\big{(}(L(1),X(1))\big{)}=i\Big{)}.$$

$$\displaystyle\mathbb{P}\Big{(}\Xi_{\mathcal{G}}\big{(}(L(1),X(1)\big{)}=i\Big{)}$$

$$=$$

$$dL(t)=L(t){\mu }^T(X(t))dX(t)$$

$$|p(x)-p(y)|\le C_S|x-y|\forall x,y\in S.$$

$$p(x)\ge {\delta }_S\forall x\in S.$$

$${\mathcal{G}}_{\text{loc}}=\{[i_1,i_1+1)\times \cdots \times [i_d,i_d+1):(i_1,\dots ,i_d)\in {\mathbb{Z}}^d\}.$$

$$\overline{p_n}(x)=[V(r_n){]}^{-1}{\int }_{B_{r_n}(x)}p(x)dx.$$

$${\hat{p}}_n(x)=[V(r_n){]}^{-1}\times I(X(1)\in B_{r_n}(x))$$

Full text

\authornames

J. BLANCHET AND F. ZHANG

\authorone

[Stanford University]Jose Blanchet \authortwo[Stanford University]Fan Zhang \addressoneHuang Engineering Center, 475 Via Ortega, Stanford, CA 94305, United States. \emailonejose.blanchet@stanford.edu \addresstwoHuang Engineering Center, 475 Via Ortega, Stanford, CA 94305, United States. \emailtwofzh@stanford.edu

\ams

34K50,65C05,82B8097K60

1 Introduction

Consider a probability space $(\Omega,\mathcal{F},\mathbb{P})$ and an Itô Stochastic Differential Equation (SDE)

[TABLE]

where $W(\cdot)$ is a $d^{\prime}$-dimensional Brownian motion under $\mathbb{P}$, and $\mu(\cdot)=(\mu_{i}(\cdot))_{d}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ and $\sigma(\cdot)=(\sigma_{ij}(\cdot))_{d\times d^{\prime}}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d\times d^{\prime}}$ satisfy suitable regularity conditions. For instance, in order for (1) to have a strong solution, it is sufficient to assume that both $\mu\left(\cdot\right)$ and $\sigma\left(\cdot\right)$ are uniformly Lipschitz.

Under additional regularity conditions on $\mu\left(\cdot\right)$ and $\sigma\left(\cdot\right)$, this paper provides the first Monte Carlo simulation algorithm which allows sampling any discrete skeleton $X\left(t_{1}\right),...,X\left(t_{m}\right)$ exactly, without any bias.

The precise regularity conditions that we impose on $\mu\left(\cdot\right)$ and $\sigma\left(\cdot\right)$ are stated in Section 3. In particular, it is sufficient for the validity of our Monte Carlo method to assume $\mu\left(\cdot\right)$ and $\sigma\left(\cdot\right)$ to be three times continuously differentiable, both with Lipschitz continuous derivatives of order three. In addition, we must assume that $\sigma\left(\cdot\right)$ is uniformly elliptic.

Exact simulation of SDEs has generated a substantial amount of interest in the applied probability and Monte Carlo simulation communities. The landmark paper of [5], introduced what has become the standard procedure for the design of generic exact simulation algorithms for diffusions. The authors in [5] propose a clever acceptance-rejection sampler which uses Brownian motion as a proposal distribution. The authors in [7] apply a localization technique which eliminates certain boundedness assumptions which are originally present in [5]; see also [4] for the use of retrospective simulation ideas to dispense with boundedness assumptions.

The fundamental assumption underlying the work of [5] and its extensions is that the underlying (target) process has a constant diffusion coefficient, i.e., $\sigma\left(x\right)=\sigma$ for every $x$. Beskos and Roberts [5] note that in the case $d=1$, owing to Lamperti’s transformation, the constant diffusion coefficient assumption comes basically at no cost in generality.

Unfortunately, however, Lamperti’s transformation is only generally applicable in one dimension. In fact, [1] characterizes the multidimensional diffusions for which Lamperti’s transformation can be successfully applied and these models are very restrictive.

Moreover, even if Lamperti’s transformation is applicable in a multidimensional setting, another implicit assumption in the application of the Beskos and Roberts acceptance-rejection procedure is that the drift coefficient $\mu\left(\cdot\right)$ is the gradient of some function (i.e. $\mu\left(x\right)=\nabla v\left(x\right)$ for some $v\left(\cdot\right)$). This assumption, once again, comes at virtually no cost in generality in the one-dimensional setting, but it may be very restrictive in the multidimensional case.

Because of these limitations, a generic algorithm for exact simulation of multidimensional diffusions, even under the regularity conditions that we impose here, requires a completely different set of ideas.

The contribution in this paper is therefore not only the production of such a generic exact simulation algorithm, but also the development of the ideas that are behind its construction. In Section 2 we introduce an algorithm which assumes a constant diffusion coefficient, but removes the assumption of the drift coefficient being the gradient of some function. In Section 3, we eventually remove the requirement on a constant diffusion matrix and propose an algorithm applicable to general diffusions. The algorithms in Section 2 and Section 3 are different in nature. However, they share some common elements, such as the use of so-called Tolerance-Enforced Simulation techniques based on rough path estimates. Even though the algorithm in Section 3 is more general, we believe that there is significant value in developing the algorithm in Section 2 because of two reasons. The first one is pedagogical, the algorithm in Section 2 is easier to understand while building on a key idea, which involves localizing essential quantities within specific compact domains with probability one. The second reason is that the algorithm in Section 2, being simpler, may be subject to potential improvement methodologies to be pursued in future research.

Potential improvements are particularly interesting directions specially given that, unfortunately, the algorithms that we present have infinite expected termination time. We recognize that this issue should be resolved for the algorithms to be widely used in practice, and we discuss the elements which lead to infinite expected running time in Section 4. There are basically two types of elements that affect the running time of the algorithm in Section 3, one of them has to do with the types of issues that arise when trying to fully remove the bias in Tolerance-Enforced Simulation and related approximations, and the other issue has to do with the use of a density approximation coupled with Bernoulli factories. In contrast, the algorithm in Section 2 is only affected by removing the bias in Tolerance-Enforced Simulation type approximations. We must stress, however, that the present paper shows for the first time that it is possible to perform exact sampling of multidimensional diffusions in substantial generality and, in doing so, it provides a conceptual framework different from the prevailing use of Lamperti transformation, which is the only available generic approach for producing exact sampling of diffusions.

Now, despite the algorithm’s practical limitations, it is vital to recognize the advantages that unbiased samplers have over biased samplers in the context of a massive parallel computing environment, because it is straightforward to implement a parallel procedure to reduce the estimation error of an unbiased sampler.

Recently, there have been several unbiased estimation procedures which have been proposed for expectations of the form $\alpha=\mathbb{E}\left(f\left(X\left(t\right)\right)\right)$, assuming $\text{Var}\left(f\left(X\left(t\right)\right)\right)<\infty$. For example, the work of [21] shows that if $f\left(\cdot\right)$ is twice continuously differentiable (with Lipschitz derivatives) and if there exists a discretization scheme which can be implemented with a strong convergence error of order 1, then it is possible to construct an unbiased estimator for $\alpha$ with finite variance and finite expected termination time. The work of [12] shows that such a discretization scheme can be developed if $\mu\left(\cdot\right)$ and $\sigma\left(\cdot\right)$ are sufficiently smooth under certain boundedness assumptions. The paper [13] also develops an unbiased estimator for $\alpha$ using a regime-switching technique. Our work here is somewhat related to this line of research, but an important difference is that we actually generate $X\left(t\right)$ exactly, while all of the existing algorithms which apply in multidimensional diffusion settings generate $Z$ such that $\mathbb{E}\left(Z\right)=\alpha$. So, for example, if $f\left(\cdot\right)$ is positive, one cannot guarantee that $Z$ is positive using the type of samplers suggested in [21]. However, by sampling $X\left(t\right)$ directly, one maintains the positivity of the estimator.

Another instance in which direct exact samplers are useful arises in the context of stochastic optimization. For instance, consider the case in which one is interested in optimizing a convex function of the form $g\left(\theta\right)=\mathbb{E}\left(h\left(X\left(t\right),\theta\right)\right)$, where $h\left(x,\cdot\right)$ is differentiable. In this case, one can naturally construct an estimator $Z\left(\theta\right)$ such that $g\left(\theta\right)=\mathbb{E}\left(Z\left(\theta\right)\right)$ using the results in [21] and optimize the mapping $\theta\rightarrow$ $n^{-1}\sum_{i=1}^{n}Z_{i}\left(\theta\right)$, which unfortunately will typically not be convex. So, having access to a direct procedure to sample $X\left(t\right)$ in this setting is particularly convenient as convexity is preserved.

The rest of the paper is organized as follows. In Section 2, we consider the case of multidimensional diffusions with a constant diffusion coefficient and a Lipschitz continuous (suitably smooth) drift. The general case is discussed in Section 3, our development uses localization ideas which are introduced in Section 2, but also some basic estimates of the transition density of the underlying diffusion (e.g. Lipschitz continuity), these estimates are developed in Appendix A. As mentioned before we discuss the bottlenecks in the expected running time of the algorithm in Section 4.

2 Exact Simulation of SDEs with Identity Diffusion Coefficient

In case that the Lamperti’s transformation is applicable, the SDE of interest is reducible to another SDE whose diffusion matrix is the identity. As a result, it suffices to consider simulating the following SDE

[TABLE]

where $W=\left\{W(t)=(W_{1}(t),\cdots,W_{d}(t)):0\leq t<\infty\right\}$ is a $d$-dimensional Brownian motion. In this section we concentrate on the identity diffusion case (2), but the development can be immediately extended to the case of a constant diffusion matrix. However, throughout this section we must impose the following assumptions.

{assumption}

The SDE (2) has a strong solution.

{assumption}

The drift coefficient $\mu(\cdot)$ is three times continuously differentiable.

Assumption 2 is the requirement of TES, the theoretical foundation of our algorithm, which we shall introduce later.

Let us introduce some notations first. For any set $G$ and $x\in\mathbb{R}^{d}$, we use $d(x,G)=\inf\left\{\|x-y\|_{2}:y\in G\right\}$ to denote the distance between $x$ and $G$; $\mathring{G}$ denotes the interior of $G$; $\partial G$ denotes the boundary of $G$; $G^{c}$ denotes complementary of $G$.

Consider a probability space $(\Omega,\mathcal{F},\tilde{\mathbb{P}})$ endowed with a filtration $\left\{\mathcal{F}_{t}:0\leq t\leq 1\right\}$, and supporting a $d$-dimensional Brownian motion

[TABLE]

Let $\left\{L(t):0\leq t\leq 1\right\}$ to be a $\tilde{\mathbb{P}}$-local martingale defined as

[TABLE]

where $\mu^{T}(\cdot)$ denotes the transpose of the column vector $\mu(\cdot)$. Under Assumption 2, $L(\cdot)$ is a $\tilde{\mathbb{P}}$-martingale, see Corollary 3.5.16 of [15].

In this case we can define a probability measure $\mathbb{P}$ through

[TABLE]

where $I(A)$ denotes the indicator function of the set $A$ and $\mathbb{E}^{\tilde{\mathbb{P}}}\left(\cdot\right)$ is the expectation operator under $\tilde{\mathbb{P}}$.

Let

[TABLE]

be a $d$-dimensional process defined by

[TABLE]

The following theorem provides the distribution of $W(\cdot)$.

Theorem 2.1** (Girsanov Theorem).**

If Assumption 2 is satisfied, then the process $W(\cdot)$ is a $d$-dimensional Brownian motion on probability space $(\Omega,\mathcal{F},\mathbb{P}).$

Proof 2.2**.**

See, for instance, Theorem 3.5.1 of [15].

It is readily apparent from (4) that $X(\cdot)$ is a weak solution to SDE (2) on the probability space $(\Omega,\mathcal{F},\mathbb{P})$. The exact simulation problem becomes sampling $X(1)$ under measure $\mathbb{P}$. Since $X(1)$ follows a normal distribution under measure $\tilde{\mathbb{P}}$, we can attempt to use acceptance-rejection to sample $X(1)$. A direct application of acceptance-rejection may proceed by using the $\tilde{\mathbb{P}}$ distribution of $X(1)$ (which is simply normal distribution) as the proposal, which, if acceptance-rejection is applicable, would then be accepted with probability proportional to $L(1)$. However, there are two obstacles when trying to apply such a direct acceptance-rejection approach. First, the presence of the general stochastic integral appearing in the definition of $L\left(1\right)$ makes the likelihood ratio difficult to directly compute. Second, a direct application of acceptance-rejection requires the likelihood ratio, $L(1)$, to be bounded, which is unfortunately violated.

In order to deal with the first obstacle, we note that it is really not necessary to accurately evaluate the likelihood ratio. In the standard procedure of acceptance-rejection, the likelihood ratio is only used for comparison with an independent uniform random variable. Thus, to address the first obstacle, we can approximate the likelihood ratio with a deterministic error bound, and keep refining until we can decide to either accept or reject the proposal. It turns out, as we shall see in Corollary 2.7, that the same approximation technique can actually be used to localize $L\left(1\right)$ and also resolve the second obstacle. Then, we will sample the distribution of $X(1)$ conditional on the localization of $L\left(1\right)$ using acceptance-rejection, where the rejection scheme is suggested by Lemma 2.1.

The theoretical foundation for such approximation and refinement strategy is given by Tolerance-Enforced Simulation, which is presented in Theorem 2.3.

Theorem 2.3** (Tolerance-Enforced Simulation).**

Consider a probability space $(\Omega,\mathcal{F},\mathbb{P})$ and the following SDE:

[TABLE]

where $\alpha(\cdot)=(\alpha_{i}(\cdot))_{d}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$, $\nu(\cdot)=(\nu_{ij}(\cdot))_{d\times d^{\prime}}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d\times d^{\prime}}$ and $W(\cdot)$ is a $d^{\prime}$-dimensional Brownian motion. Suppose that $\alpha(\cdot)$ is continuously differentiable and that $\nu(\cdot)$ is three times continuously differentiable. Then, given any deterministic $\varepsilon>0$, there is an explicit Monte Carlo procedure that allows us to simulate a piecewise constant process $Y_{\varepsilon}(\cdot)$, such that

[TABLE]

with probability one. Furthermore, for any $m>1$ and $0<\varepsilon_{m}<\dots<\varepsilon_{1}<1$, we can simulate $Y_{\varepsilon_{m}}$ conditional on $Y_{\varepsilon_{1}},\dots,Y_{\varepsilon_{m-1}}$.

Proof 2.4**.**

See Theorem 2.1, Theorem 2.2 and Section 2.1 of [6], where a detailed procedure of Tolerance-Enforced Simulation is also provided.

Remark 2.5**.**

Tolerance-Enforced Simulation is based on the Lévy-Ciesielski Construction of the driving Brownian motion $W(\cdot)$ up to a random level. Consequently, $W(1)$ is obtained for free when we run TES in which a skeleton of the driving Brownian motion $W(\cdot)$ is simulated. In particular, for any $m>1$ and $0<\varepsilon_{m}<\dots<\varepsilon_{1}<1$, we can simulate $Y_{\varepsilon_{m}}$ conditional on $Y_{\varepsilon_{1}},\dots,Y_{\varepsilon_{m-1}}$ and $W(1)$.

As a straightforward consequence of Theorem 2.3, we develop a localization procedure of SDE in Corollary 2.7. Before moving forward to state the result, we define some notations that will be therein used.

Definition 2.6**.**

A family of (Borel measurable) sets $\mathcal{G}=\{G_{i}\subset\mathbb{R}^{d}:i\in\mathbb{N}\}$ is said to be a countable continuous partition for a $d$-dimensional random vector $Y$, if and only if

The sets in $\mathcal{G}$ are mutually disjoint, i.e. $G_{i}\cap G_{j}=\emptyset$ for $i\neq j$; 2. 2.

$Y$* is concentrated on $\mathcal{G}$, that is $\mathbb{P}(Y\in\cup_{i\in\mathbb{N}}G_{i})=1$;* 3. 3.

$\mathbb{P}(Y\in\partial G_{i})=0,\forall i\in\mathbb{N}$.

In addition, a function $\Xi_{\mathcal{G}}(x):\text{supp}(Y)\rightarrow\mathbb{N}$ is defined such that $\Xi_{\mathcal{G}}(x)=i$ if and only if $x\in G_{i}$.

Corollary 2.7**.**

Under the setting of Theorem 2.3, let $\mathcal{G}=\{G_{i}:i\in\mathbb{N}\}$ be a countable continuous partition for $Y(1)$, then there is an algorithm for simulating $\Xi_{\mathcal{G}}(Y(1))$ that terminates in finite time with probability one. In particular, for any set $G$ such that $\mathbb{P}(Y(1)\in\partial G)=0$, there is an algorithm for simulating $I(Y(1)\in G)$.

Proof 2.8**.**

Notice that $\mathbb{P}(Y(1)\in\partial G_{i})=0$, so $Y(1)\in\bigcup_{i\in\mathbb{N}}\mathring{G}_{i}$ holds almost surely. Recalling from Theorem 2.3 that $\|Y_{\varepsilon}(1)-Y(1)\|_{2}\leq\varepsilon$ a.s., which suggests that

[TABLE]

Thus, we pick $\varepsilon\in(0,1)$ and apply TES to simulate the approximation process $Y_{\varepsilon}(1)$. If

[TABLE]

then $\Xi_{\mathcal{G}}(Y(1))=\Xi_{\mathcal{G}}(Y_{\varepsilon}(1))$, which terminates the algorithm. Otherwise we keep refining the approximation of TES, by setting $\varepsilon\leftarrow\varepsilon/2$, until d\Big{(}Y_{\varepsilon}(1),G_{\Xi_{\mathcal{G}}(Y_{\varepsilon}(1))}^{c}\Big{)}>\varepsilon. The algorithm will ultimately terminate since

[TABLE]

The procedure for simulation of $I(Y(1)\in G)$ is just a particular case, by setting $\mathcal{G}=\{G,G^{c}\}$. The details of the algorithm are given in Algorithm 1.

The algorithm for simulating $X(1)$ is performed in a two-stage fashion. At first stage, the likelihood ratio $L(1)$ is localized with the help of Corollary 2.7. (The efficiency of the algorithm may be slightly improved if we localize $X(1)$ and $L(1)$ simultaneously at the first step, then applying acceptance-rejection based on localization. However, this does not solve the problem of the infinite expected running time.) Then, at second stage, $X(1)$ is sampled conditional on the result of localization.

We now illustrate how to localize $L(1)$ in detail. In order to write the dynamics of $Y(1)$ in standard form as in (5), we consider the SDE of $(L(\cdot),X(\cdot))$ under measure $\mathbb{P}$ as follows,

[TABLE]

Let $\mathcal{G}=\{G_{i}=[i,i+1)\times\mathbb{R}^{d}:i\in\mathbb{N}\}$ in the rest of this section. As (3) guarantees that $L(1)$ is non-negative, it follows immediately that $\mathcal{G}$ is a countable continuous partition for $L(1)$. Therefore, Algorithm 1 is directly applicable to sample $\Xi_{\mathcal{G}}((L(1),X(1)))$ using SDE (6). Without loss of generality, we assume $\Xi_{\mathcal{G}}((L(1),X(1)))=i$ in the rest of this section. It remains to sample $X(1)$ conditional on $\Xi_{\mathcal{G}}((L(1),X(1)))=i$ under probability measure $\mathbb{P}$.

The following lemma provides an alternative expression of the conditional distribution of $X(1)$, which facilitates the simulation of $X(1)$ conditional on localization.

Lemma 2.9**.**

Let $U\sim\textrm{Unif}\;(0,1)$ independent of everything else under probability measure $\tilde{\mathbb{P}}$, then we have

[TABLE]

Proof 2.10**.**

Due to the definition of conditional probability,

[TABLE]

Recall that $d\tilde{\mathbb{P}}=L(1)d\mathbb{P}$, we have

[TABLE]

Since on $\Xi_{\mathcal{G}}((L(1),X(1))=i$,

[TABLE]

we can rewrite the expectation into a probability by introducing $U\sim\text{Unif}(0,1)$, namely,

[TABLE]

By substitution, it follows easily that

[TABLE]

It remains to prove that

[TABLE]

By a similar argument we can deduce that

[TABLE]

which ends the proof.

As a direct implication of Lemma 2.9, in order to obtain an example sample for $X\left(1\right)$ under $\mathbb{P}$, given \Xi_{\mathcal{G}}\big{(}(L(1),X(1)\big{)}=i, we can simply simulate $X(1)$ conditional on $\max(i,(i+1)U)<L(1)<i+1$ under probability measure $\tilde{\mathbb{P}}$. In order to do this sampling under $\tilde{\mathbb{P}}$, we can sample $U$ first and denote the value by $u$. Then, observing that $X(\cdot)$ is the driving Brownian motion under the probability measure $\mathbb{P}$, Algorithm 1 is applied to the SDE

[TABLE]

to simulate the indicator function $I(\max(i,u)<L(1)<i+1)$. In addition, according to Remark 2.5, when TES is employed in Algorithm 1, a sample of $X(1)$ is also produced simultaneously. Thereafter, the value of $X(1)$ is accepted if and only if $I(\max(i,u)<L(1)<i+1)=1$; otherwise we repeat the procedure in this paragraph, but we fix the parameter $i$, because $\Xi_{\mathcal{G}}((L(1),X(1)))=i$ has already been sampled under the correct distribution $\mathbb{P}$. The output of the algorithm, once the value $X\left(1\right)$ is finally accepted, follows the distribution of $X(1)$ under $\mathbb{P}$ without any bias.

We summarize the discussion in this section in the following theorem.

Theorem 2.11**.**

If Assumption 2 and 2 are satisfied, then there is an exact simulation algorithm for $X(1)$ that terminates with probability one, see Algorithm 2.

3 Exact Simulation for General SDEs

In this section, we will develop an exact simulation algorithm for the SDE (1). We shall fix $X\left(0\right)=x_{0}$ and the dependence of $x_{0}$ in some objects (such as the transition density of $X\left(1\right)$ will be omitted).

We are still going to construct an exact simulation algorithm based on acceptance-rejection in this section. However, for SDEs with non-constant diffusion matrix, applying Girsanov’s theorem no longer provides a Brownian type proposal distribution for acceptance-rejection, so we will construct an acceptance-rejection algorithm based on the density of $X(1)$.

Throughout the rest of this section, we shall assume the following assumptions and conditions.

{assumption}

The drift coefficient $\mu(\cdot)$ is continuously differentiable, and the diffusion coefficient $\sigma(\cdot)$ is three times continuously differentiable. Moreover, a strong solution to SDE (1) exists.

{condition}

The probability distribution of $X(1)$ is absolutely continuous with respect to Lebesgue measure. In other words, $X(1)$ has a density function denoted by $p(\cdot)$ with respect to the Lebesgue measure.

{condition}

For any relatively compact set $S$, the density $p(\cdot)$ is Lipschitz continuous with Lipschitz constant $C_{S}$, i.e.

[TABLE]

{condition}

For any relatively compact set $S$, there exist $\delta_{S}>0$ such that

[TABLE]

As we have seen in the previous section, Assumption 3 is the necessary condition for the applicability of the TES result introduced in Theorem 2.3, which enables us to strongly approximate $X(1)$. Condition 3 will eventually be used to apply the acceptance-rejection technique using an absolutely continuous (with respect to the Lebesgue measure) proposal distribution. Conditions 3 and 3, as we shall see, will allow us to control the bound of the likelihood ratio when applying acceptance-rejection.

It is important to ensure that the constants $C_{S}$ and $\delta_{S}$ are explicitly computable in terms of $\mu\left(\cdot\right)$ and $\sigma\left(\cdot\right)$ only, but we should also emphasize that we are not assuming that the density $p(\cdot)$ is known.

There are many ways in which the computability of $C_{S}$ and $\delta_{S}$ can be enforced. For instance, in Appendix A we discuss a set of assumptions involving classical estimators of the fundamental solutions of parabolic equations, which we review in order to compute $C_{S}$ and $\delta_{S}$ explicitly.

The standard use of the acceptance-rejection algorithm requires knowing the density $p(x)$, which seems hopeless for the general SDE problem that we study. An alternative approach is constructing a non-negative, bounded, and unbiased estimator of $p(x)$. While the density $p(x)$ is unknown, an unbiased estimator of $p(x)$, denoted by $\Lambda_{N}(x)$ in Section 3.1, can be constructed by means of a local approximation of the density. However, the unbiased estimator $\Lambda_{N}(x)$ may be negative, so it cannot be directly used in acceptance-rejection. To remedy this problem, in Lemma 3.5 we construct a non-negative and unbiased estimator $\Lambda^{+}_{N}(x)$ of $p(x)$ using a random walk and a suitable Bernoulli factory. However, the estimator $\Lambda^{+}_{N}(x)$ is unbounded, so we propose to sample enough information about the SDEs (the ancillary variable $N^{\prime}$), such that the estimator $\Lambda^{+}_{N}(x)$ conditional on the sampled information is locally bounded. Consequently, conditional on the localization of $X(1)$ and the additional information $N=N^{\prime}$, the estimator $\Lambda^{+}_{N}(x)$ is bounded and non-negative.

We now state the outline of our exact simulation algorithm. First of all, we apply a localization technique on the countable continuous partition $\mathcal{G}_{\text{loc}}$ defined as

[TABLE]

Since $\mathcal{G}_{\text{loc}}$ has countable components, we can enumerate $\mathcal{G}_{\text{loc}}$ and rewrite it in terms of $\mathcal{G}_{\text{loc}}=\{G_{i}:i\in\mathbb{N}\}$, where $G_{i}$ is a unit hypercube. Obviously, Algorithm 1 is applicable to $X(1)$ with respect to the countable continuous partition $\mathcal{G}_{\text{loc}}$.

Then, we introduce an ancillary random variable $N^{\prime}$ coupled with $X(1)$, and simulate $(N^{\prime}|X(1)\in G_{i})$. As we shall see, the random variable $N^{\prime}$ will play an important role after we introduce a suitable family of random variables whose expectations converge to the density of $X\left(1\right)$ at a given point. In the end, we will be able to sample $X(1)$ conditional on $N^{\prime}$ and $X(1)\in G_{i}$, using the estimator $\Lambda^{+}_{N}(x)$, which is bounded and non-negative for $x\in G_{i}$ and $N=N^{\prime}$.

The following theorem provides the main contribution of this paper.

Theorem 3.1**.**

If Assumption 3 and Condition 3-3 are satisfied, then there is an algorithm for exactly simulating $X(1)$ which terminates in finite time with probability one; see Algorithm 3.

The rest of this section is organized as follows. Section 3.1 applies a technique borrowed from Multilevel Monte Carlo to construct the unbiased density estimator and the ancillary random variable $N^{\prime}$. Section 3.2 explains how to sample $N^{\prime}$ using acceptance-rejection and a suitable Bernoulli factory [19, 17, 14] conditional on localization. Section 3.3 demonstrates how to sample $X(1)$ conditional on $N^{\prime}$, once again using a suitable localization.

3.1 A Multilevel Representation of the Density

In this section, we borrow an idea from Multilevel Monte Carlo [11] to construct an unbiased estimator for $p(\cdot)$, and we also introduce the ancillary random variable $N^{\prime}$.

In order to illustrate our idea, first we need to introduce some notations. For any $x$ in $G_{i}$, we define $\{B_{r_{n}}(x):n\geq 1\}$ as a sequence of open balls centered at $x$, whose radii $\{r_{n}:n\geq 1\}$, form a decreasing sequence and $r_{n}\rightarrow 0$ as $n\rightarrow\infty$.

Let $V(r)$ denote the volume of a $d$-dimensional ball with radius $r$ (i.e. the volume of $B_{r}\left(0\right)$). We define $\overline{p_{n}}(x)$ to be the average density over the ball $B_{r_{n}}(x)$, namely,

[TABLE]

Let $\hat{p}_{n}(x)$ denote a nonnegative unbiased estimator for $\overline{p_{n}}(x)$, i.e.

[TABLE]

for $n\geq 1$, where $\hat{p}_{n}(x)$ is defined using the same realization $X(1)$ for all $n$ and $x$. We define $\hat{p}_{0}(x):=0$ and $\overline{p_{0}}:=0$ for notational simplicity. It follows immediately that $\mathbb{E}[\hat{p}_{n}(x)]=\overline{p_{n}}(x)$ for $n\geq 0.$

The density $p(x)$ is first decomposed into an infinite telescoping sums,

[TABLE]

Then, we introduce a randomization technique inspired by Randomized Multilevel Monte Carlo (see [18, 22]). The density $p(x)$ can be decomposed as expectation of an infinite sum of estimators, which is truncated to a finite but random level so that the expectation is invariant. The idea is to introduce an integer-valued random variable $N$, which is independent of everything else. Then $p(x)$ can be expressed as

[TABLE]

where the third equality follows from Fubini’s theorem, which can be justified if

[TABLE]

We will show $\sum_{n=0}^{\infty}r_{n}<\infty$ in the sequel. Moreover, by the tower property we have

[TABLE]

Therefore, if we define

[TABLE]

it follows easily that

[TABLE]

We now are interested in obtaining bounds for $\Lambda_{n}(x)$ and its expectation for $x\in G_{i}$. To this end we first define $G_{i,r_{1}}$ as $r_{1}$-neighborhood of set $G_{i}$, which consists of all points that at a distance less than $r_{1}$ from $G_{i}$, i.e.

[TABLE]

It is not hard to observe that $G_{i,{r_{i}}}$ is a relatively compact set, to which Conditions 3-3 are applicable. In the following lemma, we will demonstrate that under such conditions, one can judiciously pick the distribution of $N$ and the radii $\{r_{n}:n\geq 1\}$ in order to establish explicit bounds for $\Lambda_{n}(x)$ and $\mathbb{E}[\Lambda_{n}(x)]$, respectively.

Lemma 3.2**.**

Suppose that $x\in G_{i}$ and Condition 3-3 are satisfied. If we pick

[TABLE]

for $n\geq 1$, then we have

[TABLE]

and

[TABLE]

Proof 3.3**.**

Let us construct the lower bound of $\mathbb{E}[\Lambda_{n}(x)]$ first. By triangle inequality,

[TABLE]

On the one hand, from the definition of $\Lambda_{0}(x)$ and Condition 3, we can conclude that

[TABLE]

On the other hand, using the triangle inequality

[TABLE]

Then, from Condition 3 we have

[TABLE]

Recall that $\overline{p_{k}}(x)$ is the average density over the ball $B_{r_{k}}(x)$ and $B_{r_{k}}(x)\subseteq B_{r_{1}}(x)\subseteq G_{i,r_{1}}$ for $x\in G_{i}$. It then follows form the (11) that

[TABLE]

Consequently,

[TABLE]

Combining the above inequality with $\mathbb{E}[\Lambda_{0}(x)]\geq\delta_{G_{i,r_{1}}}$ yields

[TABLE]

Similarly, observing that $\mathbb{E}[\Lambda_{0}(x)]=\mathbb{E}[\hat{p}_{1}(x)]=[V(r_{1})]^{-1}\times\mathbb{P}(X(1)\in B_{r_{1}}(x))\leq[V(r_{1})]^{-1}$, for the upper bound we have

[TABLE]

We can also derive an upper bound of $|\Lambda_{n}(x)|$:

[TABLE]

which ends the proof.

In the rest of this paper, we will adopt the value of $r_{n}$ and the distribution of $N$ in Lemma 3.2.

Even though we have constructed an unbiased estimator $\Lambda_{N}(x)$ for $p(x)$, acceptance-rejection is not applicable because $\Lambda_{N}(x)$ may be negative and unbounded. In order to apply acceptance-rejection, we need a nonnegative unbiased estimator for $p(x)$, which will be constructed in Lemma 3.5. The idea of such construction borrows an idea from [8]. Let $\{\Lambda_{n,k}(x):k\geq 1\}$ be i.i.d. copies of $\Lambda_{n}(x)$. We then define

[TABLE]

and

[TABLE]

By Wald’s first equation,

[TABLE]

Note that $S_{n,\tau_{n}(x)}(x)\geq 0$, but now we have an additional contribution $1/\mathbb{E}[\tau_{n}(x)]$, which can be interpreted as a probability. In order to sample a Bernoulli random variable with such probability, we will need the following result which we refer to as the Bernoulli Factory.

Theorem 3.4** (Bernoulli factory [19, 14]).**

Assume that $\epsilon\in(0,1/2]$ and $\alpha>0$ are two known constants and that we have an oracle that outputs i.i.d. Bernoullies with parameter $p\in(0,\left(1-\epsilon\right)/\alpha]$. Then, there is an algorithm which takes the output of the oracle and produces a Bernoulli random variable with parameter $f(p)=\min\left(\alpha p,1-\epsilon\right)=\alpha p$. Moreover, if $\bar{N}\left(\alpha,\epsilon\right)$ is the number of Bernoulli$\left(p\right)$ random variables required to output Bernoulli$\left(f\left(p\right)\right)$ then $.004\cdot\alpha/\epsilon\leq\mathbb{E}\left(\bar{N}\left(\alpha,\epsilon\right)\right)\leq 10\cdot\alpha/\epsilon$.

We now can explain how to construct $\Lambda_{n}^{+}(x)\geq 0$ such that $\mathbb{E}[\Lambda_{n}^{+}(x)]=\mathbb{E}[\Lambda_{n}(x)]$.

Lemma 3.5**.**

There exists a family of random variables $\{\Lambda_{n}^{+}(x):n\in\mathbb{N},x\in G_{i}\}$, such that the following properties hold:

$0\leq\Lambda_{n}^{+}(x)\leq m_{n}$. 2. 2.

$\mathbb{E}[\Lambda_{n}^{+}(x)]=\mathbb{E}[\Lambda_{n}(x)]$. 3. 3.

Given $n$ and $x$, there is an algorithm for simulating $\Lambda_{n}^{+}(x)$.

Proof 3.6**.**

Let $\bar{\Gamma}_{n}(x)$ be a Bernoulli random variable with parameter $\left(\mathbb{E}[\tau_{n}(x)]\right)^{-1}$, and independent of everything else. It follows that

[TABLE]

We write $\Lambda_{n}^{+}(x):=\bar{\Gamma}_{n}(x)S_{n,\tau_{n}(x)}(x)$. Property 1 follows from the facts that $0\leq S_{n,\tau_{n}(x)}(x)\leq\Lambda_{n,\tau_{n}(x)}(x)\leq m_{n}$, and that $0\leq\bar{\Gamma}_{n}(x)\leq 1$. Property 2 is justified directly by equation (12). To show that $\Lambda_{n}^{+}(x)$ can be simulated, we just need to provide an algorithm for simulating $\bar{\Gamma}_{n}(x)$.

Recall that $\mathbb{E}[\Lambda_{n}(x)]\geq\delta_{G_{i,r_{1}}}/2$, we have

[TABLE]

Consider Wald’s second equation

[TABLE]

which implies

[TABLE]

Now we shall proceed to simulate the random variable $\bar{\Gamma}_{n}(x)$. Consider a random variable $T_{n}(x)$ with distribution

[TABLE]

Since $I(T_{n}(x)=1)$ is the desired Bernoulli random variable with parameter $(\mathbb{E}[\tau_{n}(x)])^{-1}$, it then suffices to simulate $T_{n}(x)$. Towards this end, we apply acceptance-rejection again using another random variable $T^{\prime}$ as proposal, whose distribution is

[TABLE]

Since $\tau_{n}(x)$ is nonnegative, Markov’s inequality asserts that

[TABLE]

where the last inequality follows from (13). Also, from the definition of $\tau_{n}(x)$ we know that $\mathbb{E}[\tau_{n}(x)]\geq 1$. Consequently, the likelihood ratio between $T_{n}(x)$ and $T^{\prime}$ is given by

[TABLE]

From the above inequality we see that the likelihood ratio is bounded, so the acceptance-rejection procedure is applicable.

Conditional on the proposal $T^{\prime}=k$, we then introduce a new Bernoulli random variable $\tilde{\Gamma}_{n,k}(x)$ to decide whether or not the proposal is accepted as $T_{n}(x)$. The distribution of $\tilde{\Gamma}_{n,k}(x)$ is defined as

[TABLE]

Hence it follows from (14) that

[TABLE]

Observe that the indicator $I(\tau_{n}(x)\geq k)$ is simulable, but its distribution is not explicitly accessible, so it is natural to sample $\tilde{\Gamma}_{n,k}(x)$ via Bernoulli factory introduced in Theorem 3.4. Due to (15), the function $f(\cdot)$ involved in Bernoulli factory is a linear function as following

[TABLE]

We summarize the procedure of simulating $\Lambda_{n}^{+}(x)$ in Algorithm 4.

We now introduce an ancillary random variable $N^{\prime}$ coupled with $X(1)$ in the following way,

[TABLE]

Assuming $(N^{\prime}|X(1)\in G_{i})$ can be simulated, $(X(1)|N^{\prime},X(1)\in G_{i})$ can be easily simulated by acceptance-rejection as well due to the convenient density representation given by (18). The algorithm for sampling $(N^{\prime}|X(1)\in G_{i})$ will be explained in the next section.

3.2 Conditional Sampling of $N^{\prime}$

In this section we will focus on the procedure for simulating $N^{\prime}$ conditional on $X(1)\in G_{i}$.

First we derive from equation (16) the probability mass function of $(N^{\prime}|X(1)\in G_{i})$, namely

[TABLE]

Due to the upper bound of $\mathbb{E}[\Lambda_{n}]$ given by Lemma 3.2, we have the following inequality

[TABLE]

Simulation of $(N^{\prime}|X(1)\in G_{i})$ can be achieved by acceptance-rejection. Consider a Bernoulli random variable $\Gamma_{n}(x)$ defined as

[TABLE]

Then the outline of the acceptance-rejection algorithm for simulating $N^{\prime}$ would be:

Step 1

Sample $n$ from random variable $N$.

Step 2

Sample $x$ from uniform distribution $U_{G_{i}}\sim\text{Unif}(G_{i})$.

Step 3

Simulate $\Gamma_{n}(x)$. If $\Gamma_{n}(x)=0$, go to Step

1. Otherwise accept $n$ as a sample of $N^{\prime}$.

The only difficult step in the above procedure is Step 3, namely, sample $\Gamma_{n}(x)$, which we will discuss now.

Lemma 3.2 implies that $0\leq\Lambda_{n}^{+}(x)\leq m_{n}$. Let $U\sim\text{Unif}(0,m_{n})$, which is independent of everything else, then $I(U\leq\Lambda_{n}^{+}(x))$ is a Bernoulli random variable with parameter $(m_{n})^{-1}\mathbb{E}[\Lambda_{n}^{+}(x)]$. Due to (17), Bernoulli factory given in Theorem 3.4 is applicable to simulate $\Gamma_{n}(x)$, using $I(U\leq\Lambda_{n}^{+}(x))$ as input and using the function

[TABLE]

To conclude, we synthesize the complete steps for simulating $N^{\prime}$ in the following algorithm.

3.3 Sampling of $(X(1)|N^{\prime},X(1)\in G_{i})$

In this section, we will focus on sampling $(X(1)|N^{\prime},X(1)\in G_{i})$.

Without loss of generality, let us assume $N^{\prime}=n$ throughout the the rest of this section. According to equation (16) and Lemma 3.5, the conditional distribution of $X(1)$ can be written as

[TABLE]

where $C_{n,G_{i}}$ is a constant independent of $x$. Once again, as we shall see $(X(1)|N^{\prime}=n,X(1)\in G_{i})$ can be be simulated by acceptance-rejection. We use the uniform distribution $U_{G_{i}}$ as the proposal distribution, and accept the proposal with probability $m_{n}^{-1}\times\mathbb{E}[\Lambda_{n}^{+}(x)]$, so we can accept if and only if $I\left(U\leq\Lambda_{n}^{+}(x)\right)=1$, where $U\sim\text{Unif}(0,m_{n})$ is independent of everything else. The output of the acceptance-rejection follows the desired distribution. The explicit procedure for simulating $(X(1)|N^{\prime},X(1)\in G_{i})$ is given in the following algorithm.

4 Conclusion

The main contribution of this paper is the construction of the first generic exact simulation algorithm for multidimensional diffusions. The algorithm extensively uses several localization ideas and the application of TES techniques. But it also combines ideas from multilevel Monte Carlo in order to construct a sequence of random variables which ultimately provides an unbiased estimator for the underlying transition density.

Although the overall construction can be implemented with a finite termination time almost surely, the expected running time is infinite. Thus, the contribution of the paper is of theoretical nature, showing that it is possible to perform exact sampling of multivariate diffusions without applying Lamperti’s transformation. However, more research is needed in order to investigate if the algorithm can be modified to be implemented in finite expected time, perhaps under additional assumptions. Alternatively, perhaps some controlled bias can be introduced while preserving features such as positivity and convexity in the applications discussed at the end of the Introduction. To this end, we conclude with a discussion of the elements in the algorithm which are behind the infinite expected termination time.

There are three basic problems that cause the algorithm to have infinite expected termination time. Two of them can be appreciated already from the constant diffusion discussion and involves the use of the localization techniques that we have introduced. The third problem has to do with the application of the Bernoulli factory.

$\bullet\quad$Problem 1: The first problem arises when trying to sample a Bernoulli of the form $I\left(X\left(1\right)\in G\right)$. Given $\epsilon_{n}>0$, sampling $X_{\epsilon_{n}}\left(1\right)$ such that $\left\|X_{\epsilon_{n}}\left(1\right)-X\left(1\right)\right\|\leq\epsilon_{n}$ takes $O(\epsilon_{n}^{-\left(2+\delta\right)})$ computational cost for any $\delta>0$. So, if $G$ is a unit hypercube in $d$ dimensions, using the density estimates for $X\left(1\right)$, we obtain

[TABLE]

for some $c_{G}>0$. Therefore, if $N_{0}$ is the computational cost required to sample $I\left(X\left(1\right)\in G\right)$ we have that for some $\delta_{0}>0$

[TABLE]

Therefore,

[TABLE]

which yields that $\mathbb{E}\left(N_{0}\right)=\infty$.

$\bullet\quad$Problem 2: The second problem arises in the acceptance-rejection step applied in Lemma 2.9, which requires sampling $X\left(1\right)$ under $\tilde{\mathbb{P}}$ conditional on $\max(i,(i+1)U)<L(1)<i+1$. Directly sampling from this conditional law might be inefficient if $i$ is large. However, this problem can be mitigated using rare-event simulation techniques, which might be available in the presence of additional structure on the drift because under $\tilde{\mathbb{P}}\left(\cdot\right)$, $X\left(\cdot\right)$ follows a Brownian motion.

$\bullet\quad$Problem 3: Arises because, as indicated in Theorem 3.4, the computational complexity of the Bernoulli factory of a linear function of the form $f\left(p\right)=\min\left(\alpha p,1-\epsilon\right)$ scales in order $O\left(\alpha/\epsilon\right)$. In our case, we are able to select $\epsilon=1/2$ and we invoke the Bernoulli factory in Algorithms 4 and 5. In Algorithm 4, $\alpha=O\left(k^{2}\right)$, given $T^{\prime}=k$ and $\mathbb{E}\left(T^{\prime}\right)=\infty$. In Algorithm 5, $\alpha=O\left(m_{n}\right)$, given $N=n$. Although the bound which is used to define $m_{n}$ in Lemma 3.5 is far from optimal, in its current form, $m_{n}=O\left(n^{3d+2}\right)$, we have that $\mathbb{E}\left(N\right)=\infty$.

\acks

We gratefully acknowledge support from the following NSF grants 1915967, 1820942, 1838676 as well as AFOSR.

Appendix A Transition Density Estimates

In the appendix, we will discuss some assumptions which are sufficient for the applicability of Conditions 3, 3, and 3. In addition, we also give explicit procedures for computing the constants which appears in such Conditions.

Consider a matrix valued function $a(\cdot)=(a_{ij}(\cdot))_{d\times d}:\mathbb{R}^{d}\mapsto\mathbb{R}^{d\times d}$ defined as

[TABLE]

{assumption}

Every component of $\mu$ and $a$ are three times continuously differentiable. Moreover, there exist a constant $M$ such that $\|\mu^{(i)}\|_{\infty}\leq M$ and $\|a^{(i)}\|_{\infty}\leq M$ for $i=0,1,2,3$.

{assumption}

There exist constants $0<\lambda_{\downarrow}<\lambda_{\uparrow}<\infty$, such that for all $x\in\mathbb{R}^{d}$ and $\xi=(\xi_{i})_{d}\in\mathbb{R}^{d}$, we have

[TABLE]

Under Assumption A and A, it is proved in [10] that the SDE (1) possesses a transition density denoted by $p(x,t;y,\tau)$, which satisfies

[TABLE]

for $\tau<t$. Therefore, Condition 3 is proved given assumptions A and A.

In the following Proposition, we will establish Condition 3 via Kolmogorov forward equation.

Proposition A.1**.**

Suppose Assumptions A and A are satisfied, then for any relatively compact set $S$, the density $p(\cdot)=p(\cdot,T;x_{0},0)$ is Lipschitz continuous with Lipschitz constant $C_{S}$, i.e.

[TABLE]

Furthermore, $C_{S}$ can be computed by Algorithm 7.

Proof A.2**.**

Our methodology is closely related to parametrix method introduced in [10]. Following the same scheme, we focus on explicit computation of the constants.

Under Assumptions A and A, $p(\cdot,\cdot;y,\tau)$ is a solution of Kolmogorov forward equation, namely,

[TABLE]

We shall rewrite equation (19) into its non-divergence form as

[TABLE]

where

[TABLE]

and $\mathcal{L}_{f}$ is a uniform parabolic operator. By Assumption A, it follows that

[TABLE]

We denote by $a^{-1}(x)=(a_{ij}^{-1}(x))_{d\times d}$ the inverse matrix of $(a_{ij}(x))_{d\times d}$, and define

[TABLE]

Assumption A implies that for all $\xi\in\mathbb{R}^{d}$

[TABLE]

Following the idea of the parametrix method, we also define a partial differential equations with constant coefficients, namely,

[TABLE]

The fundamental solution of function $\mathcal{L}_{0}^{y}u=0$ is given by

[TABLE]

where

[TABLE]

[TABLE]

According to Theorem 1.3 and Theorem 1.10 in [10], $p(x,t;\xi,\tau)$ can be represented by the parametrix method as

[TABLE]

where

[TABLE]

for $k\geq 1$. Furthermore, the partial derivatives of the fundamental solution admit the following expression,

[TABLE]

Let us pick $\epsilon\in(0,1)$, then we can derive a bound for $Z$ as

[TABLE]

For the bound of $\partial Z(x,t;\xi,\tau)/\partial x_{i}$, note that

[TABLE]

and that

[TABLE]

Combining the definition of $C_{Z}(\cdot)$ and the above two equations imply that

[TABLE]

Applying the inequalities

[TABLE]

and

[TABLE]

we obtain

[TABLE]

by setting

[TABLE]

Similarly, we can derive a bound for $\partial^{2}Z(x,t;\xi,\tau)/\partial x_{i}\partial x_{j}$ and $\partial^{2}Z(x,t;\xi,\tau)/\partial x_{i}^{2}$. For $i\neq j$,

[TABLE]

where

[TABLE]

By definition of $Z(\cdot)$ we can observe that

[TABLE]

Suppose $0\leq t-\tau\leq T$ in the sequel. By considering the upper bounds of partial derivatives of $Z$, as well as (20) and Assumption A, we obtain

[TABLE]

where

[TABLE]

Now, in order to find a bound for $\Phi(x,t;\xi,\tau)$, we need to introduce a technical lemma.

Lemma A.3** (Lemma 1.3 of [10]).**

If $\beta$ and $\gamma$ are two constants in $(-\infty,\frac{d}{2}+1)$, then

[TABLE]

where Beta$(\cdot)$ is Beta function.

Due to (22) and Lemma A.3, we can derive

[TABLE]

where

[TABLE]

By induction we can show that, for any positive integer $m$,

[TABLE]

It turns out that

[TABLE]

where

[TABLE]

Recalling equation (21), we can apply Lemma A.3 again and conclude that

[TABLE]

where

[TABLE]

Therefore, we obtain an upper bound for $|\nabla_{x}p(x,t;\xi,\tau)|$, by considering

[TABLE]

Therefore, for all $x,y\in S$ we have

[TABLE]

where

[TABLE]

Next, we will propose a computational procedure for lower bounds of transition density. There is a substantial amount of literature that studies lower bounds for the transition density of diffusions, through analytical approaches or probabilistic approaches. For instance, Aronson [2] develops estimates of lower bounds of fundamental solutions of second order parabolic PDEs in divergence form. Using Malliavin calculus, Kusuoka and Stroock [16] derived a lower bound for the transition density of uniformly elliptic diffusions. Bally [3] generalized the idea of [16] to locally elliptic diffusions. We follow the approach suggested by Sheu [23] and review it in order to find explicit expressions to obtain a computable lower bound.

In order to keep our paper self-contained, first we need to introduce some notations used later.

Let $\mathcal{L}_{b}$ be the generator of Komolgorov backward equation:

[TABLE]

The transition density as a function of $(x,t)\mapsto p(y,t;x,0)$ coincides with the fundamental solution of Komolgorov backward equation:

[TABLE]

Throughout the rest of this section, we suppose that Assumptions A and A are in force.

Recall that $a^{-1}(x)$ is the inverse matrix of $a(x)$, and define

[TABLE]

For a fixed $y_{0}\in\mathbb{R}^{d}$, we define

[TABLE]

and

[TABLE]

The continuity of the density implies

[TABLE]

For simplicity, we also define the logarithmic transform of $p$ and $p^{\beta}$ as

[TABLE]

To prepare the analysis, which is based on stochastic control, we introduce the space of control functions by $\mathcal{F}_{T,x}$. The class $\mathcal{F}_{T,x}$ is defined as a family of measurable functions $\psi:[0,T]\times\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ such that the SDE

[TABLE]

has a weak solution $\eta(\cdot)$ that satisfies

[TABLE]

Now we state a lemma that is crucial for proving the main result of this section.

Lemma A.4**.**

Recall the definition of $\mathcal{F}_{T,x}$ and $\eta(\cdot)$ from previous paragraph, then we have

[TABLE]

Together with (23), we see that

[TABLE]

Proof A.5**.**

See [9].

Theorem A.6**.**

Suppose Assumptions A and A are satisfied. Then, for any relatively compact set $S$, the density $p(\cdot)=p(\cdot,T;x_{0},0)$ has a uniform lower bound $\delta_{S}>0$ in $S$, i.e.

[TABLE]

Furthermore, $\delta_{S}$ can be computed by Algorithm 8.

Proof A.7**.**

Finding a lower bound of density $p(y_{0},T;x_{0},0)$ is equivalent to finding an upper bound for $J(T,x_{0})$. Towards this end, it suffices to find an upper bound for $J^{\beta}(T,x_{0})$ that is uniform in $\beta$. We shall define $\phi(\cdot)$ as a linear function, such that $\phi(0)=x_{0},\phi(T)=y_{0}$. Write

[TABLE]

It is not hard to see that $\psi\in\mathcal{F}_{T-\beta,x_{0}}$. Therefore,

[TABLE]

according to lemma A.4. Notice that

[TABLE]

It follows that

[TABLE]

and

[TABLE]

We now apply a Taylor expansion of $k(\eta(t),\psi(t))$ around $(\phi(t),\dot{\phi}(t))$, where $\dot{\phi}(\cdot)$ denotes the derivative of $\phi(\cdot)$. For notational simplicity, we define

[TABLE]

We also define

[TABLE]

*and similarly for $D_{\psi_{i}}k$, $D_{x_{i},x_{j}}k$, $D_{x_{i},\psi_{j}}k$ and $D_{\psi_{i},\psi_{j}}k$. The Taylor expansion with remainders of third order is given as following *

[TABLE]

where

[TABLE]

Now we integrate all the above terms from [math] to $T-\beta$ with respect to variable $t$, then take expectations, and analyze the upper bounds of the result term by term.

$\bullet$** Zeroth Order Term: Notice that $k$ is in quadratic form, with matrix $(a_{ij}^{-1}(x))$, so**

[TABLE]

$\bullet$** First Order Terms: We treat first order term $k_{1}(t)$ first. Noting that $\Delta_{2,i}(t)$ is a martingale due to (26), the first order terms**

[TABLE]

$\bullet$** Second Order Terms: We then treat the second order terms. As $D_{\psi_{i},\psi_{j}}k(0)=a_{ij}^{-1}(\phi(t))$,**

[TABLE]

Writing

[TABLE]

and noticing that $(a_{ij}(t))$ is symmetric, we see that

[TABLE]

Assumption A implies the Lipschitz continuity of $a(\cdot)$, which gives,

[TABLE]

Due to (27) and Jensen’s inequality,

[TABLE]

Observe that $\sum_{i,j=1}^{d}|a_{ij}^{-1}(\phi(t))|\leq d\lambda_{\downarrow}^{-1}$, so we have

[TABLE]

By the Lipschitz continuity of $(a_{ij}(\cdot))$, it follows that

[TABLE]

Therefore,

[TABLE]

Combining (28), (29) and (30) yields

[TABLE]

By the chain rule and Assumption A, we obtain

[TABLE]

where $\Psi_{i}(\cdot):\mathbb{R}\rightarrow\mathbb{R};i=1,2$ are defined as

[TABLE]

Taking (27) into consideration, we obtain,

[TABLE]

$\bullet$** Third Order Terms: We proceed to analyze the third order remainder terms. Let us consider $k_{3,3}(t)$ first. Notice that**

[TABLE]

thus,

[TABLE]

Then, by Burkholder-Davis-Gundy inequality,

[TABLE]

where $C_{BDG}(3)$ is the explicit constant in the Burkholder-Davis-Gundy inequality. We can pick $C_{\text{BDG}}(p)=\left(\frac{p(p-1)}{2}(\frac{p}{p-1})^{p}\right)^{p/2}$ (See Proposition 4.4.3 of [20].). To summarize, we obtain

[TABLE]

Next, we consider the other two remainders $k_{3,2}(t),k_{3,1}(t)$. Observe that

[TABLE]

and

[TABLE]

Thus, by a similar argument, we can also derive

[TABLE]

and

[TABLE]

where $\Psi_{3}(\cdot)$ and $\Psi_{4}(\cdot)$ are defined as

[TABLE]

[TABLE]

Finally, let us consider $\mathbb{E}(J^{\beta}(0,\eta(T-\beta)))$. Since

[TABLE]

It follows that

[TABLE]

To conclude, let us summarize all the intermediate results, and substitute them into (25), we have

[TABLE]

where $J_{\uparrow}(\cdot;T)$ is defined as

[TABLE]

Therefore, if we pick $D_{S}=\sup_{x\in S}\|x-x_{0}\|_{2}$, it follows that

[TABLE]

which ends the proof.