<figure><img src="image_3.png" /><figcaption>Figure 4.1: On the left: Weak discretization error $n \mapsto nh^n(\theta^*)$. On the right: Implicit discretization error $n \mapsto n(\theta^{*,n} - \theta^*)$, $n = 100, \dots, 500$.</figcaption></figure>

where for all $\theta \in \mathbb{R}^d$, $\Lambda_1^l(\theta), \dots, \Lambda_{k-l}^l(\theta)$ and $o(n^{-\alpha(k-l)})$ are multilinear maps from $(\mathbb{R}^d)^l$ to $\mathbb{R}^d$.

3. For all $1 \le l \le k$, $(D^l h^n)_{n \ge 1}$ converges locally uniformly towards $D^l h$.

4. $Dh(\theta^*)$ is invertible.

Under (H-1) we have seen in Proposition 4.1.1 that a first order expansion of $\theta^{*,n} - \theta^*$ holds, that is, $\theta^{*,n} - \theta^* = C_1 n^{-\alpha} + o(n^{-\alpha})$. In [10], we obtained a generalisation of this first result.

**Theorem 4.3.1.** Assume that $\theta^{*,n} \to \theta^*$, $n \to +\infty$, and that [H-p] holds for some $p \in \mathbb{N}^*$. Then, $\theta^{*,n} - \theta^*$ has an expansion up to order $p$, that is, the following expansion holds:

$$ \theta^{*,n} - \theta^* = \frac{C_1}{n^\alpha} + \dots + \frac{C_p}{n^{\alpha p}} + o\left(\frac{1}{n^{\alpha p}}\right). $$

From the previous expansion, one is able to develop a multi-step Richardson-Romberg extrapolation method for stochastic approximation algorithm. The idea is the same as for the Richardson-Romberg method in the context of Monte Carlo linear estimator for the computation of $\mathbb{E}[f(X_T)]$, see Section 1.2.3 for a brief description and Pageès [109] for a complete overview.

We proceed as follows. Let $R \ge 2$ be an integer. To devise a SA estimator whose target has an implicit discretisation error of order $n^{-\alpha R}$ as $n \to +\infty$, we introduce a sequence of $R$ random vectors $\{U^{rn}, 1 \le r \le R\}, n \in \mathbb{N}^*$. It is assumed that this sequence satisfies $U^{rn} \xrightarrow{\mathbb{P}} U^r$ with $U^r \stackrel{d}{=} U$, $1 \le r \le R$, all variables being defined on the same probability space. If assumption [H-R] holds then for $r = 1, \dots, R$ one gets

$$ \theta^{*,rn} = \theta^* + \sum_{p=1}^{R-1} \frac{C_p}{r^{\alpha p} n^{\alpha p}} + \frac{C_R}{r^{\alpha R} n^{\alpha R}} (1 + \varepsilon_r(n)), \quad \varepsilon_r(n) \to 0. $$

By solving a Vandermonde system, one finds an $Rd \times d$ weight matrix $\mathbf{w} = (\mathbf{w}_1, \dots, \mathbf{w}_R)^T$, $\mathbf{w}_i$ being a $d \times d$ matrix $i = 1, \dots, R$ such that $\sum_{r=1}^R \mathbf{w}_r = I_d$ and

$$ \sum_{r=1}^{R} \mathbf{w}_r \theta^{*,rn} = \theta^* + C_R \frac{1}{n^{\alpha R}} \frac{(-1)^{R-1}}{R!^{\alpha}} (1 + \varepsilon_{R+1}(n)) \quad (4.6) $$

Hence, we removed the $R-1$ first terms of the expansion by considering a suitable linear combination of the $(\theta^{*,rn})_{1 \le r \le R}$. We now approximate the new target $\sum_{r=1}^R \mathbf{w}_r \theta^{*,rn}$, by means of $M \in \mathbb{N}^*$ steps of $R$ stochastic algorithms with dynamics

$$ \theta_{p+1}^{rn} = \theta_p^{rn} - \gamma_{p+1} H(\theta_p^{rn}, U^{rn,(p+1)}), \quad 0 \le p \le M-1, \ r = 1, \dots, R, \quad (4.7) $$
<figure><img src="image_4.png" /><figcaption>Figure 4.2: Complexity with respect to RMSE.</figcaption></figure>
where $(U^{rn,(p)}, r = 1, \dots, R)_{1 \le p \le M}$ is an i.i.d sequence with the same law as $(U^{rn}, r = 1, \dots, R)$, $\theta_0^{rn}$, $r = 1, \dots, R$ are the initial conditions independent of the innovation sequence satisfying $\sup_{n \ge 1} |\mathbb{E}|\theta_0^{rn}|^2 < +\infty$ and the sequence $(\gamma_p)_{p \ge 1}$ satisfies (1.9). Now the new statistical error of the Richardson-Romberg extrapolation estimator writes

$$ \mathcal{E}_S^{R-R}(n, M) := \sum_{r=1}^{R} \mathbf{w}_r (\theta^{*,rn} - \theta_M^{rn}). $$

Under mild assumptions, we prove that this statistical error satisfies

$$ \mathbb{E}[|\mathcal{E}_S^{R-R}|] \le \nu_R \gamma^{1/2}(M) (1 + \phi_1^R(n) + \phi_2^R(M)) $$

where $\nu_R := C \mathbb{E} \left[ \left| \sum_{r=1}^R \mathbf{w}_r H(\theta^*, U^r) \right|^2 \right]^{1/2}$, $\phi_1^R, \phi_2^R$ are two positive functions satisfying: $\phi_1^R(n) \to 0$ and $\phi_2^R(M) \to 0$ respectively as $M \to +\infty$, $n \to +\infty$ and $\phi_2^R$ is non-increasing. Hence, the global error of the Richardson-Romberg estimator is bounded by $\mathbb{E}[|\mathcal{E}_{glob}^{R-R}|] \le \mu_R n^{-\alpha R} (1 + |\epsilon_{R+1}(n)|) + \nu_R \gamma^{1/2}(M) (1 + \phi_1^R(n) + \phi_2^R(M))$ with $\mu_R := \frac{|C_R|}{(R!n^R)^\alpha}$.

In practical situations, one is looking for an efficient estimator among the family $\{\sum_{r=1}^R \mathbf{w}_r \theta_M^{rn}, (n, M) \in (\mathbb{N}^*)^2\}$. To be more precise, one is interested in minimising the computational cost for a given $L^1(\mathbb{P})$-error $\varepsilon > 0$. We assume that the cost of a single simulation of $U^n$ is proportional to $n$ and is given by $K \times n$, where $K$ is a generic positive constant independent of $n$. It notably corresponds to the case of discretisation schemes of a stochastic process. In the case of the Richardson-Romberg method for SA, at each step $p = 1, \dots, M$ of the procedure, for $r = 1, \dots, R$, one has to simulate the random vector $(U^n, U^{2n}, \dots, U^{Rn})$ so that the global computational cost is given by

$$ \text{Cost}(\text{R-R}) := KM \sum_{r=1}^{R} rn = KMn \frac{R(R+1)}{2}. $$

Hence the problem of interest writes

$$ (n(\epsilon), M(\epsilon)) = \arg \min_{\mu_R n^{-\alpha R} (1 + |\epsilon_{R+1}(n)|) + \nu_R \gamma^{1/2}(M) (1 + \phi_1^R(n) + \phi_2^R(M)) \le \epsilon} \text{Cost}(\text{R-R}). $$

We solve the above computational optimisation problem asymptotically. More precisely, if $\gamma(p) = \gamma_0/p^\beta$, $\gamma_0 > 0$, $p > 0$, $\beta \in (1/2, 1]$. The multi-step Richardson-Romberg estimator of order $R$ satisfies

$$ \inf_{\mu_R n^{-\alpha R} (1 + |\epsilon_{R+1}(n)|) + \nu_R \gamma^{1/2}(M) (1 + \phi_1^R(n) + \phi_2^R(M)) \le \epsilon} \text{Cost}(\text{R-R}) \sim K \frac{R(R+1)}{2} \gamma_0^{\frac{1}{\beta}} \nu_R^{\frac{2}{\beta}} \mu_R^{\frac{1}{\alpha R}} \frac{1}{\varepsilon^{\frac{2}{\beta} + \frac{1}{\alpha R}}} \left(1 + \frac{2\alpha R}{\beta}\right)^{\frac{1}{\alpha R}} \left(1 + \frac{\beta}{2\alpha R}\right)^{\frac{2}{\beta}} $$

as $\varepsilon \to 0$. Eventually this asymptotically optimal bound may be achieved with parameters satisfying:

$$ n(\varepsilon) \sim \left( \frac{2\alpha R}{\beta} + 1 \right)^{\frac{1}{\alpha R}} \mu_R^{\frac{1}{\alpha R}} \varepsilon^{-\frac{1}{\alpha R}} \quad \text{and} \quad M(\varepsilon) \sim \gamma_0^{\frac{1}{\beta}} \nu_R^{\frac{2}{\beta}} \left(1 + \frac{\beta}{2\alpha R}\right)^{\frac{2}{\beta}} \varepsilon^{-\frac{2}{\beta}} \quad \text{as } \varepsilon \to 0. \quad (4.8) $$

Let us note that when one decides to implement the Richardson-Romberg extrapolation scheme with an innovation satisfying $U^r = U$ a.s. $r = 1, \dots, R$ then one has $H(\theta^*, U^r) = H(\theta^*, U)$ a.s. for every $r = 1, \dots, R$ which in turn yields

$$ \mathbb{E} \left[ \left| \sum_{r=1}^{R} \mathbf{w}_r H(\theta^*, U^r) \right|^2 \right] = \mathbb{E} \left[ \left| \left( \sum_{r=1}^{R} \mathbf{w}_r \right) H(\theta^*, U) \right|^2 \right] = \mathbb{E} \left[ |H(\theta^*, U)|^2 \right]. $$

Hence we clearly see that this choice leads to a control in the $L^1$-norm of the statistical error of the multi-step Richardson-Romberg SA estimator. On the opposite considering mutually independent innovations $U^r$ lead to an explosion of the previous control with respect to $R$. Using the exact expression of the weights, one has

$$ \begin{aligned} \mathbb{E} \left[ \left| \sum_{r=1}^{R} \mathbf{w}_r H(\theta^*, U^r) \right|^2 \right] &= \left( \sum_{r=1}^{R} \frac{r^{2\alpha R}}{\prod_{j=0}^{r-1} (r^\alpha - j^\alpha)^2 \prod_{j=r+1}^{R} (j^\alpha - r^\alpha)^2} \right) \mathbb{E} \left[ |H(\theta^*, U)|^2 \right] \\ &\ge \left( \frac{R^R}{R!} \right)^{2\alpha} \mathbb{E} \left[ |H(\theta^*, U)|^2 \right] \\ &\sim \left( \frac{e^R}{\sqrt{2\pi\sqrt{R}}} \right)^{2\alpha} \mathbb{E} \left[ |H(\theta^*, U)|^2 \right] \quad \text{as } R \to +\infty, \end{aligned} $$
For instance when one is concerned with the discretisation of a Brownian diffusion, the first aforementioned case consists in implementing the Richardson-Romberg method with R Euler schemes devised with the same Brownian motion W namely $W^r = W$, $r = 1, \dots, R$ whereas the second case consists in implementing the method with mutually independent Brownian motions $W^r$. The optimality of this choice is discussed in Pagès [109].

Finally, we illustrate our results by applying the Richardson-Romberg method to the estimation of the $\ell$-quantile ($\ell \in (0,1)$) $\theta^* = \inf\{\theta \in \mathbb{R} : \mathbb{P}_x(X_1^d \le \theta) \ge \ell\}$ where X is a d-dimensional stochastic process solution to a stochastic differential equation (1.17) where the underlying process Z is a symmetric $\alpha$-stable, for $\alpha \in (0,2]$. Note that the case $\alpha = 2$ corresponds to the standard Brownian motion.

The mean function h is defined by $h(\theta) = \mathbb{E}[H(\theta, U)]$, $H(\theta, u) = \mathbf{I}_{\{u \le \theta\}}$ and $U = X_1^d$. We approximate U by $U^n = X_1^{n,d}$ where $X^n$ is the continuous Euler scheme associated to X.

In order to prove an expansion in powers of $1/n$ for the implicit discretisation error we first remark that $(h-h^n)(\theta) = \mathbb{P}^x(X_1^d \le \theta) - \mathbb{P}^x(X_1^{n,d} \le \theta)$ so that $\frac{d^k}{d\theta^k}h(\theta) - \frac{d^k}{d\theta^k}h^n(\theta) = \frac{\partial^{k-1}}{\partial\theta^{k-1}}p^x(0, 1, x, \theta) - \frac{\partial^{k-1}}{\partial\theta^{k-1}}p_n^{x,n}(0, 1, x, \theta)$. Hence, one has to rely on the smoothness of the law of $X_1^d$ and $X_1^{n,d}$ in order to verify [H-k].

We rely on the parametrix expansion for the error $p^X(0, 1, x, \theta) - p_n^{X,n}(0, 1, x, \theta)$ between the densities of the process X and its Euler approximation $X^n$ obtained by Konakov and Mammen [76] in the case of Brownian SDEs and by Konakov and Menozzi [77] in the case of stable SDEs (see Section 1.2.3 for a brief presentation) when the coefficients b, $\sigma$ of the SDE are assumed to be smooth and bounded. Moreover, we assume that the diffusion matrix is uniformly elliptic and the spectral measure of Z when $\alpha < 2$ is non-degenerated. When $\alpha \le 1$, we put $b = 0$.

Note that we require more than an expansion of the error between the two densities since we also need an expansion of their derivatives in order to verify [H-k]. As a simple extension of the results of Theorem 1.1 in Konakov and Mammen in [76], for the Brownian case, and Theorem 21 in Konakov and Menozzi in [77], for the Stable case, we derive in [10] an expansion in powers of $n^{-1}$ for $\partial_\theta p^X(0, 1, x, \theta) - \partial_\theta p_n^{X,n}(0, 1, x, \theta)$. Hence, [H-k] is satisfied and we obtain an expansion of $\theta^{*,n} - \theta^*$ in powers of $n^{-1}$ at any order. Our result notably extend the error bound for the computation of the quantile of a diffusion process obtained by Talay and Zheng [123].

To illustrate the method we consider a geometric Brownian motion $(X_t)_{t \in [0,T]}$. Its dynamics and $\ell$-quantile are given by

$$ \forall t \in [0, T], \quad X_t = x_0 \exp\left((r - \sigma^2/2)t + \sigma W_t\right), \quad \theta^* = x_0 \exp\left((r - \sigma^2/2)t + \sigma \sqrt{t} \phi^{-1}(\ell)\right) $$

where W is a one dimensional standard Brownian motion, $\phi$ is the distribution function of the standard normal distribution $N(0, 1)$. We use the following values for the parameters: $x_0 = 100$, $r = 0.05$, $\sigma = 0.4$, $T = 1$, $\ell = 0.7$. The reference Black-Scholes quantile is $\theta^* = 119.69$. We set $\gamma(p) = \gamma_0/p$ with $\gamma_0 = 60$.

We plot in Figure 4.3 the behaviours of $\sum_{r=1}^R \mathbf{w}_r \theta^{*,rn} - \theta^*$ for $R = 2, 3, 4$ and $n = 2, \dots, 15$. We estimate $\theta^{*,rn}$ by $\theta_M^{rn}$, with $M = 10^6$ samples for $R = 2$ and $M = 10^8$ samples for $R = 3, 4$ using consistent Brownian increments, that is, the same Brownian motion for each R. We clearly see that the Richardson-Romberg estimator efficiency increases with R and the method gives satisfying results with $R = 3, 4$ for small values of n.

<figure><img src="image_5.png" /><figcaption>Figure 4.3: Richardson Romberg SA estimators: $\sum_{r=1}^R \mathbf{w}_r \theta^{*,rn} - \theta^*$ with respect to $n = 2, \dots, 15$ for $R = 2, 3, 4$.</figcaption></figure>

In order to illustrate the computational efficiency of the Richardson-Romberg estimator we set the optimal param-
<table><caption>Table 4.1: Richardson-Romberg SA estimators for the quantile at level $\ell$ of a geometric Brownian motion with a target accuracy $\varepsilon = 2^{-p}$, $p = 1, \dots, 4$.</caption><thead><tr><th>Target accuracy: ε</th><th>L<sup>1</sup>-error</th><th>time (s)</th><th>R</th><th>n</th><th>M</th></tr></thead><tbody><tr><td>5.00 × 10<sup>-1</sup></td><td>3.21 × 10<sup>-1</sup></td><td>0.9 × 10<sup>1</sup></td><td>2</td><td>14</td><td>8.69 × 10<sup>5</sup></td></tr><tr><td>2.50 × 10<sup>-1</sup></td><td>4.80 × 10<sup>-2</sup></td><td>5.15 × 10<sup>1</sup></td><td>2</td><td>20</td><td>3.48 × 10<sup>6</sup></td></tr><tr><td>1.25 × 10<sup>-1</sup></td><td>4.32 × 10<sup>-2</sup></td><td>1.70 × 10<sup>2</sup></td><td>3</td><td>8</td><td>1.21 × 10<sup>7</sup></td></tr><tr><td>6.25 × 10<sup>-2</sup></td><td>3.48 × 10<sup>-2</sup></td><td>7.92 × 10<sup>2</sup></td><td>3</td><td>10</td><td>4.85 × 10<sup>7</sup></td></tr></tbody></table>

<table><caption>Table 4.2: Crude SA estimators for the quantile at level $\ell$ of a geometric Brownian motion with a target accuracy $\varepsilon = 2^{-p}$, $p = 1, \dots, 4$.</caption><thead><tr><th>Target accuracy: ε</th><th>L<sup>1</sup>-error</th><th>time (s)</th><th>n</th><th>M</th></tr></thead><tbody><tr><td>5.00 × 10<sup>-1</sup></td><td>2.09 × 10<sup>-1</sup></td><td>1.09 × 10<sup>2</sup></td><td>235</td><td>1.25 × 10<sup>6</sup></td></tr><tr><td>2.50 × 10<sup>-1</sup></td><td>3.84 × 10<sup>-2</sup></td><td>8.18 × 10<sup>2</sup></td><td>469</td><td>5.01 × 10<sup>6</sup></td></tr><tr><td>1.25 × 10<sup>-1</sup></td><td>3.48 × 10<sup>-2</sup></td><td>7.09 × 10<sup>3</sup></td><td>938</td><td>2.00 × 10<sup>7</sup></td></tr><tr><td>6.25 × 10<sup>-2</sup></td><td>2.91 × 10<sup>-2</sup></td><td>5.25 × 10<sup>4</sup></td><td>1876</td><td>8.01 × 10<sup>7</sup></td></tr></tbody></table>

eters according to (4.8), namely

$$n(\varepsilon) = \left[ \left( \frac{2\alpha R}{\beta} + 1 \right)^{\frac{1}{\alpha R}} \mu_R^{\frac{1}{\alpha R}} \varepsilon^{-\frac{1}{\alpha R}} \right] \quad \text{and} \quad M(\varepsilon) = \left[ \gamma_0^{\frac{1}{\beta}} \nu_R^{\frac{2}{\beta}} \left( 1 + \frac{\beta}{2\alpha R} \right)^{\frac{2}{\beta}} \varepsilon^{-\frac{2}{\beta}} \right].$$

The target accuracy $\varepsilon$ for the L$^1$-error has been set at $\varepsilon = 2^{-p}$, $p = 1, \dots, 4$. The L$^1$-error is estimated using 400 runs of the algorithm. The results are summarized in Table 4.1 for the Richardson-Romberg extrapolation SA method and in Table 4.2 for the crude SA method.<fnref n="1" /> Note that as expected the L$^1$-error is always lower than the specified $\varepsilon$ for our estimators. Using the Richardson-Romberg SA scheme instead of the crude SA method leads to a gain in terms of CPU-time varying from 12 (for $\varepsilon = 5.00 \times 10^{-1}$) to 66 (for $\varepsilon = 6.25 \times 10^{-2}$).

## 4.4 Perspective

From a theoretical point of view it would be interesting to obtain a CLT for the multi-level Robbins-Monro algorithm when one uses the averaging principle of Ruppert-Polyak since we only obtain the optimal computational cost under a constraint on $\gamma_0$ when $\gamma(p) = \gamma_0/p$. One expects to remove this constraint by devising a Robbins-Monro algorithm devised with a low decreasing step sequence $\gamma(p) = \gamma_0/p^a$, $\frac{1}{2} < a < 1$ and computing the empirical mean along its trajectories.

From a practical point of view, we investigated multi-level and multi-step Richardson-Romberg methods. It would be interesting to combine them in order to obtain a multi-level Richardson-Romberg method for stochastic approximation as it is done by Pagès and Lemaire [93] for the estimation of $\mathbb{E}[f(U)]$ by the Monte Carlo method. One expects to achieve a complexity of order $O(n^2 \log(n))$ with a target accuracy of order $O(n^{-1})$, that is, to remove a $\log(n)$ factor in the multi-level optimal complexity.

<footnote n="1">¹The computations were performed on a computer with 4 multithreaded(16) octo-core processors (Intel(R) Xeon(R) CPU E5-4620 @ 2.20GHz).</footnote>
# Part II : Markovian perturbation of stochastic processes and applications
# Chapter 5

## Weak approximation and asymptotic expansion of some Markov semigroups

In this chapter we present [11] and [12]. The first work, in collaboration with A. Kohatsu-Higa, deals with asymptotic expansions of Markov semigroups with respect to a small parameter. The key tools are the parametrix method and Malliavin's calculus. We revisit heat kernel expansions of some hypo-elliptic diffusion processes. We also discuss some extensions to skew diffusions with non-regular coefficients.

In the second work, we are interested in the weak approximation error of a skew diffusion with bounded measurable drift and Hölder continuous diffusion coefficient by an Euler-type scheme which consists of iteratively simulating skew Brownian motions with constant drift. Here, we rely on a discrete version of the parametrix method “à la McKean and Singer” [95]. This approach has been successfully applied in the Brownian case by Konakov and Mammen [75] in order to derive some local limit theorems for Markov chain approximations and to establish weak approximation expansions for the Euler scheme of a diffusion [76] when the coefficients of the underlying scheme are smooth. We also refer to Konakov and Menozzi [77] for an extension to stable driven SDEs still in the case of smooth coefficients. However, little attention has so far been given to the case of SDEs with non-smooth coefficients.

We first briefly present in Section 5.1.1 a probabilistic point of view of the parametrix method “à la McKean and Singer”. For simplicity of the presentation, we consider the case of time-homogeneous diffusion process. In Section 5.1.2, we present a discrete version for the approximation of Lévy driven SDEs.

In Section 5.2, we illustrate the effectiveness of the parametrix method for the weak approximation of a skew diffusion with non-smooth coefficients. In Section 5.3, we present a general approach for asymptotic expansion of a Markov semigroup $(P_t^\epsilon)_{t \ge 0}$ with respect to the small parameter $\epsilon$ by means of the parametrix method and Malliavin's calculus and illustrate it on two examples.

### 5.1 A brief presentation of the parametrix method “à la McKean and Singer”

#### 5.1.1 The case of diffusion processes

Let $T > 0$. We consider the (strong) time-homogeneous Markov process $(X_t)_{t \in [0,T]}$ taking values in $\mathbb{R}^d$ with infinitesimal generator defined by:

$$ \mathcal{L}f(x) = \langle b(x), \nabla f(x) \rangle + \frac{1}{2} \mathrm{Tr}(a(x)D_x^2 f(x)), \quad x \in \mathbb{R}^d \qquad (5.1) $$

where $b: \mathbb{R}^d \to \mathbb{R}^d$ and $\sigma: \mathbb{R}^d \to \mathbb{R}^d \otimes \mathbb{R}^q$ are bounded measurable functions, $a(\cdot) = (\sigma\sigma^*)(\cdot)$ is $\eta$-Hölder continuous and uniformly elliptic, that is, $\forall (x, \xi) \in (\mathbb{R}^d)^2$, $\langle a(x)\xi, \xi \rangle \ge a|\xi|^2$. Under these assumptions, it is known that the related martingale problem is *well-posed*, see Stroock and Varadhan [120], in other words, the SDE with $\mathcal{L}$ as infinitesimal generator admits a unique weak solution. We denote by $(P_t)_{t \ge 0}$ the Markov semigroup defined on bounded measurable functions $f: \mathbb{R}^d \to \mathbb{R}$ by $P_t f(x) = \mathbb{E}[f(X_t)|X_0 = x]$ where $(X_t)_{t \ge 0}$ is the unique weak solution to the SDE with $\mathcal{L}$ as infinitesimal generator. The idea of the parametrix method is to approximate the dynamics of $X$ by a simple process $\bar{X}^z$, referred as parametrix process, obtained by removing the drift and freezing the diffusion coefficient to a fixed point $z \in \mathbb{R}^d$. More precisely, one has $\bar{X}_t^z = x + \sigma(z)W_t$, where $(W_t)_{t \in [0,T]}$ stands for a $q$-dimensional Brownian motion, with infinitesimal generator $\bar{\mathcal{L}}^z f(x) = \frac{1}{2} \mathrm{Tr}(a(z)D_x^2 f(x))$. We denote by $(\bar{P}_t^z)_{0 \le t \le T}$ the associated Markov semigroup. For $t > 0$, the (transition) density of the parametrix process taken at time $t$ starting from $x$ at time 0 is $y \mapsto \bar{p}^z(t, y - x)$, where $\bar{p}^z(t, x)$ is the Gaussian density with mean 0 and covariance matrix $a(z)t$ taken at $x$. We also
denote by $(\delta_{\varepsilon}^z)_{\varepsilon>0}$ a sequence of positive mollifier converging weakly to the Dirac mass at $z$. For instance, one may consider $\delta_{\varepsilon}^z(y) = g(\varepsilon, y-z)$ with the notation $g(c, x) := (2\pi c)^{-\frac{d}{2}} \exp(-|y-z|^2/(2c))$. Let $r, \varepsilon > 0$ and $f \in C_b(\mathbb{R}^d)$. By Itô's lemma, for $s \in [0, t]$, one has

$$
\begin{align}
\bar{P}_{t-s}^z \bar{P}_r^z \delta_{\varepsilon}^z f(X_s) &= \bar{P}_{t+r}^z \delta_{\varepsilon}^z f(x) + \int_0^s (\partial_v + \mathcal{L}) \bar{P}_{t-v+r}^z \delta_{\varepsilon}^z f(X_v) dv + M_s \nonumber \\
&= \bar{P}_{t+r}^z \delta_{\varepsilon}^z f(x) + \int_0^s (\mathcal{L} - \bar{\mathcal{L}}^z) \bar{P}_{t-v+r}^z \delta_{\varepsilon}^z f(X_v) dv + M_s \tag{5.2}
\end{align}
$$

where $(M_s)_{s \in [0,t]}$ is a square integrable martingale. Taking expectations in both sides of the previous equality and letting $s \to t$ and then integrating with respect to $z$, one gets

$$
\int P_t \bar{P}_r^z \delta_{\varepsilon}^z f(x) dz = \int \bar{P}_{t+r}^z \delta_{\varepsilon}^z f(x) dz + \int_0^t \mathbb{E}[\int \int \delta_{\varepsilon}^z(y) f(y) S_{t-v+r}^z(X_v, y) dz dy] dv. \quad (5.3)
$$

where we introduced the notations:

$$
S_t^z f(x) := \int f(y)S_t^z(x, y)dy, \quad S_t^z(x, y) := \left\{ \frac{1}{2}(a_{i,j}(x) - a_{i,j}(z))H_2^{i,j}(a(z)t, y-x) - b_i(x)H_1^i(a(z)t, y-x) \right\} g(a(z)t, y-x)
$$

with $H_2^{i,j}(a,x) = (g(a,x))^{-1} \partial_{x_i x_j} g(a,x)$, $H_1^i(a,x) = (g(a,x))^{-1} \partial_{x_i} g(a,x)$. Passing to the limit as $\varepsilon \downarrow 0$ then in $r \downarrow 0$
and omitting technical details, we get

$$
P_t f(x) = \hat{P}_t f(x) + \int_0^t P_v \hat{S}_{t-v} f(x) dv, \quad t \in [0, T] \quad (5.4)
$$

with $\hat{P}_t f(x) := \int f(y)g(a(y)t, y-x)dy$ and $\hat{S}_t f(x) := \int f(y)\hat{S}_t(x, y)dy$, $\hat{S}_t(x, y) := \hat{S}_t^y(x, y)$. Importantly, using the
$\eta$-Hölder regularity of $a$, one has the following control: $\forall t \in (0, T]$, $|\hat{S}_t(x, y)| \le Ct^{-1+\frac{\eta}{2}}g(ct, y-x)$ for some positive
constants $C := C(a, b, T)$, $c := c(a)$.

The above lines of reasoning allows to select the point $z$ where the coefficients of the parametrix process are fixed
as the terminal point where the transition density is evaluated. This argument will be referred as the *diagonalisation argument*.
In Section 6.2 of Chapter 6, we will see that this argument still works when the probability measure of the
parametrix process is not absolutely continuous with respect to the Lebesgue measure. Assuming that $b$ is continuous,
one may iterate the first step expansion (5.4):

$$
P_t f(x) = \sum_{n=0}^{N-1} I_t^n f(x) + \mathcal{R}_t^N f(x) \quad (5.5)
$$

where $I_t^n f(x) := \int_{\Delta_n(t)} \hat{P}_{s_n} \hat{S}_{s_{n-1}-s_n} \cdots \hat{S}_{t-s_1} f(x) ds_n$, $I_t^0 f(x) := \hat{P}_t f(x)$ and $\mathcal{R}_t^N f(x) := \int_{\Delta_N(t)} P_{s_N} \hat{S}_{s_{N-1}-s_N} \cdots \hat{S}_{t-s_1} f(x) ds_N$.
Using repeatedly the previous control on the kernel $\hat{S}_t(x, y)$, one gets

$$
|\mathcal{R}_t^N f(x)| \le |f|_\infty \int_{\Delta_N(t)} ds_N \prod_{n=0}^{N-1} C_T(s_n - s_{n+1})^{-(1-\frac{\eta}{2})} = |f|_\infty C_T^N t^{N\eta/2} \frac{\Gamma(\eta/2)^N}{\Gamma(1+N\eta/2)},
$$

where $\Gamma$ denotes the Gamma function. From the asymptotics of the Gamma function, one may let $N$ goes to infinity
in (5.5) and deduces that the series $\sum_{n \ge 0} I_t^n f(x)$ converges absolutely and uniformly for $t \in [0, T]$. From the previous
expansion, one may prove the existence of the transition density for the process $(X_t)_{t \in [0,T]}$. Indeed, from Fubini's
theorem, one has

$$
P_t f(x) = \sum_{n \ge 0} I_t^n f(x) = \int f(z_0) \sum_{n \ge 0} \left( \int_{\Delta_n(T)} ds_n \int_{(\mathbb{R}^d)^n} dz_n g(a(z_n)s_n, z_n - x) \hat{S}_{t-s_1}^{z_0}(z_1, z_0) \prod_{i=1}^{n-1} \hat{S}_{s_i-s_{i+1}}^{z_i}(z_{i+1}, z_i) \right) dz_0
$$

so that the transition density $(t, x, y) \mapsto p(t, x, y)$ associated to $(P_t)_{t \ge 0}$ exists and is given by

$$
p(t, x, y) = \sum_{n \ge 0} \int_{\Delta_n(T)} ds_n \int_{(\mathbb{R}^d)^n} dz_n g(a(z_n)s_n, z_n - x) \hat{S}_{t-s_1}^{z_0}(z_1, y) \prod_{i=1}^{n-1} \hat{S}_{s_i-s_{i+1}}^{z_i}(z_{i+1}, z_i). \quad (5.6)
$$

Moreover, from the asymptotics of the Gamma function, one deduces the Gaussian upper-bound: ∃C := C(a, |b|∞, T), c :=
c(a) > 0 such that ∀(t, x, y) ∈ (0, T] × (ℝd)2, p(t, x, y) ≤ Cg(ct, y − x). The case of a bounded measurable drift can be
handled by an approximation argument: take (b_n)_{n≥0} such that b_n is continuous, sup_n≥1 |b_n|_∞ ≤ |b|_∞ and b_n → b a.e.
Denote by (P_n^t)_{t≥0} the semigroup associated to this approximation. By a compactness argument and uniqueness of
the martingale problem, { (X_n^t)_{t∈[0,T]}, n ≥ 0 } converges weakly to (X_t)_{t≥0} and one can pass to the limit in the corre-
sponding parametrix series. When b is Hölder continuous, one can prove that the series expansion for (t, x) ↦→ p(t, x, y)
satisfies the Backward Kolmogorov equation, namely: (∂_t + L)p(t, x, y) = 0 with p(t, ., y) → δ_y(.) as t ↓ 0, see e.g.
Friedman [44].
5.1.2 A discrete version of the parametrix method

In the previous section, we have exposed a probabilistic point of view of the parametrix method of McKean and Singer [95] for a time-homogeneous diffusion process. We illustrate now its robustness by showing that it can be easily adapted for approximation schemes. In this section, we consider an $\mathbb{R}^d$-valued process $(X_t)_{t \in [0,T]}$ with generator $\mathcal{L}$ given by (5.1).

For a positive integer $N$, we define the time step $h = T/N$ and the corresponding uniform time grid $t_i = ih$, $i = 0, \dots, N$. A natural approximation scheme of the dynamics of $(X_t)_{t \in [0,T]}$ by a Markov chain is given by the following scheme:

$$X_{t_{i+1}}^N = X_{t_i}^N + b(X_{t_i}^N)h + \sigma(x')\zeta_{i+1}^N, \quad i = 0, \dots, N-1, \quad X_0^N = x \quad (5.7)$$

where $(\zeta_i^N)_{1 \le i \le N}$ are i.i.d. random variables such that their common law approximate the law of the true increments $W_h$ and such that $\zeta_1^N$ has a density. If $\sigma$ is non-degenerated then the one step transition admits a density and so does $X_{t_i}^N$ for $i = 1, \dots, N$. The density of the Markov chain approximation scheme taken at time $t_{j'}$ starting from $x$ at time $t_j$ will be denoted $y \mapsto p_N(t_j, t_{j'}, x, y)$. In the spirit of the previous section, we introduce a parametrix process denoted by $(\tilde{X}_{t_i}^{N,x'})_{j \le i \le j'}$ associated to $(X_{t_i}^N)_{j \le i \le j'}$ which consists in removing the drift part and freezing the diffusion coefficient in the dynamics of (5.7) at the terminal point $x'$ where the density of $(X_{t_i}^N)_{j \le i \le j'}$ is evaluated. Its dynamics is given by

$$\tilde{X}_{t_{i+1}}^{N,x'} = \tilde{X}_{t_i}^{N,x'} + \sigma(\tilde{X}_{t_i}^{N,x'})\zeta_{i+1}^N, \quad i = j, \dots, j' - 1, \quad \tilde{X}_j^{N,x'} = x \quad (5.8)$$

and its transition density is denoted by $\tilde{p}_N^{x'}$. To simplify the notations, we will write $\tilde{p}_N(t_j, t_{j'}, x, x')$ for $\tilde{p}_N^{x'}(t_j, t_{j'}, x, x')$. We also introduce the discrete counterpart of the infinitesimal generators considered so far. For $g: C_b(\mathbb{R}^d)$, $j = 0, \dots, j' - 1$, we define

$$\mathcal{L}_{t_j}^N g(x) = \frac{\mathbb{E}[g(X_{t_j+h}^N) | X_{t_j}^N = x] - g(x)}{h}, \quad \text{and} \quad \tilde{\mathcal{L}}_{t_j}^N g(x) = \frac{\mathbb{E}[g(\tilde{X}_{t_j+h}^{N,x'}) | \tilde{X}_{t_j}^N = x] - g(x)}{h} \quad (5.9)$$

for $0 \le j < j' \le N$, the discrete kernel $H_N(t_j, t_{j'}, x, x') = ((\mathcal{L}_{t_j}^N - \tilde{\mathcal{L}}_{t_j}^N) \tilde{p}_N(t_j + h, t_{j'}, \dots, x'))(x)$ and finally the discrete time-space convolution type operator $\otimes_N$ as follows

$$(g \otimes_N f)(t_j, t_{j'}, x, x') = \sum_{i=j}^{j'-1} h \int_{\mathbb{R}^d} g(t_j, t_i, x, z) f(t_i, t_{j'}, z, x') dz$$

with the convention that $\sum_{i=j}^{j'-1} \dots = 0$ if $j \ge j'$. From (5.9), it is easily seen that

$$H_N(t_j, t_{j'}, x, x') = h^{-1} \int_{\mathbb{R}^d} (p_N - \tilde{p}_N^{t_{j'}, x'}) (t_j, t_{j+1}, x, z) \tilde{p}_N(t_{j+1}, t_{j'}, z, x') dz.$$

After some simple computations that we omit, see e.g. Lemma 3.6 in Konakov and Mammen [75], for $0 \le j < j' \le N$, the following expansion holds

$$\forall (x, x') \in (\mathbb{R}^d)^2, p_N(t_j, t_{j'}, x, x') = \sum_{r=0}^{j'-j} (\tilde{p}_N \otimes_N H_N^{(r)}) (t_j, t_{j'}, x, x') \quad (5.10)$$

where we define the convolution as follows: $g \otimes_N H_N^{(0)} = g$ and for $r \ge 1$, $g \otimes_N H_N^{(r)} = (g \otimes_N H_N^{(r-1)}) \otimes_N H_N$ and use the convention $\tilde{p}(t_{j'}, t_{j'}, \dots, x') = \tilde{p}(t_{j'}, t_{j'}, \dots, x') = \delta_{x'}(\cdot)$ in the computation of $\tilde{p} \otimes_N H_N^{(r)}$. From (5.10) and (5.6), the densities of $X_t$ and of its approximation scheme share a similar representation. In order to control the weak approximation error $p - p_N$ between the two densities, one has to quantify the following errors:

1. the error coming from the approximation of the “true” Brownian increment by the sequence $(\zeta_i^N)_{1 \le i \le N}$, in other words the error $\tilde{p} - \tilde{p}_N$.

2. the error coming from the discretisation of the space-time convolution operator appearing in (5.6).

3. the contribution when one replaces the continuous kernel $\hat{S}_{t-s}(\cdot, \cdot)$ by its discrete counterpart $H_N(s, t, \cdot, \cdot)$.

When the coefficients are smooth and time-inhomogeneous, following these steps, Konakov and Mammen [75] proved an error bound for the difference $p - p_N$, namely $|(p - p_N)(0, t, x, x')| \le C h^{1/2} t^{-d/2} (1 + (|x' - x|/t^{1/2})^S)^{-1}$ for some $S$ depending of the integrability properties of the underlying innovation sequence $(\zeta_i^N)_{1 \le i \le N}$. It is the approximation coming from step 1 that dominates in the error bound.
When one directly samples the sequence $(\zeta_i^N)_{1 \le i \le N}$ according to the true law $h^{\frac{1}{2}}\mathcal{N}(0, I_d)$, that is, when one removes the error coming from step 1, one obtains an error bound of order $h$. An expansion of the error in powers of $h$ has been established by Konakov and Mammen [76], as mentioned in equation (1.30) in Section 1.2.3. All these results are obtained under enough regularity on the coefficients $b$ and $\sigma$.

When the coefficients $b$ and $\sigma$ are not smooth, there are only few results. Very recently, in the diffusion case, Konakov and al. [78] derived an upper bound for $p - p_N$ of order $h^{\frac{\eta}{2} - C\psi(h)}$, where $\psi(h)$ is a slowly varying factor that goes to zero as $h \to 0$ under the assumption that the coefficients $b$ and $\sigma$ are $\eta$-Hölder continuous. Their strategy consists in introducing perturbed dynamics of the considered SDE and its scheme by suitably mollifying the coefficients of both dynamics and to quantify the distance between the densities and their respective perturbations.

We also mention the work of Mikulevicius and Platen [101] who established an error bound for the weak approximation error $E[f(X_T) - f(X_T^N)]$ of order $h^{\eta/2}$ provided $f \in C^{2+\eta}([0, T] \times \mathbb{R}^d)$ and the coefficients $b$ and $\sigma$ are $\eta$-Hölder continuous in space and $\eta/2$-Hölder continuous in time. We also refer the reader to [100] and [102] for some recent extensions of this result to the case of Lévy driven SDEs.

In the next section we present the results obtained in [12] where we proved an error bound for the difference between the densities of the skew diffusion and its Euler approximation of order $h^{\eta/2}$ under the assumption that $b$ is bounded measurable and $\sigma$ is $\eta$-Hölder continuous and uniformly elliptic.

## 5.2 Weak approximation of a skew diffusion by an Euler-type scheme

We consider the unique weak solution of the following $\mathbb{R}$-valued SDE $(X_t)_{t \ge 0}$ with dynamics

$$X_t = x + \int_0^t b(X_s)ds + \int_0^t \sigma(X_s)dW_s + (2\alpha - 1)L_t^0(X), \quad \alpha \in (0, 1) \quad (5.11)$$

where $W = (W_t)_{t \ge 0}$ is a one dimensional Brownian motion defined on a filtered probability space $(\Omega, \mathcal{F}, (\mathcal{F}_t)_{t \ge 0}, \mathbb{P})$ satisfying the usual assumptions and $L^0(X)$ is the symmetric local time of $X$ at the origin.

When $b = 0$ and $\sigma = 1$, the solution to (5.11) is called the skew Brownian motion. Harrison and Shepp [60] proved that if $|2\alpha - 1| \le 1$ then there is a unique strong solution and if $|2\alpha - 1| > 1$, there is no solution. The case $\alpha = 1$ corresponds to reflecting Brownian motion.

Here we assume that the coefficients of the SDE (5.11) satisfies the following assumptions:

(HR) The drift $b$ is bounded measurable and $a = \sigma^2$ is $\eta$-Hölder continuous, for some $\eta \in (0, 1]$. That is, there exists a positive constant $L$ such that

$$\sup_{x \in \mathbb{R}} |b(x)| + \sup_{(x,y) \in \mathbb{R}^2, x \neq y} \frac{|a(x) - a(y)|}{|x - y|^{\eta}} < L.$$

(HE) The diffusion coefficient is uniformly elliptic that is there exists $\lambda > 1$ such that for every $x \in \mathbb{R}$, $\lambda^{-1} < a(x) < \lambda$. Since $\sigma$ is continuous, without loss of generality, we may assume that $\sigma$ is positive.

The previous assumptions guarantee the existence of a unique weak solution to (5.11). Moreover, for any $(t, x) \in \mathbb{R}_+^* \times \mathbb{R}$, $X_t$ admits a density $y \mapsto p(0, t, x, y)$, which is continuous on $\mathbb{R}^*$ and satisfies a Gaussian upper-bound. We refer to Kohatsu-Higa and al. [73] for more details, see also Le Gall [88]. We also refer the interested reader to the recent survey of Lejay [89] and the references therein for various applications of such equation.

In particular, in [73], a representation of the transition density $p(0, t, x, y)$ by means of the parametrix method “à la McKean and Singer” similar to the expansion (5.6) is obtained. However, the parametrix process is chosen to be the skew Brownian motion, that is, $\bar{X}_t^z = x + \sigma(z)W_t + (2\alpha - 1)L_t^0(\bar{X}^z)$ for which the transition density $\bar{p}^z(0, t, x, y)$ is defined by

Case 1: For $x \ge 0$, one has

$$\bar{p}^z(0, t, x, y) = \{g(a(z)t, y - x) + (2\alpha - 1)g(a(z)t, y + x)\} \mathbf{I}_{\{y \ge 0\}} + 2(1 - \alpha)g(a(z)t, y - x) \mathbf{I}_{\{y < 0\}}.$$

Case 2: For $x < 0$, one has

$$\bar{p}^z(0, t, x, y) = \{g(a(z)t, y - x) + (1 - 2\alpha)g(a(z)t, y + x)\} \mathbf{I}_{\{y < 0\}} + 2\alpha g(a(z)t, y - x) \mathbf{I}_{\{y \ge 0\}}.$$

Note also that the infinitesimal generators $\mathcal{L}$ of (5.11) and $\bar{\mathcal{L}}^z$ of $\bar{X}^z$ are respectively given by $\mathcal{L} = \frac{1}{2}a(x)\partial_x^2 + b(x)\partial_x$ and $\bar{\mathcal{L}}^z f(x) = \frac{1}{2}a(z)\partial_z^2$ acting on $\mathcal{D}(\alpha)$ the set of bounded continuous functions $f: \mathbb{R} \to \mathbb{R}$ with bounded continuous derivatives $f'$ and $f''$ on $\mathbb{R}^*$ such that $f'(0+)$ and $f'(0-)$ exists, are finite and satisfies $\alpha f'(0+) = (1 - \alpha)f'(0-)$. Instead of Itô's lemma, one has to use the Itô-Tanaka formula in (5.2).
To approximate equation (5.11) on the time interval $[0, T]$, $T > 0$, we introduce the Euler type scheme $(X_t^N)_{t \in [0,T]}$ associated to the time step $h = T/N$, $N \in \mathbb{N}^*$ and time grid $t_i = ih$, $i = 0, \dots, N$, defined by $X_0^N = x$ and for all $t_i \le t \le t_{i+1}$

$$X_t^N = X_{t_i}^N + b(X_{t_i}^N)(t - t_i) + \sigma(X_{t_i}^N)(W_t - W_{t_i}) + (2\alpha - 1)L_{t-t_i}^0(X^N(t_i, X_{t_i}^N)). \quad (5.12)$$

Observe that the above scheme does not correspond to a standard Euler-Maruyama approximation scheme since we do not discretise the part corresponding to the local time in the dynamics (5.11). However, its computation only requires to be able to simulate exactly the skew Brownian motion with a constant drift at time $t - t_i$. This process is known to be exactly simulatable and we refer to Etoré and Martinez [40] for the exact expression of its density. In particular, under assumptions (HR) and (HE), the transition between $t_i$ and $t \in (t_i, t_{i+1}]$ admits a positive density so that the discretisation scheme (5.12) admits a positive transition density between two times $t_i$ and $t_j$, that we will denote $p_N(t_i, t_j, x, y)$, for any $0 \le t_i < t_j \le T$, $(x, y) \in \mathbb{R} \times \mathbb{R}^*$. In particular, a Gaussian upper-bound has been established in [73] for the transition density of the skew diffusion (5.11) under (HR) and (HE). Our first result is to prove similar Aronson's estimate, that is a Gaussian upper estimate but also a lower bound hold for the discretisation scheme (5.12).

**Theorem 5.2.1.** (Two sided Gaussian estimates for the scheme) Under (HR) and (HE), there exist two constants $C := C(T, b, \sigma)$, $c := c(\lambda, \eta) > 1$ such that for every $0 \le j < i \le N$,

$$\forall (x, y) \in \mathbb{R} \times \mathbb{R}^*, \quad C^{-1}g(c^{-1}(t_i - t_j), y - x) \le p_N(t_j, t_i, x, y) \le Cg(c(t_i - t_j), y - x).$$

The proof relies on a parametrix expansion of the density $p_N(t_j, t_i, x, y)$ as explained in the previous sections. In particular, one proves that the series representation (5.10) holds for the transition density of the scheme (5.12) with $\tilde{p}_N(t_j, t_{j'}, x, x') = \tilde{p}(t_j, t_{j'}, x, x') = \bar{p}^{x'}(t_j, t_{j'}, x, x')$. The lower bound is obtained by a chaining argument as usually done in this context, see Chapter VII in Bass [14], Kusuoka and Stroock [84] and Lemaire and Menozzi [91] in the case of discretisation schemes of uniformly elliptic diffusions and of some degenerate Kolmogorov processes.

Our second result is an error bound for the difference between the densities of the skew diffusion (5.11) and its approximation scheme (5.12).

**Theorem 5.2.2** (Error bound on the difference between the densities). Assume that (HR) and (HE) hold. Then, there exists a constant $c := c(\lambda, \eta) > 1$ such that, for all $0 \le t_j < t_i \le N$, one has

$$\forall (x, y) \in \mathbb{R} \times \mathbb{R}^*, \quad |(p - p_N)(t_j, t_i, x, y)| \le C(T, b, \sigma)h^{\eta/2}g(c(t_i - t_j), y - x)$$

where $T \mapsto C(T, b, \sigma)$ is a non-decreasing positive function.

Observe that the weak rate $h^{\eta/2}$ is coherent with previous results obtained in the literature for the weak approximation error $\mathbb{E}[f(X_T) - f(X_T^N)]$ where $X$ is given by a Brownian SDE and $X^N$ is its Euler scheme when $b$ and $\sigma$ are $\eta$-Hölder continuous in space and $\eta/2$-Hölder continuous in time, see [101]. In Konakov and al. [78], the error bound on the densities is proved to be of order $h^{\eta/2 - C\psi(h)}$, $\psi(h) = \log_3(h^{-1})/\log_2(h^{-1})$ when the coefficients $b$ and $\sigma$ are $\eta$-Hölder continuous. Moreover, in their results, the singularity in time is given by $C(t_i - t_j)^{-(1-\eta/2)\eta/2}$ whereas this quantity does not appear in the previous theorem so it is tighter in this sense. For $\alpha = 1/2$, which corresponds to the case of (time homogeneous) diffusion process (since the local time part vanishes), compared to [78], our result notably removes the slowly varying factor $\psi(h)$ and shows that the drift plays no role in the approximation of the density since we only require $b$ to be a bounded measurable function. This phenomenon is not surprising since one of the advantages of the parametrix method is the removal of the drift part in the analysis. Eventually, it should be possible to extend our strategy of proof to the case of multi-dimensional Brownian diffusion.

In order to prove the error bound of Theorem 5.2.2, our strategy is the following. The main point is to compare the parametrix series of the transition $p(t_j, t_i, x, y)$ obtained in [73] and the parametrix series of the Euler scheme (5.12) which differs on account of the discrete nature of time-space convolution operator $\otimes_N$ and the discrete smoothing kernel $H_N$. In order to do this, we introduce for $0 \le j < j' \le N$,

$$\forall (x, x') \in \mathbb{R} \times \mathbb{R}^*, \quad p_N^d(t_j, t_{j'}, x, x') = \sum_{r \ge 0} \tilde{p} \otimes_N H^{(r)}(t_j, t_{j'}, x, x'). \quad (5.13)$$

Arguments similar to those used in Section 5.1.1 show that the series (5.13) converge absolutely and uniformly on $\mathbb{R} \times \mathbb{R}^*$ and that $p_N^d$ satisfies the following Gaussian upper-bound:

$$\forall (x, x') \in \mathbb{R} \times \mathbb{R}^*, \quad p_N^d(t_j, t_{j'}, x, x') \le C(T, |b|_\infty, a)g(c(t_{j'} - t_j), x' - x). \quad (5.14)$$

Indeed, by induction on $r$, one proves that

$$\forall (x, x') \in \mathbb{R} \times \mathbb{R}^*, \quad |\tilde{p} \otimes_N H^{(r)}(t_j, t_{j'}, x, x')| \le C^r (t_{j'} - t_j)^{r\eta/2} \prod_{i=1}^{r} B \left( 1 + \frac{(i-1)\eta}{2}, \frac{\eta}{2} \right) g_{c(t_{j'}-t_j)}(x' - x) \quad (5.15)$$
where $C := C(\lambda, \eta)(|b|_\infty T^{\frac{1-n}{2}} + 1)$ and $B(m, n) = \int_0^1 dv (1-v)^{m-1}v^{n-1}$ is the Beta function. Then, (5.14) follows from the previous bound. We omit technical details.

The idea is now to decompose the global error as follows:

$$ (p - p_N)(t_j, t_i, x, y) = (p - p_N^d)(t_j, t_i, x, y) + (p_N^d - p_N)(t_j, t_i, x, y). $$

A similar decomposition has been used in [76] when the coefficients $b$ and $\sigma$ are smooth. The smoothness of the coefficients notably allows to use Taylor expansions to express the transition density $p$ as the fundamental solution of the underlying parabolic PDE and to use integration by parts (that may be expressed as the duality relation satisfied by the infinitesimal generator when seen as a differential operator) in order to equilibrate time singularities. Obviously, these arguments do not work here under the mild smoothness assumption (HR) so that computations become more delicate. In a first step, we express the difference $p_N^d - p_N$ in an infinite parametrix series that involves the difference between the two kernels $H$ and $H_N$. Then, the symmetric Itô-Tanaka formula allows to express the difference $H - H_N$ as the difference of the kernel $H$ between two consecutive discretisation times plus a remainder term $\mathcal{R}_N$. We then study the weak approximation rate induced by $\mathcal{R}_N$. The delicate point here lies in the non-differentiability of $x \mapsto \partial_x \tilde{p}(t_j, t_j', x, x')$ at zero caused by the presence of the local time part in the dynamics of $X^N$ and $\tilde{X}^N$ which prevents us to use (again) the Itô-Tanaka formula. In a second step, in order to study $p - p_N^d$, we express this difference as an infinite parametrix series that involves the difference between the two convolution operators $\otimes$ and $\otimes_N$. Then, we notably use a (kind of) Lipschitz property in time for the transition density $p$ and a smoothing procedure for the drift part.

## 5.3 Asymptotic expansion of Markov semigroups and application to diffusion processes

In this section, we present the results obtained in [11] in collaboration with A. Kohatsu-Higa. We develop a general framework for asymptotic expansion of a Markov semigroup $(P_t^\varepsilon)_{t \ge 0}$ with respect to a small parameter $\varepsilon$ based on the parametrix technique and Malliavin's calculus. The main idea can be summed up as follows. Assume that a first order expansion for the infinitesimal generator $\mathcal{L}_t^\varepsilon$ of $(P_t^\varepsilon)_{t \ge 0}$ holds, that is $\mathcal{L}_t^\varepsilon = \mathcal{L}_t^0 + \varepsilon \mathcal{L}_t^{1, \varepsilon}$, $t \ge 0$, where $\mathcal{L}_t^0$ is the infinitesimal generator of some strongly continuous Feller semigroup $(P_{s,t}^0)_{0 \le s \le t}$, which will be referred to as the proxy semigroup. The basic argument in order to obtain the expansion in powers of $\varepsilon$ consists in writing

$$ P_{0,t}^{\varepsilon} f(x) - P_{0,t}^{0} f(x) = \int_{0}^{t} \partial_{s} (P_{0,s}^{\varepsilon} P_{s,t}^{0}) f(x) ds = \int_{0}^{t} P_{0,s}^{\varepsilon} (\mathcal{L}_{s}^{\varepsilon} - \mathcal{L}_{s}^{0}) P_{s,t}^{0} f(x) ds = \varepsilon \int_{0}^{t} P_{0,s}^{\varepsilon} \mathcal{L}_{s}^{1, \varepsilon} P_{s,t}^{0} f(x) ds \quad (5.16) $$

where we assumed that $P_{s,t}^0 f \in \text{Dom}(\mathcal{L}_s^\varepsilon)$ for $f \in \text{Dom}(\mathcal{L}_s^0)$ and $s \in [0, t]$. In order to iterate the above expansion, one needs an expansion of $\mathcal{L}^\varepsilon$, namely we assume $\mathcal{L}^\varepsilon = \mathcal{L}^0 + \varepsilon \mathcal{L}^1 + \dots + \varepsilon^{\ell+1} \mathcal{L}^{\ell+1, \varepsilon}$. More precisely, for every $f \in \mathcal{C}_c^\infty(\mathbb{R}^d) \subset \cap_{t \ge 0} \text{Dom}(\mathcal{L}_t^\varepsilon) \cap \text{Dom}(\mathcal{L}_t^0) \dots \cap \text{Dom}(\mathcal{L}_t^{\ell+1, \varepsilon})$, we assume that

$$ \mathcal{L}_t^\varepsilon f(x) = \mathcal{L}_t^0 f(x) + \varepsilon \mathcal{L}_t^1 f(x) + \dots + \varepsilon^{\ell+1} \mathcal{L}_t^{\ell+1, \varepsilon} f(x). \quad (5.17) $$

Then, under some integrability conditions that we omit, for all $(t, x) \in \mathbb{R}_+ \times \mathbb{R}^d$, the following expansion holds

$$ P_{0,t}^{\varepsilon} f(x) = P_{0,t}^{0} f(x) + \sum_{p=1}^{\ell} \varepsilon^{p} I_{p}(f)(x) + \varepsilon^{\ell+1} (R_{\ell}^{1, \varepsilon} + R_{\ell}^{2, \varepsilon}) \quad (5.18) $$

with

$$ I_p(f)(x) := \sum_{k=1}^{p} \sum_{\substack{\alpha = (\alpha_1, \dots, \alpha_k) \\ \|\alpha\| = p}} \int_{\Delta_k(t)} ds_k P_{0,s_k}^0 \mathcal{L}_{s_k}^{\alpha_k} P_{s_k, s_{k-1}}^0 \cdots P_{s_2, s_1}^0 \mathcal{L}_{s_1}^{\alpha_1} P_{s_1, t}^0 f(x) $$

and $\limsup_{\varepsilon \to 0} \{|R_\ell^{1, \varepsilon}| + |R_\ell^{2, \varepsilon}|\} < \infty$. The norm of the index $\alpha = (\alpha_1, \dots, \alpha_\ell)$, $\alpha_i \in \{1, \dots, \ell\}$, is defined as $\|\alpha\| = \sum_{j=1}^\ell \alpha_j$.

The main application discussed in this work concerns semigroups generated by diffusion processes which are defined as solutions of SDEs. In order to extend the expansion of $P_{0,t}^\varepsilon f$ to the class of bounded Borel functions, we will rely on some integration by parts formula on the Wiener space. We now explain the main idea. We will use the standard notations of Malliavin calculus. We refer to Ikeda and Watanabe [64] and Nualart [108] for a complete exposition of this topic. We consider the unique strong solution $Z$ of the following SDE

$$ Z_t = Z_s + \int_s^t b(u, Z_u) du + \int_s^t \sigma_\ell(u, Z_u) dW_u^\ell; \quad Z_0 = z \in \mathbb{R}^d, (s, t) \in [0, T]^2 \quad (5.19) $$
with coefficients $b(t, \cdot), \sigma_\ell(t, \cdot) \in C^\infty(\mathbb{R}^d, \mathbb{R}^d)$, $t \in [0, T]$, $\ell = 1, \dots, q$, with bounded derivatives of all orders greater than one. We also assume that $b(t, 0), \sigma_1(t, 0), \dots, \sigma_q(t, 0)$ are bounded on $[0, T]$. The corresponding time inhomogeneous Markov semigroup is given by $P_{s,t}f(x) := \mathbb{E}[f(Z_t(s, x))], 0 \le s \le t \le T$, for any given bounded measurable function $f$.

It is well known that the above non-homogeneous diffusion defines a $C^\infty$ non-homogeneous stochastic flow, see e.g. Kunita [80]. Therefore we may define flow derivatives of the type $\partial_{Z_{s_1}^{\alpha_1}, \dots, Z_{s_k}^{\alpha_k}}^k Z_t$ for any $\mathbf{s}_k \in [0, t]^k$ and any multi-index $\alpha \in \{1, \dots, d\}^k$, $k \in \mathbb{N}^*$. Similarly, we also define accordingly the multilinear form $\partial_{Z_s}^i Z_t$, $i \in \mathbb{N}$. Notice that we can therefore understand $\partial_{Z_{s_1}^{\alpha_1}, \dots, Z_{s_k}^{\alpha_k}}^k Z_t$ as a (random) function of the vector $(Z_{s_1}, \dots, Z_{s_k})$.

For $j \in \mathbb{N}$, $s \in [0, T]$, $\varepsilon \in [0, 1]$, we consider coefficients $B_s^j(\varepsilon, \cdot) \in C^\infty(\mathbb{R}^d, \mathbb{R}^d)$ and $A_s^j(\varepsilon, \cdot) \in C^\infty(\mathbb{R}^d, S^d(\mathbb{R}))$ such that for any multi-index $\alpha$ of elements of $\{1, \dots, d\}$ and any $j \in \mathbb{N}$,

$$ \sup_{s \in [0, T], \varepsilon \in [0, 1]} \{|\partial_{\alpha} B_s^j(\varepsilon, x)| + |\partial_{\alpha} A_s^j(\varepsilon, x)|\} \le Q_{\alpha}(|x|) \quad (5.20) $$

for some polynomial function $Q_\alpha$. We also define the corresponding non-homogeneous derivative operators

$$ \mathcal{L}_s^{j, \varepsilon} f(x) = B_s^j(\varepsilon, x) \nabla_x f(x) + \frac{1}{2} A_s^j(\varepsilon, x) \partial_x^2 f(x) \quad (5.21) $$

and note that we are using product notation: $b\nabla_x f(x) = \sum_{i=1}^d b_i \partial_{x_i} f(x)$ and $A\partial_x^2 f(x) = \sum_{i,j=1}^d A_{i,j} \partial_{x_i x_j}^2 f(x)$ for $f \in C^2(\mathbb{R}^d)$, $b \in \mathbb{R}^d$ and $A \in S^d(\mathbb{R})$. With this notation together with the existence of smooth flows the interchange of derivatives and expectations is straightforward, that is for any $f \in C_b^\infty(\mathbb{R}^d)$, one has

$$ \mathcal{L}_s^{j, \varepsilon} P_{s,t} f(x) = \mathbb{E} \left[ \mathcal{L}_{Z_s}^{j, \varepsilon} f(Z_t) \middle| Z_s = x \right], \quad \text{with} \quad \mathcal{L}_{Z_s}^{j, \varepsilon} f(Z_t) = B_s^j(\varepsilon, Z_s) \nabla_{Z_s} f(Z_t) + \frac{1}{2} A_s^j(\varepsilon, Z_s) \partial_{Z_s}^2 f(Z_t) \quad (5.22) $$

Clearly, this formula transfers derivatives operators applied to semigroups into differentiation along the flow generated by the diffusion. From here, we can also obtain the following generalised transfer formula:

$$ y \mapsto \mathcal{L}_{s_p}^{\alpha_p, \varepsilon} P_{s_p, s_{p-1}} \cdots P_{s_2, s_1} \mathcal{L}_{s_1}^{\alpha_1, \varepsilon} P_{s_1, t} f(y) = \mathbb{E} \left[ \mathcal{L}_{Z_{s_p}}^{\alpha_p, \varepsilon} \cdots \mathcal{L}_{Z_{s_1}}^{\alpha_1, \varepsilon} f(Z_t) \middle| Z_{s_p} = y \right] $$

with

$$ \mathcal{L}_{Z_{s_p}}^{\alpha_p, \varepsilon} \cdots \mathcal{L}_{Z_{s_1}}^{\alpha_1, \varepsilon} f(Z_t) = \sum_{\beta \in A_0(2p)} \partial_{\beta} f(Z_t) \Theta_{\varepsilon}^{\beta, \alpha}(\mathbf{s}_p, Z_{s_p}). $$

We use the notation $A_0(p)$ for the set of multi-indices $\alpha = (\alpha_1, \dots, \alpha_\ell) \in \{1, \dots, d\}^\ell$, with length $\ell \in \{1, \dots, p\}$. Also $\{y \mapsto \Theta_{\varepsilon}^{\beta, \alpha}(\mathbf{s}_p, y), \beta \in A_0(2p), \varepsilon \in [0, 1]\}$ is a family of smooth functions such that for any $\beta \in A_0(2p)$, $\varepsilon \in [0, 1]$ and any tuple $\gamma$ of elements of $\{1, \dots, d\}$, $\Theta_{\varepsilon}^{\beta, \alpha}: \Delta_p(t) \times \mathbb{R}^d \to \mathbb{D}^\infty$ satisfies

1. $\Theta_{\varepsilon}^{\beta, \alpha}(\mathbf{s}_p, \cdot) \in C^\infty(\mathbb{R}^d)$ and $[\partial^\gamma \Theta_{\varepsilon}^{\beta, \alpha} / \partial y^\gamma](\cdot, \cdot)$ is continuous on $\Delta_p(t) \times \mathbb{R}^d$ a.s.

2. $\forall (k, p) \in \mathbb{N} \times [1, \infty)$, there exists a polynomial $Q_{k,p}$ function such that

$$ \forall y \in \mathbb{R}^d, \quad \sup_{\mathbf{s}_p \in \Delta_p(t), \varepsilon \in [0,1]} \left\| [\partial^\gamma \Theta_{\varepsilon}^{\beta, \alpha} / \partial y^\gamma](\mathbf{s}_p, y) \right\|_{k,p} < Q_{k,p}(|y|). $$

Moreover, if one assumes that the Malliavin covariance matrix $\Gamma_{Z_T}$ of $Z_T$ is invertible and that $\det(\Gamma_{Z_T})^{-1}$ has moments of all orders. Then, making use of the Malliavin integration by parts formula, for any multi-index $\alpha$ of length $p$ and any $\mathbf{s}_p \in \Delta_p(T)$, one has

$$ P_{0,s_p} \mathcal{L}_{s_p}^{\alpha_p, \varepsilon} P_{s_p, s_{p-1}} \cdots P_{s_2, s_1} \mathcal{L}_{s_1}^{\alpha_1, \varepsilon} P_{s_1, T} f(x) = \sum_{\beta \in A_0(2p)} \mathbb{E}[f(Z_T(x)) H_\beta(Z_T(x), \Theta_{\varepsilon}^{\beta, \alpha}(\mathbf{s}_p, Z_{s_p}(x)))]. \quad (5.23) $$

We apply this approach to a $d$-dimensional diffusion processes $(X_t^\varepsilon)_{t \ge 0}$ solution of an SDE with dynamics:

$$ X_t^\varepsilon = x_0 + \int_0^t b(\varepsilon, X_s^\varepsilon) ds + \varepsilon \int_0^t \sigma_\ell(X_s^\varepsilon) dW_s^\ell, \quad (\varepsilon, t) \in [0, 1] \times [0, T], \quad (5.24) $$

where $(W_t = W_t^1, \dots, W_t^q)_{t \in [0, T]}$ is a $q$-dimensional Brownian motion, $\sigma = (\sigma_1, \dots, \sigma_q)$ is a system of $\mathbb{R}^d$-valued functions all defined on $\mathbb{R}^d$ and $b: [0, 1] \times \mathbb{R}^d \to \mathbb{R}^d$ and $\varepsilon \in [0, 1]$. Our aim is to derive an asymptotic expansion of the density $p^\varepsilon(0, T, x_0, y)$ of (5.24) starting from $x_0$ at time 0 and evaluated at terminal point $y$ at time $T$ as $\varepsilon \downarrow 0$.
Small noise as well as short time asymptotics can be derived from the heat kernel expansion of (5.24). As the first part of our main results, we obtain the expansion at any order of $p^\epsilon(0, T, x_0, y)$ by combining the parametrix approach exposed above with Malliavin's calculus techniques.

Short time expansions, that is, the case $b(\epsilon, \cdot) = \epsilon^2 b(\cdot, \cdot)$, have been investigated by Molchanov [106] and Azencott [6] for elliptic operators. An application of the Malliavin calculus to this problem was first discussed by Bismut [21] under the so-called (H2) assumption which is weaker than the ellipticity assumption. Then, Watanabe [127] obtained a similar expansion through the asymptotics of generalized Wiener functionals. Following [6] and [21], Ben Arous [16], [17], combined the Laplace method on Wiener space with the Malliavin calculus to derive short-time asymptotics for hypo-elliptic diffusion processes.

A first order expansion for the transition density of (5.24) has already been considered in Deuschel and al. [37] using the Laplace method on Wiener space, as initially developed by Ben Arous in [16] and [17] in the small time case. Here, following our approach, we derive the full expansion of $p^\epsilon(0, t, x_0, y)$ at any order under mild conditions, including unicity of the minimiser and a non-degeneracy condition (ND) (see below for precise statement) in the spirit of Bismut's hypothesis.

In the sequel we will suppose that the following assumption is in force:

(H1) $b_i \in C^\infty([0, 1] \times \mathbb{R}^d, \mathbb{R})$ and $\sigma_{i,j} \in C^\infty(\mathbb{R}^d, \mathbb{R})$, $(i,j) \in \{1, \dots, d\} \times \{1, \dots, q\}$. Moreover, for all $\epsilon \in [0, 1]$ and for all integer $n$, the functions $\partial_{\epsilon}^n b_i(\epsilon, \cdot)$, $\sigma_{i,j}(\cdot)$ have bounded derivatives of all orders greater or equal to one.

(WH) $\text{span}\{\sigma_i(x_0), i = 1, \dots, q; [\sigma_i, \sigma_j](x_0), 0 \le i, j \le q; \dots\} = \mathbb{R}^d$, where we set $\sigma_0(\cdot) = b(0, \cdot)$ for convenience.

We consider the process $(X_t^{\epsilon,h})_{t \in [0,T]}$ with dynamics

$$ \forall (\epsilon, t) \in [0, 1] \times [0, T], \quad X_t^{\epsilon,h} = x_0 + \int_0^t b(\epsilon, X_s^{\epsilon,h}) ds + \epsilon \int_0^t \sigma_\ell(X_s^{\epsilon,h}) dW_s^\ell + \int_0^t \sigma_\ell(X_s^{\epsilon,h}) \dot{h}_s^\ell ds \quad (5.25) $$

where $h$ is an element of the Cameron-Martin space $H$. Under (H1), for fixed $x$ and $h \in H$, for all integer $N$, $X_t^{\epsilon,h}$ has the asymptotic expansion

$$ X_t^{\epsilon,h} = X_t^{0,h} + \sum_{k=1}^{N} \epsilon^k \Delta_0^k X^{\cdot,h}(t) + \epsilon^{N+1} \Delta_{\epsilon}^{N+1} X^{\cdot,h}(t) \quad (5.26) $$

with $\Delta_0^k X^{\cdot,h}(t) = \frac{1}{k!} \partial_{\epsilon}^k X_t^{\epsilon,h}|_{\epsilon=0}$, residue $\Delta_{\epsilon}^{N+1} X^{\cdot,h}(t) = \frac{1}{N!} \int_0^1 d\lambda (1-\lambda)^N \partial_{\alpha}^{N+1} X_t^{\alpha,h}|_{\alpha=\lambda\epsilon}$ and the convention $\sum_{k=1}^0 = 0$.

Moreover, from the previous definition, for all integer $N$, for all integer $p \ge 1$, one has

$$ \sup_{r \in [0,1]} \mathbb{E} \left[ \sup_{0 \le t \le T} |\Delta_r^{N+1} X_t^{\cdot,h}|^p \right] < \infty. \quad (5.27) $$

In particular, the limit equation as $\epsilon$ goes to zero of the previous system is given by

$$ \begin{aligned} X_t^{0,h} &= x_0 + \int_0^t b(0, X_s^{0,h}) ds + \int_0^t \sigma_\ell(X_s^{0,h}) \dot{h}_s^\ell ds, \\ X_t^{1,h} &= \int_0^t Y_t^h (Y_s^h)^{-1} \partial_\epsilon b(0, X_s^{0,h}) ds + \int_0^t Y_r^h (Y_s^h)^{-1} \sigma_\ell(X_s^{0,h}) dW_s^\ell, \\ Y_t^h &= I_d + \int_0^t \partial_x b(0, X_s^{0,h}) Y_s^h ds + \int_0^t \partial_x \sigma_\ell(X_s^{0,h}) Y_s^h \dot{h}_s^\ell ds. \end{aligned} $$

A crucial assumption to derive the expansion of the density is that the set $K_{x_0,y} := \{h \in H : X_T^{0,h} = y\}$ is non-empty. It is satisfied under the strong Hörmander condition (see e.g. Coron [31])

(SH) $\forall x \in \mathbb{R}^d$, $\text{Lie}\{\sigma_1, \dots, \sigma_q\}(x) = \mathbb{R}^d$,

Now let $\bar{h} \in K_{x_0,y}$ which minimises the norm that is $\|\bar{h}\|_H^2 = \inf_{h \in K_{x_0,y}} \|h\|_H^2$. A standard weak-compactness argument shows that such an $\bar{h}$ always exists (see e.g. Bismut [21] or Watanabe [127]). In the spirit of Bismut [21], we introduce the following non-degeneracy assumption:

(ND) For every $u \in \mathbb{R}^q$,

$$ \mathcal{L}\{\sigma_1(x_0), \dots, \sigma_\ell(x_0), [\sigma_1, V](x_0), \dots, [\sigma_\ell, V](x_0)\} = \mathbb{R}^d, \quad \text{for } V(x_0) = \sum_{\ell=1}^q u_\ell \sigma_\ell(x_0) + b(0, x_0) $$
Assuming that $(H_1)$ and $(ND)$ are satisfied and that $h \in H \cap C^2([0, T])$ then the Malliavin covariance matrix $\Gamma_{\Delta_{\bar{\epsilon}}^1 X_T^{\bar{\epsilon},h}}$ is uniformly non degenerate in the sense that

$$ \limsup_{\bar{\epsilon} \to 0} \| \det(\Gamma_{\Delta_{\bar{\epsilon}}^1 X_T^{\bar{\epsilon},h}})^{-1} \|_p < \infty. \qquad (5.28) $$

Once we have fixed the minimiser on the set $\mathcal{K}_{x_0,y}$ of $\|\bar{h}\|_H^2$ and in the interests of simplifying notations we omit the variable denoting this element of the Cameron-Martin space in all processes. Therefore, we will have, e.g., $X^{\epsilon, \bar{h}} \equiv \bar{X}^{0, \epsilon} \equiv \bar{X}^{\epsilon}, \Delta_{\bar{\epsilon}}^i \bar{X}^{\cdot, h} \equiv \bar{X}^{i, \epsilon}$, for $\epsilon \in [0, 1]$, $X^{0, \bar{h}} \equiv \bar{X}^0$, $\Delta_0^i \bar{X}^{\cdot, \bar{h}} \equiv \bar{X}^i$ and $\bar{Y}^{\bar{h}} \equiv \bar{Y}$.

In particular, the “deterministic” Malliavin covariance matrix $\Gamma_{\bar{X}_T^1}$ of $\bar{X}_T^1$ defined by

$$ (\Gamma_{\bar{X}_T^1})_{i,j} = \sum_{\ell=1}^{q} \int_{0}^{T} (\bar{Y}_T(\bar{Y}_s)^{-1} \sigma_{\ell}(\bar{X}_s^0))_i (\bar{Y}_T(\bar{Y}_s)^{-1} \sigma_{\ell}(\bar{X}_s^0))_j ds, \quad (i,j) \in \{1, \dots, d\}^2 \qquad (5.29) $$

is invertible under (ND).

Moreover, in order to characterise some properties of the minimiser, one uses the Lagrange multiplier method which ensures that

$$ \exists \lambda \in \mathbb{R}^d \text{ s.t. for a.a. } t \in [0, T], \quad \dot{h}_t^\ell = \langle \lambda, \bar{Y}_T(\bar{Y}_t)^{-1} \sigma_\ell(\bar{X}_t^0) \rangle, \quad \ell = 1, \dots, q. \qquad (5.30) $$

Using the dynamics of $(X_t^{1,h})_{t \in [0,T]}$ with (5.30), we deduce the following crucial identity

$$ \int_0^T \langle \dot{h}_s, dW_s \rangle = \langle \lambda, \int_0^T \bar{Y}_T(\bar{Y}_s)^{-1} \sigma_\ell(\bar{X}_s^0) dW_s^\ell \rangle = \langle \lambda, \bar{X}_T^1 \rangle - \langle \lambda, \int_0^T \bar{Y}_T(\bar{Y}_s)^{-1} \partial_\epsilon b(0, \bar{X}_s^0) ds \rangle. \qquad (5.31) $$

Finally we introduce the following condition

$$ (H_2) \quad \mathbb{E} [\exp\langle\lambda, \bar{X}_T^2\rangle | \bar{X}_T^1 = 0] < \infty. $$

We now give a guideline of the arguments to prove the main result (Theorem 5.3.1) of this section. Let $\delta > 0$ to be fixed later on. Under (WH), the density of $X_T^\epsilon$ may be decomposed as follows

$$ p^{\epsilon}(0, T, x_0, y) = \mathbb{E}[\delta_y(X_T^{\epsilon})] = \mathbb{E}[\psi_{\delta}(Z_T^{1, \epsilon})\delta_y(X_T^{\epsilon})] + \mathbb{E}[(1 - \psi_{\delta}(Z_T^{1, \epsilon}))\delta_y(X_T^{\epsilon})] \qquad (5.32) $$

where $\psi \in C^\infty(\mathbb{R})$ is a bump function such that $0 \le \psi(x) \le 1$, $\psi(x) = 1$ on $|x| \le 1/2$ and $\psi(x) = 0$ for $|x| > 1$, for $\delta > 0$, $\psi_\delta(x) := \psi(x/\delta)$ and we introduced the following process $Z_t^{1, \epsilon} = \int_0^t |X_s^\epsilon - \bar{X}_s^0|^2 ds, t \in [0, T]$. The above decomposition is necessary in order to deal with the remainder terms in the expansion in powers of $\epsilon$ of the first term appearing in the right-hand side of (5.32).

Using Malliavin's integration by parts and a large deviation arguments, one can prove that $0 \le \mathbb{E}[(1 - \psi_\delta(Z_T^{1, \epsilon}))\delta_y(X_T^\epsilon)] \le C \exp(-(\|\bar{h}\|_H^2 + c)/(2\epsilon^2))$, for some positive constant $c$, $C$, so that, as it will be clear later on, this term only appears in the remainder of the asymptotic expansion (at any order).

We now focus on the first term appearing in the right hand side of (5.32). By the Girsanov Theorem with drift equals to $\bar{h}/\epsilon$, $\bar{h} \in \mathcal{K}_{x_0,y}$, (5.26) with $N = 0, 1$ and (5.31), one successively deduces

$$ \begin{aligned} \mathbb{E}[\psi_{\delta}(Z_T^{1, \epsilon})\delta_y(X_T^{\epsilon})] &= \exp\left(-\frac{\|\bar{h}\|_H^2}{2\epsilon^2}\right) \mathbb{E}\left[\exp\left(-\frac{1}{\epsilon}\int_0^T \langle \bar{h}_s, dW_s\rangle\right) \psi_{\delta}(\bar{Z}_T^{1, \epsilon})\delta_y(\bar{X}_T^{\epsilon})\right] \\ &= \epsilon^{-d} \exp\left(-\frac{\|\bar{h}\|_H^2}{2\epsilon^2}\right) \mathbb{E}\left[\exp\left(-\frac{1}{\epsilon}\int_0^T \langle \bar{h}_s, dW_s\rangle\right) \psi_{\delta}(\bar{Z}_T^{1, \epsilon})\delta_0(\bar{X}_T^{1, \epsilon})\right] \\ &= \epsilon^{-d} \exp\left(-\frac{\|\bar{h}\|_H^2}{2\epsilon^2}\right) \mathbb{E}\left[\exp\left(-\frac{1}{\epsilon}\int_0^T \langle \lambda, \bar{Y}_T(\bar{Y}_s)^{-1} \sigma_{\ell}(\bar{X}_s^0) dW_s^{\ell}\rangle\right) \psi_{\delta}(\bar{Z}_T^{1, \epsilon})\delta_0(\bar{X}_T^{1, \epsilon} + \epsilon \bar{X}_T^{2, \epsilon})\right] \\ &= \epsilon^{-d} \exp\left(-\frac{\|\bar{h}\|_H^2}{2\epsilon^2}\right) \exp\left(\frac{1}{\epsilon}\int_0^T \bar{Y}_T(\bar{Y}_s)^{-1} \partial_{\epsilon} b(0, \bar{X}_s^0) ds\right) \mathbb{E}\left[e^{-\frac{\langle \lambda, \bar{X}_T^1 \rangle}{\epsilon}} \psi_{\delta}(\bar{Z}_T^{1, \epsilon})\delta_0(\bar{X}_T^{1, \epsilon} + \epsilon \bar{X}_T^{2, \epsilon})\right] \\ &= \epsilon^{-d} \exp\left(-\frac{\|\bar{h}\|_H^2}{2\epsilon^2}\right) \exp\left(\frac{1}{\epsilon}\int_0^T \bar{Y}_T(\bar{Y}_s)^{-1} \partial_{\epsilon} b(0, \bar{X}_s^0) ds\right) \mathbb{E}\left[e^{\langle \lambda, \bar{X}_T^{2, \epsilon} \rangle} \psi_{\delta}(\bar{Z}_T^{1, \epsilon})\delta_0(\bar{X}_T^{1, \epsilon})\right]. \end{aligned} \qquad (5.33) $$

In the above formula, we extracted the main term of the asymptotic expansion. To conclude we prove that $\mathbb{E}[e^{\langle \lambda, \bar{X}_T^{2, \epsilon} \rangle} \psi_{\delta}(\bar{Z}_T^{1, \epsilon}) \delta_0(\bar{X}_T^{1, \epsilon})]$ admits an expansion in powers of $\epsilon$ using the framework developed above. Roughly
speaking we introduce the Markov semigroup $(P_t^\varepsilon)_{t \in [0,T]}$ associated to the trivariate process $(\bar{X}_t^{1,\varepsilon}, \bar{Z}_t^{1,\varepsilon}, \bar{X}_t^{2,\varepsilon})_{t \in [0,T]}$. We prove that its infinitesimal generator admits an expansion in powers of $\varepsilon$. We use the integration by parts formula (5.23) w.r.t the non-degenerate process $\bar{X}^1$ in order to get rid of the derivatives on the test function. In order to deal with the remainder term in the expansion, one has to use a localised version of the usual Malliavin's integration by parts and the asymptotic non-degeneracy (5.28). Finally we obtain the following expansion.

**Theorem 5.3.1.** Let $x_0, y \in \mathbb{R}^d$. Assume that $(\mathbf{H}_1)$, $(\mathbf{WH})$, $(\mathbf{H}_2)$ and $(\mathbf{ND})$ hold. Suppose that $\mathcal{K}_{x_0,y} \neq \emptyset$ (which is guaranteed under (SH)) and that the minimizer $\bar{h} = \arg \min_{h \in \mathcal{K}_{x_0,y}} \|h\|_{H}^2$ is unique. Then, for all $\ell \in \mathbb{N}$, one has

$$p^{\varepsilon}(0, T, x, y) = \varepsilon^{-d} \exp \left( -\frac{\|\bar{h}\|_{H}^{2}}{2\varepsilon^{2}} \right) \exp \left( \frac{1}{\varepsilon} \langle \lambda, \int_{0}^{T} \bar{Y}_{T}(\bar{Y}_{s})^{-1} \partial_{\varepsilon} b(0, \bar{X}_{s}^{0}) ds \rangle \right) \left( c_{0}(x, y) + \varepsilon c_{1}(x, y) + \dots + \varepsilon^{\ell} c_{\ell}(x, y) + \mathcal{O}(\varepsilon^{\ell+1}) \right).$$

We can compute the exact expressions of the coefficients. One notably has

$$c_0(x, y) = \mathbb{E} \left[ \delta_0(\bar{X}_T^1) e^{\langle \lambda, \bar{X}_T^2 \rangle} \right], \\ c_1(x, y) = \mathbb{E} \left[ \delta_0(\bar{X}_T^1) e^{\langle \lambda, \bar{X}_T^2 \rangle} \sum_{\beta_1=1}^{2d} \left\{ \lambda^{\beta_1} \tilde{H}_{\tilde{\beta}_1}(\bar{X}_T^1, \Theta^{\beta_1}(s_1, \bar{X}_{s_1}), \langle \lambda, \bar{X}_T^2 \rangle) + \sum_{\beta_2=1}^{2d} \lambda^{\beta_2} \tilde{H}_{\tilde{\beta}_2}(\bar{X}_T^1, \Theta^{\beta}(s_1, \bar{X}_{s_1}), \langle \lambda, \bar{X}_T^2 \rangle) \right\} \right].$$

where, $(s, x) \mapsto \Theta^{\beta_1}(s, x)$, $\Theta^{\beta}(s, x)$ are two explicit smooth random functions and we interpret $\tilde{\beta}_1 = \emptyset$ if $\beta_1 \in \{d+1, \dots, 2d\}$ and that $H_{\emptyset} = 1$.

As a second application we consider a one dimensional Markov semigroup with singular coefficients for which we derive an asymptotic expansion with respect to a small parameter $\varepsilon$. We consider the following skew diffusion with small skew

$$X_t^\varepsilon = x + \int_0^t b(X_s^\varepsilon) ds + \int_0^t \sigma(X_s^\varepsilon) dW_s + \varepsilon L_t^0(X^\varepsilon), \quad \varepsilon \in [0, 1]. \qquad (5.34)$$

where $W = (W_t)_{t \ge 0}$ is a one-dimensional standard Brownian motion and $L^0(X^\varepsilon) = (L_t^0(X^\varepsilon))_{t \ge 0}$ is the symmetric local time of $X^\varepsilon$ at the origin. We work under assumptions (HR) and (HE) presented in Section 5.2 which guarantees the existence of a unique weak solution to (5.34). The first order approximation is the semigroup $(P_t^0)_{t \ge 0}$ associated to $\mathcal{L}_s^0 = \mathcal{L}^0 = \frac{1}{2}a\partial_x^2 + b\partial_x$. For $f \in C_c^\infty(\mathbb{R})$, we prove the following first step expansion

$$P_t^\varepsilon f(x) = P_t^0 f(x) + \varepsilon \int_0^t \int_{\mathbb{R}} \frac{p^\varepsilon(0, s, x, 0+) + p^\varepsilon(0, s, x, 0-)}{2} \partial_x p^0(s, t, 0, y) f(y) dy ds$$

or equivalently,

$$p^{\varepsilon}(0, t, x, y) = p^{0}(0, t, x, y) + a(0)\varepsilon \int_{0}^{t} ds_{1} \frac{p^{\varepsilon}(0, s_{1}, x, 0+) + p^{\varepsilon}(0, s_{1}, x, 0-)}{2} \partial_{x}p^{0}(s_{1}, t, 0, y).$$

We understand that the expansion (5.17) is satisfied in a generalised sense with $\mathcal{L}^{1,\varepsilon} = \mathcal{L}^1 = a(0)\delta_0\partial_x$. Then, using some sharp estimates on $\lim_{\eta \downarrow 0} \frac{1}{2}(\partial_x p^0(s_1, t, 0, \eta) + \partial_x p^0(s_1, t, 0, -\eta))$, we prove that we can iterate the above first step expansion.

**Theorem 5.3.2.** Under (HR) and (HE), for all $\varepsilon \in [0, 1)$, the transition density of (5.34) can be expressed as

$$\forall (t, x, y) \in \mathbb{R}_+^* \times \mathbb{R} \times \mathbb{R}^*, \quad p^\varepsilon(0, t, x, y) = p^0(0, t, x, y) + \sum_{p \ge 1} (a(0)\varepsilon)^p \int_{\Delta_p(t)} ds_p \prod_{\ell=1}^{p-1} \partial_x p^0(s_{\ell+1}, s_\ell, 0, 0) p^0(0, s_p, x, 0) \partial_x p^0(s_1, t, 0, y)$$

with $s_0 = t$. Alternatively, for all positive integer $n$, one has

$$\forall (t, x, y) \in \mathbb{R}_+^* \times \mathbb{R} \times \mathbb{R}^*, \quad p^\varepsilon(0, t, x, y) = p^0(0, t, x, y) + \sum_{p=1}^{n} (a(0)\varepsilon)^p \int_{\Delta_p(t)} ds_p \prod_{\ell=1}^{p-1} \partial_x p^0(s_{\ell+1}, s_\ell, 0, 0) p^0(0, s_p, x, 0) \partial_x p^0(s_1, t, 0, y) \\ + \varepsilon^{n+1} t^{(n+1)\frac{\eta}{2}-\frac{1}{2}} R_n^{1,\varepsilon}$$

with $\limsup_{\varepsilon \to 0} |R_n^{1,\varepsilon}| < +\infty$.
# Chapter 6

## One-dimensional elliptic diffusions and some of their path-functionals

In this last chapter we present the results of [13] in collaboration with A. Kohatsu-Higa from Ritsumeikan university and L. Li from the university of New South Wales and [14] in collaboration with L. Li. In both works we are interested in Markov perturbation problems related to some path-functionals of a one-dimensional elliptic diffusion.

In Section 6.1, we expose the results obtained in [13] related to the joint law of the first hitting time of a threshold by a one-dimensional uniformly elliptic diffusion process and to the associated process stopped at the threshold. Our methodology relies on the parametrix method that we apply to the associated Markov semigroup. In the spirit of Bally and Kohatsu-Higa [7], see Section 1.2.5 for a brief account, we obtain a probabilistic representation that may be useful for the construction of an unbiased Monte Carlo path simulation method. Some integrations by parts formula are also established under mild assumptions on the coefficients.

In Section 6.2, we present the results obtained in [14] related to the weak existence and uniqueness as well as density estimates for SDEs with coefficients depending on some path-functionals of the process. We consider two examples: a diffusion process with coefficients depending on its running symmetric local time and a diffusion process with coefficients depending on its running maximum. In each example, we also prove the existence of the associated transition density and establish some Gaussian upper-estimates.

### 6.1 Hitting time and killed process associated to a one-dimensional diffusion

In this section, we expose the results obtained in [13] in collaboration with A. Kohatsu-Higa and L. Li related to the first hitting time of a one-dimensional diffusion process with dynamics

$$X_t^x = x + \int_0^t b(X_s^x)ds + \int_0^t \sigma(X_s^x)dW_s, \quad x \in \mathbb{R} \qquad (6.1)$$

where $(W_t)_{t \ge 0}$ stands for a one-dimensional Brownian motion on a given filtered probability space $(\Omega, \mathcal{F}, (\mathcal{F}_t)_{t \ge 0}, \mathbb{P})$. We denote by $\tau^x := \inf\{t \ge 0, X_t^x \ge L\}$ the first hitting time of the level $L$ (or equivalently the exit time of the open set $(-\infty, L)$) by the one-dimensional process $X$ and by $(X_{\tau_t^x})_{t \ge 0}$ the associated killed process. Here we introduced the notation $\tau_t^x := \tau^x \wedge t$. In order to study the joint law of $(\tau_t^x, X_{\tau_t^x})_{t \ge 0}$, we introduce the collection of linear maps $(P_t)_{t \ge 0}$, acting on $B_b(\mathbb{R})$, as follows

$$\forall (u, x) \in \mathbb{R}_+ \times \mathbb{R}, \quad P_t h(u, x) = \mathbb{E} \left[ h \left( u + \tau_t^x, X_{\tau_t^x}^x \right) \right]. \qquad (6.2)$$

In a multi-dimensional setting and for a domain $D$ such that $\partial D$ is smooth and noncharacteristic, Cattiaux [28] developed a Malliavin's calculus approach to prove that the semigroup associated to a process killed when it hits the boundary $\partial D$ admits an infinitely differentiable kernel under a restricted Hörmander condition on the vector fields. Gaussian bounds on this kernel are also established in small time. We also refer the reader to Ladyzenskaja and al. [85], Friedman [44] and Garroni and Menaldi [46] for constructions of Green functions related to a class of Cauchy-Dirichlet value problems in a uniformly elliptic setting using a partial differential equation framework. In particular, in Chapter IV in the latter reference, the authors introduced the notion of Poisson and Green functions related to a similar Cauchy-Dirichlet problem and in Chapter VI existence and Gaussian upper-bounds are established when the coefficients $b$ and $\sigma$ are Hölder continuous and $a = \sigma^2$ is uniformly elliptic. As a slight extension of the latter result, in
[13], we prove the existence of a transition density for $(P_t)_{t \ge 0}$, thus proving existence of a density for $(\tau_t^x, X_{\tau_t^x})$, $t > 0$, under the assumption that $b$ is bounded measurable and $\sigma$ is Hölder continuous and $a$ is uniformly elliptic. Some Gaussian upper-bounds for both densities are also established. In this section, we rather focus on the probabilistic representation for $(P_t)_{t \ge 0}$.

The fundamental idea in order to obtain a probabilistic representation for the joint law $(\tau_t^x, X_{\tau_t^x})$ is to see the operator $P_t$ as a perturbation of the operator $\bar{P}_t$ defined by $\bar{P}_t h(u, x) = \mathbb{E}[h(u + \bar{\tau}_t^x, \bar{X}_{\bar{\tau}_t^x})]$ where $(\bar{X}_t)_{t \ge 0}$ is obtained from (6.7) by removing the drift and freezing the diffusion coefficient to the initial point $x$, namely $\bar{X}_t^x = x + \sigma(x)W_t$ and $\bar{\tau}^x$ is its associated first hitting time. Note that from the reflection principle of Brownian motion, the kernel associated to $\bar{P}_t$ is given by

$$
\begin{align*}
\bar{P}_t h(u, x) &= \mathbf{I}_{\{x \ge L\}} h(u, x) + \mathbb{E}[h(u + t, \bar{X}_t^x) \mathbf{I}_{\{\bar{\tau}^x > t\}}] + \mathbb{E}[h(u + \bar{\tau}^x, L) \mathbf{I}_{\{\bar{\tau}^x \le t\}}] \\
&= \mathbf{I}_{\{x \ge L\}} h(u, x) + S_t^x h(u, x) + K_t^x h(u, x) \\
&= \mathbf{I}_{\{x \ge L\}} h(u, x) + \mathbf{I}_{\{x < L\}} \int_u^{u+t} \int_{-\infty}^L h(s, z) [f_{\bar{\tau}}^x(x, s-u) \delta_L(dz) ds + \bar{q}_t^x(x, z) dz \delta_{u+t}(ds)].
\end{align*}
$$

with

$$
\begin{align*}
K_t^x h(u, x) &= \mathbf{I}_{\{x<L\}} \int_0^t h(u+s,L) f_{\bar{\tau}}^x(x,s) ds, \quad f_{\bar{\tau}}^y(x,s) := \partial_s \mathbb{P}(\bar{\tau}^x \le s) = \frac{L-x}{\sqrt{2\pi a(y)s^{3/2}}} \exp\left(-\frac{(L-x)^2}{2a(y)s}\right) \mathbf{I}_{\{x \le L\}} \\
S_t^x h(u, x) &= \mathbf{I}_{\{x<L\}} \int_{-\infty}^L h(u+t,z) \bar{q}_t^x(x,z) dz, \quad \bar{q}_t^y(x,z) := (g(a(y)t,z-x) - g(a(y)t,z+x-2L)) \mathbf{I}_{\{x \le L\}} \mathbf{I}_{\{z \le L\}},
\end{align*}
$$

and $g(a, x) = (2\pi a)^{-\frac{1}{2}} \exp(-\frac{1}{2}a^{-1}x^2)$. A direct application of Itô's lemma allows to characterise the infinitesimal generators $\mathcal{L}$ of $P_t$ and $\bar{\mathcal{L}}^x$ of $\bar{P}_t$, namely for $h \in C_b^{1,2}(\mathbb{R}_+ \times (-\infty, L])$ and $(u, x) \in \mathbb{R}_+ \times \mathbb{R}$

$$
\begin{align*}
\frac{P_t h(u, x) - h(u, x)}{t} &\rightarrow \mathcal{L}h(u, x) := \mathbf{I}_{\{x<L\}} \left( b(x)\partial_2 h(u, x) + \frac{1}{2}a(x)\partial_2^2 h(u, x) + \partial_1 h(u, x) \right), \\
\frac{\bar{P}_t h(u, x) - h(u, x)}{t} &\rightarrow \bar{\mathcal{L}}^x h(u, x) = \mathbf{I}_{\{x<L\}} \left( \frac{1}{2}a(x)\partial_2^2 h(u, x) + \partial_1 h(u, x) \right)
\end{align*}
$$

as $t \downarrow 0$. We assume that $a \in C_b^2(\mathbb{R})$, $b \in C_b^1(\mathbb{R})$. Roughly speaking, the first step of the forward parametrix method writes as follows

$$
P_t h(u, x) = \bar{P} h(u, x) + \int_0^t \partial_s (\bar{P}_{t-s} P_s h(u, x)) ds = \bar{P} h(u, x) + \int_0^t \bar{P}_{t-s} (\mathcal{L} - \bar{\mathcal{L}}^x) P_s h(u, x)) ds
$$

where $(\mathcal{L} - \bar{\mathcal{L}}^x)h(v, y) = \frac{1}{2}(a(y) - a(x))\partial_2^2 h(v, y) + b(y)\partial_2 h(v, y)$. Now, by an integration by parts formula, one can remove the derivatives on $P_s h$ and put them on the kernel of $\bar{P}_{t-s}$. Taking care of border terms, one gets

$$
P_T h(u, x) - \bar{P}_T h(u, x) = \int_0^T (\bar{K}_{T-s} P_s h(u, x) + \bar{S}_{T-s} P_s h(u, x)) ds. \quad (6.3)
$$

with

$$
\begin{align*}
\bar{K}_t h(u, x) &= \bar{K}_t(x, L)h(u + t, L), \quad \bar{K}_t(x, L) = \mathbf{I}_{\{x<L\}} \frac{(a(L) - a(x))}{2} \frac{2(L-x)}{a(x)t} g(a(x)t, L-x) = \mathbf{I}_{\{x<L\}} \frac{(a(L) - a(x))}{a(x)} f_{\bar{\tau}}^x(x, t), \\
\bar{S}_t h(u, x) &= \int_{-\infty}^L h(u + t, z) \bar{S}_t(x, z) dz, \quad \bar{S}_t(x, z) = \mathbf{I}_{\{x<L\}} \left\{ \frac{1}{2} \partial_z^2 \left[ (a(z) - a(x)) \bar{q}_t^x(x, z) \right] - \partial_z [b(z) \bar{q}_t^x(x, z)] \right\}.
\end{align*}
$$

Moreover, one has the following controls of the two kernels: $|\bar{K}_t(x, L)| \le C \frac{|L-x|^2}{t} g(a(x)t, L-x) \mathbf{I}_{\{x<L\}} \le C \mathbf{I}_{\{x<L\}} g(ct, L-x)$ and $|\bar{S}_t(x, z)| \le C_T \mathbf{I}_{\{x<L\}} \frac{1}{t^2} g(ct, z-x)$. In order to iterate the first step expansion (6.3), we remark that $\bar{K}_s h(u, L) = \bar{S}_s h(u, L) = 0$ so that $\bar{K}_{T-s_1} \bar{K}_{s_1-s_2} h(u, x) = \bar{K}_{T-s_1} \bar{S}_{s_1-s_2} h(u, x) = 0$. We also define the following family of operators for $(s_0, u, x) \in [0, T] \times \mathbb{R}_+ \times (-\infty, L]$:

$$
I_{s_0}^n h(u, x) = \begin{cases} \int_{\Delta_n(s_0)} ds_n \left\{ \left( \prod_{i=0}^{n-1} \bar{S}_{s_i - s_{i+1}} \right) \bar{P}_{s_n} h(u, x) + \left( \prod_{i=0}^{n-2} \bar{S}_{s_i - s_{i+1}} \right) \bar{K}_{s_{n-1} - s_n} h(u, x) \right\} & \text{if } n \ge 1, \\ \bar{P}_{s_0} h(u, x) & \text{if } n = 0, \end{cases} \quad (6.4)
$$

where $ds_p = ds_1 \cdots ds_p$ and for a fixed time $t \ge 0$, we denote by $\Delta_p(t) = \{s_p \in [0, t]^p : s_{p+1} := 0 \le s_p \le s_{p-1} \le \cdots \le s_1 \le t =: s_0\}$.
In [13], we prove that the series $\sum_{n \ge 0} I_T^n h(u, x)$ converges absolutely and uniformly and iterating (6.3) one obtains

$$ \forall h \in C_b^{1,2}(\mathbb{R}_+ \times (-\infty, L]), \quad P_t h(u, x) = \sum_{n \ge 0} I_t^n h(u, x), \quad (t, u, x) \in [0, T] \times \mathbb{R}_+ \times \mathbb{R}. \quad (6.5) $$

Actually we will work with the following expression of $I_t^n h(u, x)$. By Fubini's theorem and the change of variable $t_i = t - s_i$, $i = 0, \dots, n$, one has

$$ I_t^n h(u, x) = \int_{\Delta_n(t)} d\mathbf{s}_n \left( \prod_{i=0}^{n-1} \bar{S}_{s_i - s_{i+1}} \right) \bar{P}_{s_n} h(u, x) = \int_{\Delta_n^*(t)} dt_n \left( \prod_{i=0}^{n-1} \bar{S}_{t_{i+1} - t_i} \right) \bar{P}_{t-t_n} h(u, x) $$

where we used the notation $\Delta_p^*(t) = \{s_p \in [0, t]^p : s_0 := 0 \le s_1 \le s_2 \le \dots \le s_p \le t =: s_{p+1}\}$. We begin by providing a probabilistic interpretation of each kernel appearing in (6.4). First one has,

$$ \begin{aligned} \bar{S}_t(x, z) &= \mathbf{I}_{\{x<L\}} \left\{ \partial_z^2 \left[ \frac{1}{2}(a(z) - a(x))\bar{q}_t^x(x, z) \right] - \partial_z [b(z)\bar{q}_t^x(x, z)] \right\} \\ &= \mathbf{I}_{\{x<L\}} \left\{ \left( \frac{1}{2}a''(z) - b'(z))\bar{q}_t^x(x, z) + (a'(z) - b(z))\partial_z \bar{q}_t^x(x, z) + \frac{1}{2}(a(z) - a(x))\partial_z^2 \bar{q}_t^x(x, z) \right) \right\} \\ &= \bar{\theta}_t(x, z)g(a(x)t, z - x). \end{aligned} \quad (6.6) $$

$$ \theta_t(x, z) := \mathbf{I}_{\{x<L\}} \Lambda_t(x, z) \left\{ \left( \frac{1}{2} a''(z) - b'(z) \right) + (a'(z) - b(z)) \mu_t^1(x, z) + \frac{1}{2} (a(z) - a(x)) \mu_t^2(x, z) \right\}, $$

$$ \mu_t^1(x, z) := H_1(a(x)t, z - x) - \frac{1}{a(x)t} \frac{2(L - x)}{\exp\left(-\frac{2(z-L)(L-x)}{a(x)t}\right) - 1}, $$

$$ \mu_t^2(x, z) := H_2(a(x)t, z - x) + \frac{1}{a^2(x)t^2} \frac{4(z - L)(L - x)}{\exp\left(-\frac{2(z-L)(L-x)}{a(x)t}\right) - 1}, $$

$$ \Lambda_t(x, z) := \mathbb{P} \left( \max_{0 \le v \le t} W_v \le \frac{L-x}{\sigma(x)} \mid W_t = \frac{z-x}{\sigma(x)} \right) = \left\{ 1 - \exp \left( -2 \frac{(L-x)(L-x-(z-x))}{ta(x)} \right) \right\} \mathbf{I}_{\{x \le L\}} \mathbf{I}_{\{z \le L\}}. $$

From the expression above, for $h \in \mathcal{B}_b(\mathbb{R}_+ \times (-\infty, L])$, one gets $\bar{S}_T h(u, x) = \mathbb{E}[h(u+T, \bar{X}_T^x)\bar{\theta}_T(x, \bar{X}_T^x)]$. For a given time partition $\pi: 0 = t_0 < t_1 < \dots < t_N < t_{N+1} = T$, we introduce the Euler scheme $\bar{X}^\pi = (\bar{X}_{t_i}^\pi)_{0 \le i \le N+1}$ with the following dynamics

$$ \bar{X}_{t_{i+1}}^\pi = \bar{X}_{t_i}^\pi + \sigma(\bar{X}_{t_i}^\pi)(W_{t_{i+1}} - W_{t_i}), \quad \bar{X}_{t_0}^\pi = x $$

which in turn allows us to write the following (partial) probabilistic representation

$$ \begin{aligned} \left( \prod_{i=0}^{n-1} \bar{S}_{t_{i+1}-t_i} \right) \bar{P}_{T-t_n} h(u, x) &= \mathbb{E}[h(u+T, \bar{X}_T^\pi) \mathbf{I}_{\{t_n+\bar{\tau}^{\bar{X}_{t_n}^\pi} > T\}} \bar{\theta}_{t_n-t_{n-1}} (\bar{X}_{t_{n-1}}^\pi, \bar{X}_{t_n}^\pi) \cdots \bar{\theta}_{t_1} (x, \bar{X}_{t_1}^\pi)] \\ &\quad + \mathbb{E}[h(u+t_n, \bar{\tau}^{\bar{X}_{t_n}^\pi}, L) \mathbf{I}_{\{t_n+\bar{\tau}^{\bar{X}_{t_n}^\pi} \le T\}} \bar{\theta}_{t_n-t_{n-1}} (\bar{X}_{t_{n-1}}^\pi, \bar{X}_{t_n}^\pi) \cdots \bar{\theta}_{t_1} (x, \bar{X}_{t_1}^\pi)] \\ &= \mathbb{E}[h(u+(t_n+\bar{\tau}^{\bar{X}_{t_n}^\pi}) \wedge T, \bar{X}^\pi_{(t_n+\bar{\tau}^{\bar{X}_{t_n}^\pi}) \wedge T}) \bar{\theta}_{t_n-t_{n-1}} (\bar{X}_{t_{n-1}}^\pi, \bar{X}_{t_n}^\pi) \cdots \bar{\theta}_{t_1} (x, \bar{X}_{t_1}^\pi)], \end{aligned} $$

where $\bar{\tau}^x$ is a Lévy distributed random variable with scale parameter $c = (L-x)^2/a(x)$ and location parameter $\mu = 0$ which is independent of the underlying Brownian motion $W$. Similarly, for the second term appearing in (6.4), after noting that

$$ \int_0^{s_{n-1}} \bar{K}_{s_{n-1}-s_n} h(u, x) ds_n = \mathbb{E}[h(u + \bar{\tau}^x, L) \frac{(a(L) - a(x))}{a(x)} \mathbf{I}_{\{\bar{\tau}^x \le s_{n-1}\}}] =: \tilde{K}_{s_{n-1}} h(u, x) $$

we can write

$$ \int_{\Delta_n(T)} d\mathbf{s}_n \left( \prod_{i=0}^{n-2} \bar{S}_{s_i - s_{i+1}} \right) \bar{K}_{s_{n-1}-s_n} h(u, x) = \int_{\Delta_{n-1}^*(T)} dt_n \left( \prod_{i=0}^{n-2} \bar{S}_{t_{i+1}-t_i} \right) \check{K}_{T-t_{n-1}} h(u, x). $$

Similarly to the previous term, one obtains the following partial probabilistic representation

$$ \left( \prod_{i=0}^{n-2} \bar{S}_{t_{i+1}-t_i} \right) \check{K}_{T-t_{n-1}} h(u, x) = \mathbb{E} \left[ h(u + t_{n-1} + \bar{\tau}^{\bar{X}_{t_{n-1}}^{\pi}}, L) \mathbf{I}_{\{\bar{\tau}^{\bar{X}_{t_{n-1}}^{\pi}} \le T-t_{n-1}\}} \frac{(a(L) - a(\bar{X}_{t_{n-1}}^{\pi}))}{a(\bar{X}_{t_{n-1}}^{\pi})} \right. \\ \left. \times \bar{\theta}_{t_{n-1}-t_{n-2}}(\bar{X}_{t_{n-2}}^{\pi}, \bar{X}_{t_{n-1}}^{\pi}) \cdots \bar{\theta}_{t_1}(x, \bar{X}_{t_1}^{\pi}) \right]. $$
Finally, in the spirit of Bally and Kohatsu-Higa [7], we let $(N(t))_{t \ge 0}$ be a simple Poisson process with intensity parameter $\lambda > 0$ and define $N \equiv N(t)$ for $t \in [0, T]$. Let $\zeta_1 < \zeta_2 < \cdots < \zeta_N$ be the jump times of the Poisson process and set $\zeta_0 = 0$, $\zeta_{N+1} = T$. We still denote by $\pi$ the random time partition $\pi: \zeta_0 = 0 < \zeta_1 < \cdots < \zeta_{N+1} = T$ and denote by $\bar{X}^\pi = (\bar{X}_{\zeta_i}^\pi)_{0 \le i \le N+1}$ its associated Euler scheme. In [13], we obtain the following theorem.

**Theorem 6.1.1.** Let $T > 0$. Assume that $a \in C_b^2(\mathbb{R})$, $b \in C_b^1(\mathbb{R})$ and $a$ is uniformly elliptic. Define the two sequences $(\Gamma_N(x))_{N \ge 0}$ and $(\bar{\Gamma}_N(x))_{N \ge 0}$ as follows

$$ \Gamma_N(x) = \begin{cases} \prod_{j=0}^{N-1} \lambda^{-1} \bar{\theta}_{\zeta_{j+1}-\zeta_j} (\bar{X}_{\zeta_j}^{\pi}, \bar{X}_{\zeta_{j+1}}^{\pi}) & \text{if } N \ge 1, \\ 1 & \text{if } N = 0, \end{cases} $$

and

$$ \bar{\Gamma}_N(x) = \begin{cases} \frac{(a(L)-a(\bar{X}_{\zeta_{N-1}}^{\pi}))}{a(\bar{X}_{\zeta_{N-1}}^{\pi})} \prod_{j=0}^{N-2} \lambda^{-1} \bar{\theta}_{\zeta_{j+1}-\zeta_j} (\bar{X}_{\zeta_j}^{\pi}, \bar{X}_{\zeta_{j+1}}^{\pi}) & \text{if } N \ge 1, \\ 0 & \text{if } N = 0. \end{cases} $$

Then, for all $h \in B_b(\mathbb{R}_+ \times \mathbb{R})$, for all $(u, x) \in \mathbb{R}_+ \times \mathbb{R}$, the following probabilistic representation holds

$$ \mathbb{E}[h(\tau_t^x, X_{\tau_t^x}^x)] = e^{\lambda t} \mathbb{E} \left[ h((\zeta_{N(t)} + \bar{\tau}^{\bar{X}_{\zeta_{N(t)}}^{\pi}}) \wedge t, \bar{X}_{(\zeta_{N(t)}+\bar{\tau}^{\bar{X}_{\zeta_{N(t)}}^{\pi}}) \wedge t}^{\pi}) \Gamma_{N(t)}(x) \right] \\ + e^{\lambda t} \mathbb{E} \left[ h(\zeta_{N(t)-1} + \bar{\tau}^{\bar{X}_{\zeta_{N(t)-1}}^{\pi}}, L) \mathbf{I}_{\{\bar{\tau}^{\bar{X}_{\zeta_{N(t)-1}}^{\pi}} \le t - \zeta_{N(t)-1}\}} \bar{\Gamma}_{N(t)}(x) \right]. $$

Similarly, for all $(s, x, z) \in (0, t] \times (-\infty, L]^2$, the following probabilistic representation for the probability law of $(\tau_t^x, X_{\tau_t^x}^x)$ is satisfied

$$ p_t(0, x, ds, dz) = p_t^D(x, z)dz\delta_t(ds) + p^K(x, s)\delta_L(dz)ds $$

with, for all,

$$ p_t^D(x, z) = e^{\lambda t} \mathbb{E} \left[ \bar{q}_{t-\zeta_{N(t)}}^{\bar{X}_{\zeta_{N(t)}}^{\pi}} (\bar{X}_{\zeta_{N(t)}}^{\pi}, z) \Gamma_{N(t)}(x) \right], $$

$$ p^K(x, s) = e^{\lambda t} \mathbb{E} \left[ f_{\bar{\tau}}^{\bar{X}_{\zeta_{N(t)}}^{\pi}} (\bar{X}_{\zeta_{N(t)}}^{\pi}, t - \zeta_{N(t)}) \mathbf{I}_{\{s \ge \zeta_{N(t)}\}} \Gamma_{N(t)}(x) + f_{\bar{\tau}}^{\bar{X}_{\zeta_{N(t)-1}}^{\pi}} (\bar{X}_{\zeta_{N(t)-1}}^{\pi}, s - \zeta_{N(t)-1}) \mathbf{I}_{\{s \ge \zeta_{N(t)-1}\}} \bar{\Gamma}_{N(t)}(x) \right]. $$

As a corollary, we also obtain the following integration by parts formula.

**Corollary 6.1.2.** Suppose that the assumptions of the previous theorem hold. Let $h \in C^1((-\infty, L])$ satisfying: there exist $c > 0$, such that for all $z \in (-\infty, L]$, $|h(z)| + |h'(z)| \le c \exp(c|z|)$. Then, for all $x \in (-\infty, L)$, one has

$$ \mathbb{E}[h'(X_t^x) \mathbf{I}_{\{\tau^x > t\}}] = -e^{\lambda T} \mathbb{E} \left[ h(\bar{X}_t^\pi) \Lambda_{T-\zeta_{N(t)}} (\bar{X}_{\zeta_{N(t)}}^\pi, \bar{X}_t^\pi) \mu_{t-\zeta_{N(t)}}^1 (\bar{X}_{\zeta_{N(t)}}^\pi, \bar{X}_t^\pi) \Gamma_{N(t)}(x) \right]. $$

An unbiased Monte Carlo method for evaluating $\mathbb{E}[h(\tau_t^x, X_{\tau_t^x}^x)]$, $\mathbb{E}[h'(X_t^x) \mathbf{I}_{\{\tau^x > t\}}]$ or the probability law $p_T(0, x, dt, dz)$ of $(\tau_t^x, X_{\tau_t^x}^x)$ stems from the probabilistic representations obtained in Theorem 6.1.1 and Corollary 6.1.2. The explosion of the variance may be an important issue that can induce poor convergence rate of the method as pointed out in [4] for unbiased simulation of multi-dimensional diffusions. In these situations, an importance sampling method on the time steps using a Beta or Gamma distribution may be used. In short, it would seem that this approximation will work well in the case of small parameters. A very close analysis could be carried out here but we do not intend to develop importance sampling schemes and refer the interested reader to [4] for some developments in the diffusion case.

## 6.2 Weak uniqueness and density estimates for SDEs with coefficients depending on some path-functionals

In this section, we expose the results obtained in [14] in collaboration with L. Li related to the weak existence and uniqueness as well as density estimates of SDEs with dynamics

$$ X_t = x + \int_0^t b(X_s, A_s(X))ds + \int_0^t \sigma(X_s, A_s(X))dW_s, \quad t \in [0, T] \qquad (6.7) $$

where $(W_t)_{t \ge 0}$ stands for a one-dimensional Brownian motion and $(A_t(X))_{t \ge 0}$ is an $\mathbb{R}^{d-1}$-valued functional depending on the path $X$, $d \ge 2$. Some examples include its local and occupation times, its running maximum or minimum,
its first hitting time of a level, its running average, etc. From the point of view of applications, systems of the type (6.7) appear in many fields. Let us mention stochastic Hamiltonian systems where $A_t(X) = \int_0^t F(s, X_s)ds$, see e.g. Soize [118] for a general overview, Talay [124] for convergence to equilibrium or Barucci and al. [13] for an application to the pricing of Asian options. We also mention the work of Forde [42], where the author constructs a weak solution to the SDE (6.7) with $b \equiv 0$ and $A_t(X) = \max_{0 \le s \le t} X_s$ is the running maximum of X and investigates an application in mathematical finance.

In order to prove weak uniqueness for the SDE (6.7), we follow the recent approach of Bass and Perkins [15] for proving uniqueness for the martingale problem and illustrated it in the framework of non-degenerate, non-divergence and time-homogeneous diffusion operators under the assumption that the diffusion matrix is strictly positive and Hölder continuous. It has also been recently extended by Menozzi [98] for a class of multi-dimensional degenerate Kolmogorov equations that is the case of a multi-dimensional path functional $A = (A_t^1, \dots, A_t^N)_{t \ge 0}$ given by: $A_t^1(X) = \int_0^t F_1(X_s, A_s(X))ds$, $A_t^2(X) = \int_0^t F_2(A_s^1, \dots, A_s^N)ds$, $\dots, A_t^N(X) = \int_0^t F_N(A_s^{N-1}, A_s^N)ds$, under an assumption of weak Hörmander type on the functions $(F_1, \dots, F_N)$. The main new feature added here compared to previous works on this topic is that our technique enables us to deal with a process whose probability law is absolutely continuous with respect to a $\sigma$-finite measure.

For instance, if one considers the couple $(X_t, A_t(X))_{t \ge 0}$, $A_t(X) = L_t^0(X)$ being the symmetric local time at point 0 accumulated by $X$ up to time $t$, it is easy to see that on $\{T_0 > t\}$, $T_0$ being the first hitting time of 0 by $X$, one has $L_t^0(X) = 0$ whereas on $\{T_0 \le t\}$, the process may accumulate local time so that the probability law of the couple $(X_t, L_t^0(X))$ consists in two parts, one being singular with an atom in the local time part, the other one (hopefully) being absolutely continuous with respect to the Lebesgue measure on $\mathbb{R} \times \mathbb{R}_+$. Hence, we see that for such dynamics the situation is more challenging than in the standard diffusion setting. This new difficulty will be overcome by choosing two independent parametrix kernels, one for each part, and then by proceeding with a non-trivial selection of the freezing point according to the singular measure induced by the approximation process.

We first provide a general framework in order to establish weak uniqueness for SDEs with dynamics (6.7). The contribution of this part compared to the existing literature on this topic lies in the fact that we identify the main assumptions, namely (H1) and (H2) see below, needed to establish weak uniqueness, thus allowing for a general treatment.

We assume that there exists a weak solution $(X, W)$, $\{\mathcal{F}_t\}$ to (6.7) and that the process $Y := (X_t, A_t(X))_{t \ge 0}$, starting from the initial point $x$ at time 0, lives on a closed space $\mathcal{J} \subset \mathbb{R}^d$. The process $Y$ induces a probability measure $\mathbb{P}^x$ (or simply denoted $\mathbb{P}$) on $\Omega = \mathcal{C}([0, \infty), \mathcal{J})$ which is endowed with the canonical filtration $(\mathcal{F}_t)_{t \ge 0}$. We consider the collection of linear maps $(P_t)_{t \ge 0}$ defined by $P_t f(x) := \mathbb{E}[f(X_t, A_t(X))]$ for $f \in \mathcal{B}_b(\mathcal{J})$.

We denote the approximation or the proxy process by $\bar{X}$, that is the solution of (6.7) with $b = 0$ and the diffusion coefficient $\sigma$ evaluated at some fixed point $z \in \mathcal{J}$. Without going into details at this point, the key idea is to consider the process $Y$ as a perturbation of the proxy $\bar{Y} = (\bar{X}_t, A(\bar{X}_t))_{t \ge 0}$ whose law is denoted by $\bar{p}_t(x, dy) = \bar{p}_t^z(x, dy)$. Accordingly, we define the collection of linear maps $(\bar{P}_t)_{t \ge 0}$ by $\bar{P}_t f(x) := \mathbb{E}[f(\bar{X}_t, A_t(\bar{X}))] = \int f(y)\bar{p}_t(x, dy)$ for $f \in \mathcal{B}_b(\mathcal{J})$.

We work under the following set of assumptions:

**Assumptions (H1):** Given the initial and frozen point $x, z \in \mathcal{J}$.

(i) (a) The proxy process $\bar{Y}^z$ is a Markov process with infinitesimal generator $\bar{\mathcal{L}}^z$.

(b) There exists a $\sigma$-finite measure $\nu(x, \cdot)$ such that for all $t > 0$, the law of $\bar{Y}_t^z$ is absolutely continuous with respect to $\nu(x, \cdot)$. More specifically, there exists a $\nu(x, dy)$-integrable function $(t, x, y) \mapsto \bar{p}_t^z(x, y)$ satisfying

$$ \bar{p}_t^z(x, dy) = \bar{p}_t^z(x, y)\nu(x, dy) \qquad (6.8) $$

and $\bar{P}_t^z f(x) = \int f(y)\bar{p}_t^z(x, y)\nu(x, dy)$ for all $f \in \mathcal{B}_b(\mathcal{J})$.

(ii) There exists a class of functions $\mathcal{D} \subset \text{Dom}(\bar{\mathcal{L}}^z) \cap \mathcal{C}_b(\mathcal{J})$ and a linear operator $\mathcal{L}$ acting on $\mathcal{D}$ such that:

(a) For all $g \in \mathcal{C}_b^\infty(\mathcal{J})$, $\bar{P}_t^z g \in \mathcal{D}$.

(b) For all functions $h$ such that: $s \mapsto h(s, \cdot) \in \mathcal{C}^1(\mathbb{R}_+, \mathcal{D})$, the process

$$ h(t, Y_t) - h(0, x) - \int_0^t \{\partial_1 h(s, Y_s) + \mathcal{L}h(s, Y_s)\} ds, \quad t \ge 0 $$

is a continuous square integrable martingale under $\mathbb{P}^x$.

(c) There exists a parametrix kernel $\bar{\theta}_t$ with respect to the measure $\nu$, that is a measurable map $(t, z, x, y) \mapsto \bar{\theta}_t^z(x, y)$ such that for all $g \in \mathcal{C}_b^\infty(\mathcal{J})$

$$ (\mathcal{L} - \bar{\mathcal{L}}^z) \bar{P}_t^z g(x) = \int g(y) \bar{\theta}_t^z(x, y) \nu(x, dy), \quad t > 0. \qquad (6.9) $$
(iii) For all $x, y \in \mathcal{J}$, the maps $(t, z) \mapsto \bar{p}_t^z(x, y)$ and $(t, z) \mapsto \bar{\theta}_t^z(x, y)$ are continuous on $(0, \infty) \times \mathcal{J}$.

(iv) For all $t > 0$, there exists some $\nu(x, dy)$-integrable functions $p_t^*(x, y)$, $\theta_t^*(x, y)$, a constant $\bar{\zeta} = \bar{\zeta}(z) \in \mathbb{R}$ and a positive constant $C$, eventually depending on $t$ but in a non-decreasing way, such that

$$ |\bar{p}_t^z(x, y)| \le p_t^*(x, y), \quad |\bar{\theta}_t^z(x, y)| \le \theta_t^*(x, y), \quad \text{and} \quad \int |\theta_t^*(x, y)|\nu(x, dy) \le Ct^{\bar{\zeta}}. \qquad (6.10) $$

For the case $z = y$, we assume that the parametrix kernel enjoys the following smoothing property: there exists $\zeta > -1$ and a positive constant $C$, eventually depending on $t$ in a non-decreasing way, such that

$$ \forall t > 0, \forall x \in \mathcal{J}, \quad \int |\bar{\theta}_t^y(x, y)|\nu(x, dy) \le Ct^{\zeta}. \qquad (6.11) $$

(v) For any $g \in C_b(\mathcal{J})$, one has

$$ \lim_{t \downarrow 0} \int g(y) \bar{p}_t^y(x, y) \nu(x, dy) = g(x). $$

For notational simplicity, we define for $t > 0$,

$$ \hat{p}_t(x, y) := \bar{p}_t^y(x, y), $$

$$ \hat{P}_t f(x) := \int f(y) \hat{p}_t(x, y) \nu(x, dy) = \int f(y) \bar{p}_t^y(x, y) \nu(x, dy), $$

$$ S_t g(x) := \int g(y) \hat{\theta}_t(x, y) \nu(x, dy) = \int g(y) \bar{\theta}_t^y(x, y) \nu(x, dy). $$

**Theorem 6.2.1.** Let $T > 0$. Assume that (H1) holds. Then, for any $g \in C_b(\mathcal{J})$ and any $t \in [0, T]$, one has

$$ P_t g(x) = \hat{P}_t g(x) + \int_0^t P_s S_{t-s} g(x) \, ds. $$

Moreover, if for $t > 0$, the function $x \mapsto S_t g$, $g \in C_b(\mathcal{J})$, belongs to $C_b(\mathcal{J})$. Then one may iterate the first order formula

$$ P_t g(x) = \hat{P}_t g(x) + \sum_{n \ge 1} I_t^n g(x), \quad \text{with } I_t^n g(x) = \int_{\Delta_n(t)} ds_n \hat{P}_{s_n} S_{s_{n-1}-s_n} \cdots S_{t-s_1} g(x) \qquad (6.12) $$

Moreover, the series (6.12) converges absolutely and uniformly for $x \in \mathcal{J}$. The notation $\Delta_n(t)$ has been defined in the previous section, just after equation (6.4).

We consider two weak solutions of the SDE (6.7) starting at time 0 from the same initial point $x \in \mathcal{J}$. Denote by $\mathbb{P}_1$ and $\mathbb{P}_2$ the two probability measures induced on the space $(C([0, \infty), \mathcal{J}), \mathcal{B}(C([0, \infty), \mathcal{J})))$. Define for $f \in B_b(\mathcal{J})$, $\lambda > 0$ and $i = 1, 2$

$$ S_{\lambda}^{i} f(x) := \mathbb{E}_{i} \left[ \int_{0}^{\infty} e^{-\lambda t} f(Y_{t}) dt \right] = \int_{0}^{\infty} e^{-\lambda t} \mathbb{E}_{i} [f(Y_{t})] dt, \quad S_{\lambda}^{\Delta} f(x) := (S_{\lambda}^{1} - S_{\lambda}^{2}) f(x), \quad \|S_{\lambda}^{\Delta}\| := \sup_{\|f\|_{\infty} \le 1} |S_{\lambda}^{\Delta} f|. $$

We notice that by (H1) (ii) b), $S_{\lambda}^i(\lambda - \mathcal{L})f(x) = f(x)$, $\forall f \in \mathcal{D}$, $i = 1, 2$. For $z \in \mathcal{J}$, the resolvent of the process with frozen coefficients is defined by

$$ \bar{R}_{\lambda} f(x) = \int_{0}^{\infty} e^{-\lambda t} \bar{P}_{t} f(x) dt, \quad \forall f \in \mathcal{B}_{b}(\mathcal{J}) \qquad (6.13) $$

and for $f \in \mathcal{D}$ one has $\bar{R}_{\lambda}(\lambda - \mathcal{L})f = f$. We make the following assumptions:

**Assumption (H2):** For all $\lambda > 0$, one has $\bar{R}_{\lambda} \mathcal{D} \subset \mathcal{D}$ and for $f \in \mathcal{D}$,

$$ (\lambda - \bar{\mathcal{L}}) \bar{R}_{\lambda} f = \bar{R}_{\lambda} (\lambda - \bar{\mathcal{L}}) f \quad \text{and} \quad (\mathcal{L} - \bar{\mathcal{L}}) \bar{R}_{\lambda} f(x) = \int_{0}^{\infty} e^{-\lambda t} (\mathcal{L} - \bar{\mathcal{L}}) \bar{P}_{t} f(x) \, dt. \qquad (6.14) $$

**Theorem 6.2.2.** Assume (H1) and (H2) are satisfied, then weak uniqueness holds for the SDE (6.7).
We sketch the proof and omit some technical details. Let $z \in \mathcal{J}$ and $r > 0$. We consider a sequence of non-negative mollifiers $\delta_{\varepsilon}^z$, $\varepsilon > 0$, converging to the Dirac mass at $z$ as $\varepsilon \downarrow 0$. Let us first observe that if $g \in C_b^{\infty}(\mathcal{J})$ then $\bar{P}_r g \in \mathcal{D}$, and by (H2) one has

$$ (\lambda - \mathcal{L}) \bar{R}_{\lambda} \bar{P}_{r} \delta_{\varepsilon}^{z} g(x) = (\lambda - \bar{\mathcal{L}}) \bar{R}_{\lambda} \bar{P}_{r} \delta_{\varepsilon}^{z} g(x) - (\mathcal{L} - \bar{\mathcal{L}}) \bar{R}_{\lambda} \bar{P}_{r} \delta_{\varepsilon}^{z} g(x) = \bar{P}_{r} \delta_{\varepsilon}^{z} g(x) - (\mathcal{L} - \bar{\mathcal{L}}) \bar{R}_{\lambda} \bar{P}_{r} \delta_{\varepsilon}^{z} g(x). \quad (6.15) $$

We integrate both sides of (6.15) with respect to $dz$ and apply $S_{\lambda}^{\Delta}$. For $i = 1, 2$, applying Fubini's theorem, one gets

$$ S_{\lambda}^{i} \int (\lambda - \mathcal{L}) \bar{R}_{\lambda}^{z} \bar{P}_{r}^{z} \delta_{\varepsilon}^{z} g \, dz = \int_{0}^{\infty} e^{-\lambda t} \mathbb{E}_{i} \left[ \int (\lambda - \mathcal{L}) \bar{R}_{\lambda}^{z} \bar{P}_{r}^{z} \delta_{\varepsilon}^{z} g(Y_t) \, dz \right] dt = \int S_{\lambda}^{i} (\lambda - \mathcal{L}) \bar{R}_{\lambda}^{z} \bar{P}_{r}^{z} \delta_{\varepsilon}^{z} g \, dz \\ \text{so that} \\ S_{\lambda}^{\Delta} \left( \int \bar{P}_{r}^{z} \delta_{\varepsilon}^{z} g \, dz \right) - S_{\lambda}^{\Delta} \left( \int (\mathcal{L} - \bar{\mathcal{L}})^{z} \bar{R}_{\lambda}^{z} \bar{P}_{r}^{z} \delta_{\varepsilon}^{z} g \, dz \right) = 0. \quad (6.16) $$

Passing to the limit as $\varepsilon \downarrow 0$ and $r \downarrow 0$, one gets $S_{\lambda}^{\Delta} g = S_{\lambda}^{\Delta} \left( \int_{0}^{\infty} e^{-\lambda t} S_t g \, dt \right)$ and one can pick $\lambda$ such that $|\int_{0}^{\infty} e^{-\lambda t} S_t g \, dt| \le |g|_{\infty} \int_{0}^{\infty} e^{-\lambda t} t \zeta \, dt = |g|_{\infty} \frac{\Gamma(\zeta)}{\lambda^{1+\zeta}} < \frac{1}{2}|g|_{\infty}$. Hence, from the definition of $\|S_{\lambda}^{\Delta}\|$, we find that $|S_{\lambda}^{\Delta} g| = |S_{\lambda}^{\Delta} \left( \int_{0}^{\infty} e^{-\lambda t} S_t g \, dt \right)| \le \frac{1}{2} \|S_{\lambda}^{\Delta}\| \|g|_{\infty}\$. Taking the supremum over $|g|_{\infty} \le 1$ yields $\|S_{\lambda}^{\Delta}\| \le \frac{1}{2} \|S_{\lambda}^{\Delta}\|$ and, since $\|S_{\lambda}^{\Delta}\| < \infty$, we conclude that $S_{\lambda}^{\Delta} = 0$. By the uniqueness of the Laplace transform together with continuity w.r.t the variable $t$, $\mathbb{E}_1[g(Y_t)] = \mathbb{E}_2[g(Y_t)]$ for all $t \ge 0$ if $g$ is bounded measurable. Now one can use the standard argument based on regular conditional probabilities to show that the finite dimensional distributions of the process $(Y_t)_{t \ge 0} = (X_t, A_t(X))_{t \ge 0}$ agree under $\mathbb{P}_1$ and $\mathbb{P}_2$.

In [14], we considered the two following examples: $A_t(X) = L_t^0(X)$ where $L_t^0(X)$ is the symmetric local time at $0$ accumulated by $X$ at time $t$ with state space $\mathcal{J} = \mathbb{R} \times \mathbb{R}_+$ and $A_t(X) = M_t := m_0 \vee \max_{0 \le s \le t} X_s$, $m_0 \ge x_0$, with state space $\mathcal{J} = \{(x, m) \in \mathbb{R}^2 : x \le m\}$. We work under the following mild assumptions:

(R-η) The coefficients $b$ and $a = \sigma^2$ are bounded measurable functions defined on $\mathcal{J}$. The diffusion coefficient $a$ is $\eta$-Hölder continuous on $\mathcal{J}$.

(UE) There exists some constant $a > 0$ such that $\forall z \in \mathcal{J}, a \le a(z)$.

In the following we detail the computations for the first case only, namely the diffusion process with its symmetric local time at zero. A similar analysis for the SDE with its running maximum is done in [14]. The set $\mathcal{D}$ is the class of function $f \in C_b^{2,1}(\mathbb{R}\setminus\{0\} \times \mathbb{R}_+) \cap C_b(\mathbb{R} \times \mathbb{R}_+)$ such that $\partial_1 f(0+, \ell) = \lim_{x \downarrow 0} \frac{f(x, \ell) - f(0, \ell)}{x}$ and $\partial_1 f(0-, \ell) = \lim_{x \uparrow 0} \frac{f(x, \ell) - f(0, \ell)}{x}$ exist, are finite and satisfy the following transmission condition:

$$ \forall \ell \in \mathbb{R}_+, \quad \frac{\partial_1 f(0+, \ell) - \partial_1 f(0-, \ell)}{2} + \partial_2 f(0, \ell) = 0. \quad (6.17) $$

We define the linear operator $\mathcal{L}$ by $\mathcal{L}f(x, \ell) = b(x, \ell)\partial_1 f(x-, \ell) + \frac{1}{2}a(x, \ell)\partial_1^2 f(x-, \ell)$, $(x, \ell) \in \mathbb{R} \times \mathbb{R}_+$. As we need a chain rule formula for the process $(X_t, A_t(X))_{t \ge 0}$, we rely on the following generalisation of Itô's Lemma:

$$ \begin{aligned} f(t, X_t, L_t^0(X)) = f(0, x, 0) &+ \int_0^t \{\partial_1 f(s, X_s, L_s^0(X)) + \mathcal{L}f(s, \cdot)(X_s, L_s^0(X))\} ds \\ &+ \int_0^t \left\{ \frac{\partial_2 f(s, 0+, L_s^0(X)) - \partial_2 f(s, 0-, L_s^0(X))}{2} + \partial_3 f(s, 0, L_s^0(X)) \right\} dL_s^0(X) \\ &+ \int_0^t \sigma(X_s, L_s^0(X)) \partial_2 f(s, X_s-, L_s^0(X)) dW_s \quad \text{a.s.} \end{aligned} $$

for $f \in C^{1,2,1}(\mathbb{R}_+ \times \mathbb{R}\setminus\{0\} \times \mathbb{R}_+) \cap C(\mathbb{R}_+ \times \mathbb{R} \times \mathbb{R}_+)$ such that $\partial_2 f(t, 0+, \ell) = \lim_{x \downarrow 0} (f(t, x, \ell) - f(t, 0, \ell))/x$ and $\partial_2 f(t, 0-, \ell) = \lim_{x \uparrow 0} (f(t, x, \ell) - f(t, 0, \ell))/x$ exist and are finite. This chain rule formula allows one to associate a martingale problem to the SDE (6.7) with $A_t(X) = L_t^0(X)$. This connection provides a natural way to establish weak existence by a compactness argument when $b$ and $\sigma$ are continuous, see e.g. Stroock and Varadhan [120]. Remark that one can remove the continuity assumption on $b$ when $\sigma$ is uniformly elliptic via the Girsanov transform method.

We now introduce the proxy process $\bar{X}_t := x_0 + \sigma(z_1)W_t$, $t \ge 0$, $z_1 = (x_1, \ell_1) \in \mathbb{R} \times \mathbb{R}_+$. For $f \in C_b(\mathbb{R} \times \mathbb{R}_+)$, simple computations show that $\bar{P}_t f(x_0, \ell_0) = \mathbb{E}[f(\bar{X}_t, \ell_0 + L_t^0(\bar{X}))] = \int_{\mathbb{R} \times \mathbb{R}_+} f(x, \ell) \bar{p}_t(x_0, \ell_0, x, \ell) \nu(x_0, \ell_0, dx, d\ell)$ with $\bar{p}_t(x_0, \ell_0, x, \ell) := \bar{f}_t(x_0, x) \mathbf{I}_{\{\ell=\ell_0\}} + \bar{q}_t(x_0, \ell_0, x, \ell) \mathbf{I}_{\{\ell_0<\ell\}}$ and

$$ \begin{aligned} \bar{f}_t(x_0, x) &:= g(a(z_1)t, x - x_0) - g(a(z_1)t, x + x_0), \\ \bar{q}_t(x_0, \ell_0, x, \ell) &:= \frac{(|x| + |x_0| + \ell - \ell_0)}{a^{\frac{3}{2}}(z_1)t} g(a(z_1)t, |x| + |x_0| + \ell - \ell_0), \\ \nu(x_0, \ell_0, dx, d\ell) &:= \mathbf{I}_{\{\ell_0<\ell\}} dxd\ell + \mathbf{I}_{\{x_0x \ge 0\}} dx\delta_{\ell_0}(d\ell). \end{aligned} $$
and $\hat{f}_t(x_0, x) := g(a(x, \ell_0)t, x - x_0) - g(a(x, \ell_0)t, x + x_0)$, $\hat{q}_t(x_0, \ell_0, x, \ell) := \frac{(|x|+|x_0|+\ell-\ell_0)}{a^{\frac{1}{2}}(x, \ell)t} g(a(x, \ell)t, |x| + |x_0| + \ell - \ell_0)$.

We also observe that the measure $\nu(x_0, \ell_0, dx, d\ell)$ satisfies an important convolution type property in the sense that

$$ \nu(x_0, \ell_0, dx', d\ell')\nu(x', \ell', dx, d\ell) = u(x_0, \ell_0, x, \ell, dx', d\ell')\nu(x_0, \ell_0, dx, d\ell) \quad (6.18) $$

As an application of the general methodology developed above, we obtain the following result.

**Theorem 6.2.3.** If (R-η) and (UE) are satisfied for some $\eta \in (0, 1]$ then weak uniqueness holds for the SDE (6.7) with $A_t(X) = L_t^0(X)$, $t \ge 0$.

Moreover, the first step expansion obtained in Theorem 6.2.1 holds. However, since $(x_0, \ell_0) \mapsto \partial_{x_0} \bar{q}_t(x_0, \ell_0, x, \ell)$ is not continuous at zero, the function $(x_0, \ell_0) \mapsto S_t g(x_0, \ell_0)$ shares the same property unless the drift $b$ satisfies $b(0, \ell) = 0$, $\ell \in \mathbb{R}_+$. Unfortunately, one cannot iterate the first step expansion. We first remove the drift. Now the function $S_t g$ is continuous. Applying repeatedly (6.18) and Fubini's theorem, one gets

$$ \int_{\Delta_n(T)} d\mathbf{s}_n \hat{P}_{\mathbf{s}_n} \mathcal{S}_{s_{n-1}-s_n} \dots \mathcal{S}_{t-s_1} g(x_0, \ell_0) = \int_{\mathbb{R} \times \mathbb{R}_+} g(x, \ell) p_t^n(x_0, \ell_0, x, \ell) \nu(x_0, \ell_0, dx, d\ell) $$

with

$$ p_t^n(x_0, \ell_0, x, \ell) := \begin{cases} \int_{\Delta_n(t)} d\mathbf{s}_n \int_{(\mathbb{R} \times \mathbb{R}_+)^n} \hat{p}_{\mathbf{s}_n}(x_0, \ell_0, x_1, \ell_1) \times \left[ \prod_{i=1}^n \hat{\theta}_{s_{n-i}-s_{n-i+1}}(x_i, \ell_i, x_{i+1}, \ell_{i+1}) u(x_0, \ell_0, x_{i+1}, \ell_{i+1}, dx_i, d\ell_i) \right] & n \ge 1, \\ \hat{p}_t(x_0, \ell_0, x, \ell) & n = 0 \end{cases} \quad (6.19) $$

where we omit the exact expression of $\hat{\theta}_t(x_0, \ell_0, x, \ell)$.

As our second main result, we prove that the transition density of $(X_t, \ell_0 + L_t^0(X))_{t \ge 0}$ exists and satisfies a Gaussian upper bound.

**Theorem 6.2.4.** Let $T > 0$. Assume that (R-η) and (UE) hold for some $\eta \in (0, 1]$. For $(t, x_0, \ell_0) \in (0, T] \times \mathbb{R} \times \mathbb{R}_+$, define the probability measure

$$ p_t(x_0, \ell_0, dx, d\ell) := p_t(x_0, \ell_0, x, \ell)\nu(x_0, \ell_0, dx, d\ell) = p_t^1(x_0, \ell_0, x, \ell) \mathbf{I}_{\{\ell_0 < \ell\}} dx d\ell + p_t^2(x_0, \ell_0, x) \mathbf{I}_{\{x_{\ell_0} \ge 0\}} dx \delta_{\ell_0}(d\ell) $$

with $p_t(x_0, \ell_0, x, \ell) := \sum_{n \ge 0} p_t^n(x_0, \ell_0, x, \ell)$, $p_t^1(x_0, \ell_0, x, \ell) := \sum_{n \ge 0} p_t^n(x_0, \ell_0, x, \ell)$, $p_t^2(x_0, \ell_0, x) := \sum_{n \ge 0} p_t^n(x_0, \ell_0, x, \ell_0)$.

Then, both series defining $p_t^1(x_0, \ell_0, x, \ell)$ and $p_t^2(x_0, \ell_0, x)$ converge absolutely and uniformly for $(x_0, \ell_0)$, $(x, \ell) \in (\mathbb{R} \times \mathbb{R}_+)^2$. Moreover for $h \in C_b(\mathbb{R} \times \mathbb{R}_+)$ the following representation for the semigroup holds,

$$ P_t h(x_0, \ell_0) = \int_{\mathbb{R} \times \mathbb{R}_+} h(x, \ell) p_t(x_0, \ell_0, x, \ell) \nu(x_0, \ell_0, dx, d\ell). $$

Therefore, for all $(x_0, \ell_0) \in \mathbb{R} \times \mathbb{R}_+$, the function $(x, \ell) \mapsto p_t(x_0, \ell_0, x, \ell)$ is the probability density function of the random vector $(X_t^{x_0}, \ell_0 + L_t^0(X^{x_0}))$ with respect to the $\sigma$-finite measure $\nu(x_0, \ell_0, dx, d\ell)$, where $X_t^{x_0}$ is the solution taken at time $t$ of the SDE (6.7) (with $A_t(X) = L_t^0(X)$) starting from $x_0$ at time 0, $L_t^0(X^{x_0})$ being its running symmetric local time at time $t$.

Finally, there exists some constants $C, c > 1$ such that for all $(x_0, \ell_0)$, $(x, \ell) \in \mathbb{R} \times \mathbb{R}_+$, the following Gaussian upper-bounds hold

$$ p_t^1(x_0, \ell_0, x, \ell) \le C t^{-1/2} g(ct, |x| + |x_0| + \ell - \ell_0) \quad \text{and} \quad p_t^2(x_0, \ell_0, x) \le C g(ct, x - x_0). \quad (6.20) $$

The proof is done in several steps. First, we assume that $b = 0$ and obtain the above series expansion. Then, we prove the existence of a transition density in the presence of a bounded measurable drift via Girsanov transform. This allows us to extend Theorem 6.2.1 for a non-continuous function $g$. Finally we repeat the same argument as above.

We point out that the proof of the convergence of the asymptotic expansion for the transition density is not standard in the current setting. Indeed, in the classical diffusion setting, the parametrix expansion of the transition density converges since the order of the singularity in time induced by the parametrix kernel $\hat{\theta}_t$ is of order $t^{-1+\eta/2}$, which is still integrable near 0. The main difficulty that appears here is when one wants to control the whole convolution appearing in the right-hand side of (6.19). More precisely, it lies in the cross-terms which are of a different nature, for instance when one convolutes the non-singular part in the convolution kernel $\hat{\theta}_{T-s_1}$ with the singular part in the convolution kernel $\hat{\theta}_{s_1-s_2}$. Standard arguments similar to those used in the standard diffusion case do not guarantee the convergence of the integral defining (6.19). To overcome this difficulty and show that the parametrix expansion for the transition density converges, one has to make use of the symmetry in the initial and terminal point of the density of the killed proxy process, in order to retrieve the integrability in time of the underlying convolution kernel.
## 6.3 Perspective

In [13], we considered probabilistic representations for the first hitting time and the associated killed process of a one dimensional elliptic diffusion. We obtain some integration by parts for the killed diffusion process for the first derivative of a test function. One may be interested in establishing higher order integration by parts and also integration by parts formulas for the first hitting time. This is not a simple extension of Corollary (6.1.2) where we performed an integration by parts in the last interval $[\zeta_N, T]$. If we want to keep the Malliavin weights in $L^2(\mathbb{P})$, one has to perform the integration by parts on the whole time interval $[0, T]$ using the whole underlying noise.

One may be interested in establishing a probabilistic representation as well as integration by parts formulas for other path functionals such as the local time, the running maximum, the occupation time, etc as considered in Section 6.2. Using the parametrix expansions of the transition densities of $(X_t, A_t(X))_{t \ge 0}$, one can prove that the corresponding semigroup $P_t h$ is smooth when $a$ and $b$ are Hölder continuous. Then, one may proceed using the forward parametrix method as presented in Section 6.1.

More challenging extensions could include other type of processes. One may notably consider the case of a skew diffusion with path-dependent coefficients involving its local and occupation times, see Appuhamillage & al. [5] for an expression of the trivariate density $(B_t^{(\alpha)}, L_t^0(B^{(\alpha)}), \Gamma_t^0(B^{(\alpha)}))$, $t \ge 0$, where $(B_t^{(\alpha)})_{t \ge 0}$ is an $\alpha$-skew Brownian motion or reflected SDEs. We are working on this topic with L. Li.

Another interesting point to investigate in that direction is to weaken the assumption on the coefficients, notably on $a$ by assuming only uniform ellipticity and continuity, in order to establish weak uniqueness for the martingale problem using Calderón-Zygmund estimates, see Menozzi [99] for some degenerate equations of Kolmogorov type. Another related issue concerns pathwise uniqueness for the SDE (6.7). A possible strategy could be to extend the results of Zvonkin [128] and Veretennikov [125], Krylov and Röckner [79] established in the diffusion case to the case of a diffusion with a path-functional.

On may study a class of Hörmander equations with non-regular coefficients. The structure of the irregularity has to be carefully chosen so that the approximation process satisfies a *uniform Hörmander condition* or even a *UFG condition*. We refer to the work of Kusuoka and Stroock [84] for a careful analysis of the regularity properties of the unique solution to the related linear parabolic PDE by probabilistic tools. We also refer to the work of Crisan and Delarue [32] for an extension to degenerate semi-linear PDE. For instance, one may consider an equation of the form: $V_0 + (1/2) \sum_{i=1}^N V_i^2, V_0, \dots, V_N$ being first order differential operators and the drift term $V_0$ being irregular. One simple example is $V_0 + \sum_{i=1}^N \varphi_i V_i + \sum_{\alpha \in A_0(2)} \varphi_\alpha V_{[\alpha]}$, where $V_0 \in C_\infty^\infty(\mathbb{R}^d)$ and $(\varphi_i)_{1 \le i \le N}$ are bounded measurable functions, $(\varphi_\alpha)_{\alpha \in A_0(2)}$ are $C_b^1(\mathbb{R}^d)$. Questions of interests are: pathwise uniqueness, density estimates, regularity of the solution of the PDE, etc.

A challenging problem is to investigate probabilistic representation such as the one obtained in Section 6.1 for McKean-Vlasov equations. These equations exhibit a dependence w.r.t to the law of the process in its coefficients. Unlike the two examples investigated in [14], the law of the underlying process is a deterministic object so that it is not natural to expect a regularising property w.r.t. this variable. We refer to the P.h.d. of Murray [96] and the recent work of Crisan and Murray [33] for some integration by parts formula and density estimates for such equations.
Bibliography

[1] M. Ben ALAYA et A. KEBAIER : Central limit theorem for the multilevel monte carlo euler method and applications to asian options. *Ann. Appl. Probab.*, 2012.

[2] A. ALFONSI, B. JOURDAIN et A. KOHATSU-HIGA : Pathwise optimal transport bounds between a one-dimensional diffusion and its Euler scheme. *Ann. Appl. Probab.*, 24(3):1049–1080, 2014.

[3] A. ALFONSI, B. JOURDAIN et A. KOHATSU-HIGA : Optimal transport bounds between the time-marginals of a multidimensional diffusion and its Euler scheme. *Electron. J. Probab.*, 20:no. 70, 31, 2015.

[4] P. ANDERSSON et A. KOHATSU-HIGA : Unbiased simulation of stochastic differential equations using parametrix expansions. *Bernoulli*, 23(3):2028–2057, 08 2017.

[5] T. APPUHAMILLAGE, V. BOKIL, E. THOMANN, E. WAYMIRE et B. WOOD : Occupation and local times for skew brownian motion with applications to dispersion across an interface. *Ann. Appl. Probab.*, 21(1):183–214, 02 2011.

[6] R. AZENCOTT : Densité des diffusions en temps petit: développements asymptotiques. I. *In Seminar on probability, XVIII*, volume 1059 de *Lecture Notes in Math.*, pages 402–498. Springer, Berlin, 1984.

[7] V. BALLY et Arturo KOHATSU-HIGA : A probabilistic interpretation of the parametrix method. *Ann. Appl. Probab.*, 25(6):3095–3138, 12 2015.

[8] V. BALLY et G. PAGÈS : Error analysis of the optimal quantization algorithm for obstacle problems. *Stochastic Process. Appl.*, 106(1):1–40, 2003.

[9] V. BALLY et C. REY : Approximation of markov semigroups in total variation distance. *Electron. J. Probab.*, 21:44 pp., 2016.

[10] V. BALLY et D. TALAY : The law of the Euler scheme for stochastic differential equations: I. Convergence rate of the distribution function. *Prob. Th. Rel. Fields*, 104-1:43–60, 1996.

[11] V. BALLY et D. TALAY : The law of the Euler scheme for stochastic differential equations, II. Convergence rate of the density. *Monte-Carlo methods and Appl.*, 2:93–128, 1996.

[12] F. BARTHE et C. BORDENAVE : Combinatorial optimization over two random point sets. *To appear in Séminaire de Probabilités*, 2011.

[13] E. BARUCCI, S. POLIDORO et V. VESPRI : Some results on partial differential equations and Asian options. *Math. Models Methods Appl. Sci.*, 11(3):475–497, 2001.

[14] R. F. BASS : *Diffusions and Elliptic Operators*. Springer, 1997.

[15] R.F. BASS et E.A. PERKINS : A new technique for proving uniqueness for martingale problems. *From Probability to Geometry (I): Volume in Honor of the 60th Birthday of Jean-Michel Bismut*, pages 47–53, 2009.

[16] G. BEN AROUS : Développement asymptotique du noyau de la chaleur hypoelliptique hors du cut-locus. *Ann. Sci. École Norm. Sup. (4)*, 21(3):307–331, 1988.

[17] G. BEN AROUS : Methodes de laplace et de la phase stationnaire sur l'espace de wiener. *Stochastics*, 25(3):125–153, 1988.

[18] F.E. BENTH, J.S. BENTH et S. KOEKEBAKKER : *Stochastic Modelling of Electricity and Related Markets*. Advanced series on statistical science & applied probability. World Scientific, 2008.
[19] A. BENVENISTE, M. MÉTIVIER et P. PRIOURET : *Adaptive algorithms and stochastic approximations*, volume 22 de *Applications of Mathematics (New York)*. Springer-Verlag, Berlin, 1990. Translated from the French by Stephen S. Wilson.

[20] A. BESKOS, O PAPASPILIΟPOULOS et Roberts G. O. : Retrospective exact simulation of diffusion sample paths with applications. *Bernoulli*, 12(6):1077-1098, 2006.

[21] J.-M. BISMUT : *Large deviations and the Malliavin calculus*, volume 45 de *Progress in Mathematics*. Birkhäuser Boston, Inc., Boston, MA, 1984.

[22] E. BOISSARD : Simple bounds for the convergence of empirical and occupation measures in 1-Wasserstein distance. *Electron. J. of Probab.*, 16, 2011.

[23] F. BOLLEY, A. GUILLIN et C. VILLANI : Quantitative concentration inequalities for empirical measures on non-compact spaces. *Prob. Th. Rel. Fields*, 137:541-593, 2007.

[24] F. BOLLEY et C. VILLANI : Weighted Csiszár-Kullback-Pinsker inequalities and applications to transportation inequalities. *Annales de la Faculté des Sciences de Toulouse (6)*, 14(3):331-352, 2005.

[25] O. BRANDIÈRE et M. DUFLO : Les algorithmes stochastiques contournent-ils les pièges? *Annales de l'I.H.P.*, 32(3):395-427, 1996.

[26] M. BROADIE et O KAYA : Exact simulation of stochastic volatility and other affine jump diffusion processes. *Operations Research*, 54(2):217-231, 2006.

[27] N. BRUTI-LIBERATI et E. PLATEN : Strong approximations of stochastic differential equations with jumps. *Journal of Computational and Applied Mathematics*, 205(2):982-1001, 2007. Special issue on evolutionary problems.

[28] P. CATTIAUX : Calcul stochastique et opérateurs dégénérés du second ordre. II. Problème de Dirichlet. *Bull. Sci. Math.*, 115(1):81-122, 1991.

[29] P. E. Chaudru de RAYNAL : Strong existence and uniqueness for degenerate SDE with Hölder drift. *Ann. Inst. Henri Poincaré Probab. Stat.*, 53(1):259-286, 2017.

[30] B. CHEN, C. W. OOSTERLEE et H. VAN DER WEIDE : A low-bias simulation scheme for the sabr stochastic volatility model. *International Journal of Theoretical and Applied Finance*, 15(02):1250016, 2012.

[31] J.-M. CORON : *Control and nonlinearity*, volume 136 de *Mathematical Surveys and Monographs*. American Mathematical Society, Providence, RI, 2007.

[32] D. CRISAN et F. DELARUE : Sharp gradient bounds for solutions of degenerate semi-linear partial differential equations. *Journal of functional analysis*, 263(10):3024-3101, 2012.

[33] D. CRISAN et E. MCMURRAY : Smoothing properties of mckean-vasov sdes. *Probab. Theory Relat. Fields*, 2017.

[34] F. DELARUE et S. MENOZZI : Density estimates for a random noise propagating through a chain of differential equations. *J. Funct. Anal.*, 259(6):1577-1630, 2010.

[35] C. DELLACHERIE et P.-A. MEYER : *Probabilités et potentiel*. Hermann, Paris, 1975. Chapitres I à IV, Édition entièrement refondue, Publications de l'Institut de Mathématique de l'Université de Strasbourg, No. XV, Actualités Scientifiques et Industrielles, No. 1372.

[36] S. DEREICH : Multilevel Monte Carlo algorithms for Lévy-driven SDEs with Gaussian correction. *Ann. Appl. Probab.*, 21(1):283-311, 2011.

[37] J. D. DEUSCHEL, P. K. FRIZ, A. JACQUIER et S. VIOLANTE : Marginal density expansions for diffusions and stochastic volatility I: Theoretical foundations. *Comm. Pure Appl. Math.*, 67(1):40-82, 2014.

[38] D DUFFIE et P. GLYNN : Efficient monte carlo simulation of security prices. *Ann. Appl. Probab.*, 5(4):897-905, 1995.

[39] M. DUFLO : *Algorithmes stochastiques*, volume 23 de *Mathématiques & Applications (Berlin) [Mathematics & Applications]*. Springer-Verlag, Berlin, 1996.
[40] P. ÉTORÉ et M. MARTINEZ : Exact simulation for solutions of one-dimensional stochastic differential equations with discontinuous drift. *ESAIM Probab. Stat.*, 18:686-702, 2014.

[41] H. FÖLLMER et A. SCHIED : Convex measures of risk and trading constraints. *Finance Stoch.*, 6(4):429-447, 2002.

[42] M. FORDE : A diffusion-type process with a given joint law for the terminal level and supremum at an independent exponential time. *Stoch. Process. Appl.*, 121(12):2802-2817, 2011.

[43] G. FORT : Central limit theorems for stochastic approximation with controlled markov chain dynamics. *ESAIM: PS*, 19:60-80, 2015.

[44] A. FRIEDMAN : *Partial differential equations of parabolic type*. Prentice-Hall, 1964.

[45] A. FRIEDMAN : *Stochastic differential equations*. Chapmann-Hall, 1975.

[46] M. G. GARRONI et J.-L. MENALDI : *Green functions for second order parabolic integro-differential problems*, volume 275 de *Pitman Research Notes in Mathematics Series*. Longman Scientific & Technical, Harlow; copublished in the United States with John Wiley & Sons, Inc., New York, 1992.

[47] M. B. GILES : Improved multilevel Monte Carlo convergence using the Milstein scheme. *In Monte Carlo and quasi-Monte Carlo methods 2006*, pages 343-358. Springer, Berlin, 2008.

[48] M. B. GILES : Multilevel Monte Carlo path simulation. *Oper. Res.*, 56(3):607-617, 2008.

[49] M. B. GILES, D. J. HIGHAM et X. MAO : Analysing multi-level Monte Carlo for options with non-globally Lipschitz payoff. *Finance Stoch.*, 13(3):403-413, 2009.

[50] E. GOBET : *Schémas d'Euler pour diffusion tuée. Application aux options barrière*. PhD Thesis, University Paris VII., 1998.

[51] E. GOBET : Euler schemes for the weak approximation of killed diffusion. *Stoch. Proc. Appl.*, 87:167-197, 2000.

[52] E. GOBET : Weak approximation of killed diffusion using Euler schemes. *Stoch. Process. Appl.*, 87(2):167-197, 2000.

[53] E. GOBET et C. LABART : Short time asymptotics of the density of the Euler scheme. *Electron. Comm. iProbab.*, 23, 2008.

[54] E. GOBET et S. MENOZZI : Stopped diffusion processes: boundary corrections and overshoot. *Stoch. Process. Appl.*, 2009.

[55] N. GOZLAN et C. LÉONARD : A large deviation approach to some transportation cost inequalities. *Prob. Th. Rel. Fields*, 139(1-2):235-283, 2007.

[56] N. GOZLAN et C. LÉONARD : Transport inequalities. A survey. *Markov Process. Related Fields*, 16(4):635-736, 2010.

[57] S. GRAF et H. LUSCHGY : *Foundations of Quantization for Probability Distributions*. Lect. Notes in Math. 1730, 2000.

[58] J. GUYON : Euler scheme and tempered distributions. *Stoch. Process. Appl.*, 116(6):877-904, 2006.

[59] M. HAIRER, M HUTZENTHALER et A. JENTZEN : Loss of regularity for kolmogorov equations. *Ann. Probab.*, 43(2):468-527, 2015.

[60] J. M. HARRISON et L. A. SHEPP : On skew Brownian motion. *Ann. Probab.*, 9(2):309-313, 1981.

[61] S. HEINRICH : *Multilevel Monte Carlo methods*. Springer, 2001. In Large-scale scientific computing.

[62] P. HENRY-LABORDÈRE, X. TAN et N. TOUZI : A numerical algorithm for a class of BSDEs via the branching process. *Stochastic Process. Appl.*, 124(2):1112-1140, 2014.

[63] I HONORÉ, S. MENOZZI et G. PAGÈS : Non-asymptotic gaussian estimates for the recursive approximation of the invariant measure of a diffusion. *Preprint*, 2016.

[64] N. IKEDA et S. WATANABE : *Stochastic differential equations*. North Holland, 1980.
[65] J. JACOD: Rates of convergence to the local time of a diffusion. Ann. I.H.P., section B, 34:505-544, 1998.

[66] J. JACOD, T. G. KURTZ, S. MÉLÉARD et P. PROTTER: The approximate euler method for Lévy driven stochastic differential equations. Ann. Inst. H. Poincaré Probab. Statist, 41-3:523-558, 2005.

[67] J. JACOD et A.N. SHIRYAYEV: Local martingales and the fundamental asset pricing theorems in the discrete-time case. Finance and Stochastics, 2:259-273, 1998.

[68] B. JOURDAIN et J. LELONG: Robust adaptive importance sampling for normal vectors. Ann. App. Probab., 19(5):1687-1718, 2009.

[69] M. KAINA et L. RÜSCHENDORF: On convex risk measures on $l^p$-spaces. Mathematical Methods in Operations Research (MMR), 69:475-495, 2009.

[70] A. KEBAIER: Statistical Romberg extrapolation: a new variance reduction method and applications to option pricing. Ann. Appl. Probab., 15(4):2681-2705, 2005.

[71] P. E. KLOEDEN et E. PLATEN: *Numerical solution of stochastic differential equations*, volume 23 de *Applications of Mathematics (New York)*. Springer-Verlag, Berlin, 1992.

[72] A. KOHATSU-HIGA et P. PROTTER: The Euler scheme for SDE's driven by semimartingales. *In Stochastic analysis on infinite-dimensional spaces (Baton Rouge, LA, 1994)*, volume 310 de *Pitman Res. Notes Math. Ser.*, pages 141-151. Longman Sci. Tech., Harlow, 1994.

[73] A. KOHATSU-HIGA, D. TAGUCHI et J. ZHONG: The parametrix method for skew diffusions. *Potential Analysis*, pages 1-31, 2016.

[74] A. KOHATSU-HIGA et P. TANKOV: Jump-adapted discretization schemes for lévy-driven sdes. *Stoch. Process. Appl.*, 120:2258-2285, 2010.

[75] V. KONAKOV et E. MAMMEN: Local limit theorems for transition densities of Markov chains converging to diffusions. Prob. Th. Rel. Fields, 117:551-587, 2000.

[76] V. KONAKOV et E. MAMMEN: Edgeworth type expansions for euler schemes for stochastic differential equations. Monte Carlo Methods Appl., 8-3:271-285, 2002.

[77] V. KONAKOV et S. MENOZZI: Weak error for stable driven stochastic differential equations: Expansion of the densities. Journal of Theoretical Probability, 24(2):454-478, 2011.

[78] V. KONAKOV et S. MENOZZI: Weak error for the Euler scheme approximation of diffusions with non-smooth coefficients. Electron. J. Probab., 22:Paper No. 46, 47, 2017.

[79] N.V. KRYLOV et M. RÖCKNER: Strong solutions of stochastic equations with singular time dependent drift. Probab. Th. Rel. fields, 131(2):154-196, 2005.

[80] H. KUNITA: *Stochastic flows and stochastic differential equations*, volume 24 de *Cambridge Studies in Advanced Mathematics*. Cambridge University Press, Cambridge, 1990.

[81] H. J. KUSHNER et G. G. YIN: *Stochastic approximation and recursive algorithms and applications*, volume 35 de *Applications of Mathematics (New York)*. Springer-Verlag, New York, second édition, 2003. Stochastic Modelling and Applied Probability.

[82] S. KUSUOKA: *Approximation of expectation of diffusion processes based on Lie algebra and Malliavin calculus*, pages 69-83. Springer Japan, Tokyo, 2004.

[83] S. KUSUOKA et D. STROOCK: Applications of the Malliavin calculus. II. *J. Fac. Sci. Univ. Tokyo Sect. IA Math*, 32:1-76, 1985.

[84] S. KUSUOKA et D. STROOCK: Applications of the Malliavin calculus. III. *J. Fac. Sci. Univ. Tokyo Sect. IA Math*, 34:391-442, 1987.

[85] O.A. LADYZENSKAJA, V.A. SOLONNIKOV et N.N. URAL'CEVA: *Linear and quasi-linear equations of parabolic type*. Vol.23 Trans. Math. Monog., AMS, Providence, 1968.

[86] S. LARUELLE: Analyse d'algorithmes stochastiques appliqués à la finance, Phd dissertation, 2011.
[87] S. LARUELLE et G. PAGÈS : Stochastic approximation with averaging innovation applied to finance. *Monte Carlo Methods and Applications*, 18(1):1–52, 2012.

[88] J.-F. LE GALL : One-dimensional stochastic differential equations involving the local times of the unknown process. In *Stochastic analysis and applications (Swansea, 1983)*, volume 1095 of *Lecture Notes in Math.*, pages 51–82. Springer, Berlin, 1984.

[89] A. LEJAY : On the constructions of the skew Brownian motion. *Probab. Surv.*, 3:413–466, 2006.

[90] LELONG, JÉRÔME : Asymptotic normality of randomly truncated stochastic algorithms. *ESAIM: PS*, 17:105–119, 2013.

[91] V. LEMAIRE et S. MENOZZI : On some non asymptotic bounds for the euler scheme. *Electron. J. Probab.*, 15:1645–1681, 2010.

[92] V. LEMAIRE et G. PAGÈS : Unconstrained recursive importance sampling. *Ann. Appl. Probab.*, 20(3):1029–1067, 2006.

[93] V. LEMAIRE et G. PAGÈS : Multilevel richardson-romberg extrapolation. *Bernoulli*, 23(4A):2643–2692, 11 2017.

[94] F. MALRIEU et D. TALAY : Concentration inequalities for euler schemes. *Monte Carlo and Quasi-Monte Carlo Methods 2004*, H. Niederreiter and D. Talay (Eds.), pages 355–37, 2006.

[95] H. P. MCKEAN et I. M. SINGER : Curvature and the eigenvalues of the Laplacian. *J. Differential Geometry*, 1:43–69, 1967.

[96] E. MCMURRAY : Regularity of mckean-vlasov stochastic differential equations and applications, phd dissertation, 2015.

[97] S. MENOZZI : Parametrix techniques and martingale problems for some degenerate Kolmogorov equations. *Electron. Commun. Probab.*, 16:234–250, 2011.

[98] S. MENOZZI : Parametrix techniques and martingale problems for some degenerate Kolmogorov equations. *Electron. Commun. Probab.*, 16:234–250, 2011.

[99] S. MENOZZI : Martingale problems for some degenerate kolmogorov equations. *Preprint*, 2015.

[100] R. MIKULEVICIUS : On the rate of convergence of simple and jump-adapted weak Euler schemes for Lévy driven SDEs. *Stochastic Process. Appl.*, 122(7):2730–2757, 2012.

[101] R. MIKULEVIČIUS et E. PLATEN : Rate of convergence of the Euler approximation for diffusion processes. *Math. Nachr.*, 151:233–239, 1991.

[102] R. MIKULEVIČIUS et C. ZHANG : Weak Euler approximation for Itô diffusion and jump processes. *Stoch. Anal. Appl.*, 33(3):549–571, 2015.

[103] G. N. MIL'SHTEIN : Approximate integration of stochastic differential equations. *Teor. Verojatnost. i Primenen.*, 19:583–588, 1974.

[104] G. N. MIL'SHTEIN : A method of second-order accuracy integration of stochastic differential equations. *Theory of Probability & Its Applications*, 23(2):396–401, 1979.

[105] G. N. MIL'SHTEIN : Weak approximation of solutions of systems of stochastic differential equations. *Theory of Probability & Its Applications*, 30(4):750–766, 1986.

[106] S. A. MOLČANOV : Diffusion processes, and Riemannian geometry. *Uspehi Mat. Nauk*, 30(1(181)):3–59, 1975.

[107] Syoiti NINOMIYA et Nicolas VICTOIR : Weak approximation of stochastic differential equations and application to derivative pricing. *Applied Mathematical Finance*, 15(2):107–121, 2008.

[108] D. NUALART : *Malliavin calculus and related topics*. Springer Verlag, 1995.

[109] G. PAGÈS : Multi-step Richardson-Romberg extrapolation: remarks on variance control and complexity. *Monte Carlo Methods Appl.*, 13(1):37–70, 2007.

[110] G. PAGÈS : Introduction to numerical probability for finance. *handout*, page 354.p, 2016.
[111] G. PAGÈS, H. PHAM et J. PRINTEMS : *Handbook of Numerical Methods in Finance*, chapitre Optimal quantization methods and applications to numerical problems in finance, pages 253-298. Birkhauser, 2004.

[112] G. PAGÈS et J. PRINTEMS : http://www.quantize.maths-fi.com, web site devoted to optimal quantization, 2005.

[113] M. PELLETIER : Weak convergence rates for stochastic approximation with application to multiple targets and simulated annealing. *Ann. Appl. Probab.*, 8(1):10–44, 02 1998.

[114] R. PEMANTLE : Non convergence to unstable points in urn models and stochastic approximations. *Ann. Prob.*, 18(2):698–712, 1990.

[115] P. PROTTER et D. TALAY : The Euler scheme for lévy driven stochastic differential equations. *Ann. Probab.*, 1(25):393–423, 1997.

[116] S. T. RACHEV et L. RÜSCHENDORF : *Mass transportation problems. Vol. II.* Probability and its Applications (New York). Springer-Verlag, New York, 1998. Applications.

[117] H. ROBBINS et S. MONRO : A stochastic approximation method. *Ann. Math. Statistics*, 22:400–407, 1951.

[118] C. SOIZE : *The Fokker-Planck equation for stochastic dynamical systems and its explicit steady state solutions*, volume 17 de *Series on Advances in Mathematics for Applied Sciences*. World Scientific Publishing Co., Inc., River Edge, NJ, 1994.

[119] D. W. STROOCK et S. R. S. VARADHAN : Diffusion processes with boundary conditions. *Comm. Pure Appl. Math.*, 24:147–225, 1971.

[120] D. W. STROOCK et S. R. S. VARADHAN : *Multidimensional diffusion processes*, volume 233 de *Grundlehren der Mathematischen Wissenschaften [Fundamental Principles of Mathematical Sciences]*. Springer-Verlag, Berlin-New York, 1979.

[121] D. TALAY : Discrétisation d'une équation différentielle stochastique et calcul approché d'espérances de fonctionnelles de la solution. *RAIRO Modél. Math. Anal. Numér.*, 20(1):141–179, 1986.

[122] D. TALAY et L. TUBARO : Expansion of the global error for numerical schemes solving stochastic differential equations. *Stoch. Anal. and App.*, 8-4:94-120, 1990.

[123] D. TALAY et Z. ZHENG : Approximation of quantiles of components of diffusion processes. *Stoch. Process. Appl.*, 109-1:23-46, 2004.

[124] TALAY, D. : Stochastic Hamiltonian systems: exponential convergence to the invariant measure, and discretization by the implicit Euler scheme. *Markov Process. Related Fields*, 8(2):163-198, 2002.

[125] A. Yu. VERETENNIKOV : On the criteria for existence of a strong solution to a stochastic equation. *Theory of Probability & Its Applications*, 27:441–449, 1982.

[126] C. VILLANI : *Optimal transport: Old and New*. Springer, 2008.

[127] S. WATANABE : Analysis of Wiener functionals (Malliavin calculus) and its applications to heat kernels. *Ann. Probab.*, 15(1):1–39, 1987.

[128] A. K. ZVONKIN : A transformation of the state space of a diffusion process that removes the drift. *Math. USSR (Sbornik)*, 22:129–149, 1974.