<figure><img src="image_1.png" /></figure>

HAL
open science

# Maximum likelihood estimator for skew Brownian motion: the convergence rate

Antoine Lejay, Sara Mazzonetto

<figure><img src="image_2.png" /></figure>

## To cite this version:

Antoine Lejay, Sara Mazzonetto. Maximum likelihood estimator for skew Brownian motion: the convergence rate. Scandinavian Journal of Statistics, 2023, 51 (2), pp.612-642. 10.1111/sjos.12694. hal-03975966

HAL Id: hal-03975966

https://hal.science/hal-03975966

Submitted on 6 Feb 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Maximum likelihood estimator
for skew Brownian motion:
the convergence rate

Antoine Lejay* Sara Mazzonetto†

February 6, 2023

## Abstract

We give a thorough description of the asymptotic property of the maximum likelihood estimator (MLE) of the skewness parameter of a Skew Brownian Motion (SBM). Thanks to recent results on the Central Limit Theorem of the rate of convergence of estimators for the SBM, we prove a conjecture left open that the MLE has asymptotically a mixed normal distribution involving the local time with a rate of convergence of order $1/4$. We also give a series expansion of the MLE and study the asymptotic behavior of the score and its derivatives, as well as their variation with the skewness parameter. In particular, we exhibit a specific behavior when the SBM is actually a Brownian motion, and quantify the explosion of the coefficients of the expansion when the skewness parameter is close to $-1$ or $1$.

**MSC(2020) Classification:** Primary 62F12; Secondary 62F03.

**Keywords:** Skew Brownian motion; maximum likelihood estimator (MLE); statistical estimation; null recurrent process.

Skew and other singular diffusions attract more and more interest in modeling diffusive stochastic behavior in presence of semi-permeable barriers, discontinuities, and thresholds. Beyond theoretical studies, simulation and inference are also

<footnote>*Université de Lorraine, CNRS, IECL, Inria, F-54000 Nancy, France; ORCID: 0000-0003-0406-9550; Antoine.Lejay@univ-lorraine.fr</footnote>

<footnote>†Université de Lorraine, CNRS, IECL, Inria, F-54000 Nancy, France; ORCID: 0000-0001-6187-2716; Sara.Mazzonetto@univ-lorraine.fr</footnote>
necessary tools for practical purposes. For some example of applications in various fields, see e.g. [28, 9, 31, 36, 30, 8, 10, 7, 13] among others.

The inference of skew diffusion cannot follow from a simple adaptation of known techniques for Stochastic Differential Equations (SDE) as the ones presented in [12, 15]. In fact their distributions are singular with the ones of classical SDE. The limits are usually mixed normal ones and the rate is not necessarily 1/2. The work dealing with the inference of skew diffusion is rather limited. Let us cite however [6, 21, 22, 27, 17, 26, 25, 24, 33, 34] for frequentist inference and [5, 3, 4] for Bayesian inference.

The Skew Brownian motion (SBM) is a basic brick for constructing Skew diffusion, as several transformations reduces Skew diffusions to SBM [16]. This latter process depends on a single parameter $\theta \in [-1, 1]$ — the skewness parameter — which rules out its behavior when it crosses zero. For $\theta = 0$, the SBM is a Brownian motion. For $\theta = \pm 1$, it is a Reflected Brownian motion.

A series of works considers the inference of the skewness parameter from high-frequency observations. In [21], the authors have given an asymptotic expansion of the Maximum Likelihood Estimator (MLE) $\theta_n$ of $\theta$ around $\theta = 0$ in power of $n^{-1/4}$, where $n$ is the number of observations. A heuristic explanation of the power 1/4 is given by analogy with the Skew Random Walk [18], where the MLE depends on the local time at zero of the discrete walk, which a random quantity of order $\sqrt{n}$. Indeed in [22], where the consistency of the MLE and another estimator is proved, it was empirically observed that the rate of convergence should be 1/4, meaning that the observed points that “carry the information” are those close to 0, and are of order $\sqrt{n}$. This also explains why the asymptotic limit of the MLE involves the local time at 0. The limit is of type $\mathcal{s}(\theta)G/\sqrt{L_1}$ (when the process is observed on $[0, 1]$), where $(L_t)_{t \ge 0}$ is the symmetric local time at 0 of the SBM and $G$ is a centered unit, Gaussian independent from the SBM. The value of $\mathcal{s}(\theta)$ was empirically observed as closed to $\kappa\sqrt{1 - \theta^2}$. It was conjectured in [22] that $\mathcal{s}(\theta) = \kappa\sqrt{1 - \theta^2}$ as for the Skew Random Walk.

The article [27] brought the missing result required to establish a Central Limit Theorem on the MLE and other estimators. A related result may also be found in [32]. As for the results in [22], this is based on an extension of the work of J. Jacod [14] that cannot be applied directly as the SBM has a singular distribution with respect to the one of the BM.

In this article, we first prove the asymptotic normality of the MLE of the skewness parameter $\theta$ of the SBM. More precisely, we establish that

$$n^{1/4}(\theta_n - \theta) \xrightarrow{n \to \infty} \mathcal{s}(\theta)W(L_T) \stackrel{\text{law}}{=} \mathcal{s}(\theta) \frac{W(1)}{\sqrt{L_T}}$$
for a Brownian motion W independent from the SBM. In particular, we show that $s(\theta) = \kappa(\theta)\sqrt{1 - \theta^2}$ with a pre-factor $\kappa(\theta)$ that varies slowly.

Second, using a recent asymptotic inversion formula [20], we establish a series expansion of $\theta_n$ as

$$ \theta_n = \theta + \sum_{k=1}^{+\infty} D_{k,n}(\theta) \left( \frac{s(\theta)P_n(\theta)}{T^{1/4}\eta^{1/4}} \right)^k, \quad (1) $$

for some random quantities $D_{k,n}(\theta)$ whose asymptotic behavior is also studied, and $P_n(\theta)$ is asymptotically pivotal (i.e., the asymptotic law does not depend on the parameter). All these results are based on the asymptotic behavior of the score and their derivatives at any orders. Expansion (1) provides some insight on the behavior of the MLE in function of the number of samples $n$ and the true value of $\theta$. In particular, the closer $\theta$ is to $\pm 1$, the more skewness is observed.

Furthermore, we show that $D_{2k+1,n}(\theta = 0)$ converges to 0 so that a more precise expression of the involving a multivariate mixed Gaussian distribution could be given for the MLE of $\theta = 0$. This expression is an alternative to the one already given in [21].

Moreover, we give some numerical experiments on the rate of convergence of the distribution functions of the MLE and the score and its derivatives towards their limiting distributions. The rate of convergence of this Berry-Esseen type analysis seems to depend on $\theta$. This will be subject to further study.

Finally, we study the behavior of the limiting coefficients in function of $\theta$. In particular, the expansion in (1) exhibits a boundary layer estimate for $\theta$ close to $\pm 1$ as the coefficients explode in powers of $(1 - \theta^2)^{-1/2}$.

**Outline.** We present our main results in Section 1. They are built on the asymptotic behavior of the score and its derivative which we give in Section 2. We also study numerically the rate of convergence of the scores and their derivatives towards their limits in Section 2.1. The properties of the limiting coefficients are studied in Section 3. And the proofs of our main theorems are given in Section 4.

# 1 Main results

The SBM $X$ of parameter $\theta \in [-1, 1]$ solves the SDE

$$ X_t = B_t + \theta L_t, \quad t \ge 0 \quad (2) $$

for $B$ a Brownian motion and $\{L_t\}_{t \ge 0}$ its symmetric local time of $X$ at point 0. Actually, the SDE (2) has a unique strong solution [16, 11]. No solution exists
when $|\theta| > 1$. For $\theta = \pm 1$, the SBM is a (positively if $\theta = 1$ or negatively if $\theta = -1$) Reflected Brownian motion.

We denote by $(\Omega, \mathcal{F}, \mathbb{P}_\theta)$ the underlying probability space and by $\{\mathcal{F}_t\}_{t \ge 0}$ the filtration with respect to which $X$ is adapted. This filtration may be taken as the one generated by the driving Brownian motion and may be assumed to satisfy the usual conditions (i.e. it is complete and right continuous).

We are concerned with the estimation of the parameter of the SBM observed at discrete times over a finite time window. We will establish limits in high-frequency.

**Data 1.** We observed the SBM $\{X_t\}_{t \in [0,T]}$ of parameter $\theta_0 \in (-1, 1)$ at times $\{t_i\}_{i=0,...,|nT|}$ on a time window $[0, T]$ with $t_i := i/n$. The starting point is $X_0 = 0$. Note that $n$ is the number of sample per unit of time.

*Remark 1.* We impose $X_0 = 0$ to ensure that $L_T \neq 0$. When $X_0 \neq 0$, we could still apply our methodology on a random window $[\tau_0 \wedge T, T]$ where $\tau_0$ is the first hitting time from 0. If $\tau_0 > T$, then no observation can be used to estimate $\theta$.

The density transition function of the SBM of parameter $\theta$ is $[35, 16]$

$$p_{\theta}(t, x, y) := \frac{1}{\sqrt{2\pi t}} \begin{cases} \exp\left(-\frac{(x-y)^2}{2t}\right) + \theta \exp\left(-\frac{(x+y)^2}{2t}\right) & \text{if } x \ge 0, y \ge 0, \\ \exp\left(-\frac{(x-y)^2}{2t}\right) - \theta \exp\left(-\frac{(x+y)^2}{2t}\right) & \text{if } x \le 0, y \le 0, \\ (1-\theta) \exp\left(-\frac{(x-y)^2}{2t}\right) & \text{if } x \ge 0, y \le 0, \\ (1+\theta) \exp\left(-\frac{(x-y)^2}{2t}\right) & \text{if } x \le 0, y \ge 0. \end{cases}$$

The SBM is null recurrent process with invariant measure $\mu_\theta(dx) := \mu_\theta(x) dx$ with

$$\mu_{\theta}(x) := \begin{cases} 1 + \theta & \text{if } x \ge 0, \\ 1 - \theta & \text{if } x < 0. \end{cases}$$

When observed at regular times as in Data 1, we call *likelihood* the random function:

$$\Lambda_n(\theta) = \prod_{i=0}^{\lfloor nT \rfloor - 1} p_{\theta}(\Delta t, X_{t_i}, X_{t_{i+1}}), \quad \theta \in [-1, 1] \text{ with } \Delta t = \frac{1}{n}.$$

Since $[-1, 1] \ni \theta \mapsto p_{\theta}(\Delta t, x, y)$ is analytic, $\theta \mapsto \Lambda_n(\theta)$ is also analytic. The *score* is $\partial_{\theta} \log \Lambda_n(\theta)$.

Let us define <fnref n="1" />

$$k_{\theta}(x, y) := \partial_{\theta} \log p_{\theta}(1, x, y) = \frac{\partial_{\theta} p_{\theta}(1, x, y)}{p_{\theta}(1, x, y)} = \frac{\operatorname{sgn}(y)}{\operatorname{sgn}(y)\theta + \exp(2(xy)^+)}$$

<footnote n="1">¹In the definition of $k_\theta(x, y)$, there is a mistake in [22]: $\Delta t$ has to be replaced by $T$. Note that in [22], $n$ is the number of samples while here, $n$ is the number of samples per unit. Therefore, the limit distribution in [22] is the local time at time 1 of $\sqrt{T}X$, while here it is the one of $X$.</footnote>
<figure><img src="image_4.png" /><figcaption>Figure 1: Density of the MLE for various values of $\theta$ using 10 000 samples of the SBM using $\Delta t = 10^{-4}$ for $T = 1$.</figcaption></figure>

where $(xy)^+$ stands for $\max(0, xy)$. For any $n > 0$, the following scaling holds true:

$$k_{\theta}(x\sqrt{n}, y\sqrt{n}) = \partial_{\theta} \log p_{\theta}(\Delta t, x, y) \text{ with } \Delta t = \frac{1}{n}.$$

*Remark 2.* The score rewrites $\partial_{\theta} \log \Lambda_n(\theta) = \sum_{i=0}^{[nT]^{-1}} k_{\theta}(X_{t_i}\sqrt{n}, X_{t_{i+1}}\sqrt{n})$.

**Proposition 1 ([22]).** The maximum likelihood estimator (MLE)

$$\theta_n := \arg\max_{\theta \in [-1,1]} \Lambda_n(\theta)$$

is the unique solution to $\partial_{\theta} \log \Lambda_n(\theta_n) = 0$ and is a consistent estimator of the parameter $\theta_0$ of the SBM under $P_{\theta_0}$.

The goal of this paper is to refine the latter result providing asymptotic information.

In Figure 1, we plot the empirical density of the MLE for various values of $\theta$. We use the method in [23] for the simulation of the SBM, while the MLE associated to each trajectory is obtained maximizing numerically the log-likelihood. We see that the MLE is concentrated around the true value. The more $\theta$ is closer to 1, the more the density is skewed to the left. In the next section we examine this behavior. For reasons related to symmetry, for the remainder of the document, we restrict to consider the case $\theta \ge 0$.

## 1.1 Asymptotic mixed normality of the MLE estimator

We establish, in the forthcoming Theorem 1, the asymptotic mixed normality of the MLE, with the non standard rate of $n^{1/4}$. This result has already been proven,
when $\theta = 0$, in [21]. Before to state the result, let us first introduce a definition specifying the nature of the asymptotic normality.

**Definition 1** (Class of $\mathcal{L}$-mixed normal distribution). Let $L_1$ be the local time at point 0 and time 1 of the SBM $X$ and $G \sim \mathcal{N}(0, 1)$ be independent from $\mathcal{F}_1$ (hence from $L_1$). Note that the distribution of $L_1$ does not depend on $\theta$. A random variable $M$ is said to be $\mathcal{L}$-mixed normal distributed if

$$M \stackrel{\text{law}}{=} \frac{G}{\sqrt{L_1}}.$$

Note that a $\mathcal{L}$-mixed normal distribution has an infinite second moment. It is a symmetric, unimodal distribution with heavy tails.

*Remark 3* (Simulation of the local time). The local time of the SBM is equal in distribution to the one of the Brownian motion, hence $L_1 \stackrel{\text{law}}{=} |H|$ for $H \sim \mathcal{N}(0, 1)$. Therefore, the local time and $\mathcal{L}$-mixed normal distributions are easily simulated.

*Remark 4* (Scaling). Let $L$ be the symmetric local time of a SBM $X$ at the point 0 and $W$ be a Brownian motion $W$ independent from $X$. They both satisfy a scaling property, in particular for the local time it holds $L_T \stackrel{\text{law}}{=} \sqrt{T}L_1$. Hence,

$$\frac{W(L_T)}{L_T} \stackrel{\text{law}}{=} \frac{W(1)}{\sqrt{L_T}} \stackrel{\text{law}}{=} \frac{W(1)}{T^{1/4}\sqrt{L_1}}.$$

In other words, $T^{1/4}W(L_T)/L_T$ and $T^{1/4}G/\sqrt{L_1}$ (for $G \sim \mathcal{N}(0, 1)$ independent from $X$) are $\mathcal{L}$-mixed normal distributed.

We also introduce a quantity related to the asymptotic variance of the estimator and to the Fisher information:

$$s(\theta) := \left( - \iint_{\mathbb{R}^2} \mu_{\theta}(x) p_{\theta}(1, x, y) \partial_{\theta} k(x, y) \, dx \, dy \right)^{-1/2} \in \mathbb{R}. \quad (3)$$

The quantity $s(\theta)$ will be studied in more details in Proposition 2, Remark 5, and in Section 3, in particular in Section 3.3.

We are now ready to state the main result which says that $(nT)^{1/4}s(\theta)^{-1}(\theta_n - \theta)$ is asymptotically a $\mathcal{L}$-mixed normal distribution under $\mathbb{P}_{\theta}$.

**Theorem 1** (Asymptotic $\mathcal{L}$-mixed normality of the MLE for the SBM). Let $\theta \in (-1, 1)$. The MLE estimator $\theta_n$ is asymptotically mixed normal with

$$n^{1/4}(\theta_n - \theta) \xrightarrow[n \to \infty]{\text{stably}} s(\theta) \frac{H}{\sqrt{L_T}}, \text{ under } \mathbb{P}_{\theta},$$

where $H \sim \mathcal{N}(0, 1)$ independent from $\mathcal{F}_T$ (hence of $L_T$).
The proof of Theorem 1 follows from the results of the next section, Section 1.2, which rely on the study of the score and its derivatives that we propose in Section 2.

Heuristically, the non standard rate of $n^{1/4}$, can be explained by the fact that the quality of the estimation depends mainly on the time spent by the SBM around 0, and that the fraction of observations when $\{X_{t_i}\}_{i=0,...,n}$ is of order $\sqrt{n}$. This fact is rigorously established for the Skew Random Walk where the local time is really the occupation time at the point where the bias-dynamic is perturbed [19].

Although it was conjectured in [22] from the results on the Skew Random Walk that the coefficient $\delta(\theta)$ in front of the mixed Gaussian should be proportional to $\sqrt{1 - \theta^2}$, we find a slowly varying pre-factor.

**Proposition 2.** For all $\theta \in (-1, 1)$, the function $(x, y) \mapsto \mu_{\theta}(x)p_{\theta}(1, x, y)\partial_{\theta}k(x, y)$ is integrable and its integral is negative, so $\delta(\theta) \in (0, \infty)$ is well defined. Besides there exist two real constants $0 < c_1 \le c_2 < \infty$ such that for all $\theta \in (-1, 1)$

$$c_{-}\sqrt{1 - \theta^{2}} \le \delta(\theta) \le c_{+}\sqrt{1 - \theta^{2}}.$$

We find that $0.79 \le c_{-} < c_{+} \le 0.88$, which is consistent with the numerical observations of [22].

The proof of Proposition 2 is provided in Section 3.3 where a more precise statement is formulated. Actually we show in Remark 14 that an accurate approximation of $\delta(\theta)$ is given by

$$\delta(\theta) \approx \frac{\sqrt{1 - \theta^2}}{\sqrt{1.292 + 0.232\theta^2 + 0.071\theta^4}}.$$

## 1.2 Asymptotic expansion for the MLE estimator

Let us first consider the following family of statistics of interest:

$$S_m(n, \theta) := \frac{1}{n^{1/2}} \sum_{i=0}^{[nT] - 1} \partial_{\theta}^m k_{\theta}(X_{t_i} \sqrt{n}, X_{t_{i+1}} \sqrt{n}) \text{ for } m \ge 0. \quad (4)$$

Remark 2 shows that $S_m(n, \theta) = \partial_{\theta}^m(\partial_{\theta} \log \Lambda_n(\theta))/n^{1/2}$ so that $S_m(n, /\theta)$ is the renormalized m-th order derivative of the score.

Let us also consider the statistics, on $\{S_1(n, \theta) \neq 0\}$,

$$d_{k,n}(\theta) := \frac{-1}{k!} \frac{S_k(n, \theta)}{S_1(n, \theta)} \text{ for } k \ge 0 \quad \text{and} \quad P_n(\theta) := n^{1/4} \frac{T^{1/4}}{\delta(\theta)} d_{0,n}(\theta),$$
as well as the constants

$$ \xi_m(\theta) := \iint_{\mathbb{R}^2} \mu_\theta(x)p_\theta(1, x, y) \partial_\theta^m k(x, y) \, dx \, dy, \quad \text{for } m \in \mathbb{N}. \quad (5) $$

The integrability of $(x, y) \mapsto \mu_\theta(x)p_\theta(1, x, y)\partial_\theta^m k(x, y)$ is shown in Section 3, together with other properties of $\xi_m(\theta)$.

Remark 5. The quantity $s(\theta)$ given by (3) is related to $\xi_1(\theta)$: $s(\theta)^{-2} = -\xi_1(\theta)$.

**Proposition 3.** Under $\mathbb{P}_\theta$,

(i) the statistics $\mathbb{P}_n(\theta) = n^{1/4}T^{1/4}\mathbf{d}_{0,n}(\theta)/s(\theta)$ is asymptotically $\mathcal{L}$-mixed normal distributed;

(ii) for every $k \ge 1$,

$$ \mathbf{d}_{k,n}(\theta) \xrightarrow{n \to \infty} \mathbf{d}_k(\theta) := \frac{-\xi_k(\theta)}{k!\xi_1(\theta)} = \frac{s(\theta)^2}{k!}\xi_k(\theta). \quad (6) $$

Moreover, under $\mathbb{P}_0$, for every $k \ge 1$, $\mathbf{d}_{2k,n}(0) \xrightarrow{n \to \infty} 0$ and $\mathbf{d}_{2k}(0)$ is, up to a multiplicative constant, $\mathcal{L}$-mixed normal distributed. Furthermore $\{n^{1/4}\mathbf{d}_{2k,n}(0)\}_{k=0,...,m}$ converges stably for any $m \ge 1$. The limit is identified in Proposition 4 in Section 2.

*Proof.* This is an immediate consequence of Proposition 4 in Section 2 on the asymptotic behavior of the score combined with Remark 5. $\square$

In Theorem 2, we give an expansion of the MLE in term of $\mathbf{d}_{0,n}(\theta)$ (and so of $s(\theta)\mathbb{P}_n(\theta)/T^{1/4n^{1/4})}$. It follows from applying Theorem 3 in [20]. For $\theta = 0$, such a type of expansion was already given in [21] in the form provided in equation (8) below. We also provide an alternative expansion based on a finer analysis of the coefficients' asymptotic behavior.

Let us introduce some notation: for any $m \in \mathbb{N} \cup \{\infty\}$, let

$$ \Phi^{[m]}(\delta, x) := \sum_{k=1}^{m} \delta_k x^k \quad (7) $$

be the formal power series in $x$ of coefficients $\delta = \{\delta_k\}_{k \ge 1}$.

**Theorem 2** (Asymptotic expansion, cf. [20, Theorem 3]). For any integer $n \ge 1$, any $N \ge 1$, and any $\theta \in (-1, 1)$, the MLE satisfies under $\mathbb{P}_\theta$,

$$ \theta_n = \theta + \Phi^{[\infty]}(\mathbf{D}_{\cdot,n}(\theta), \mathbf{d}_{0,n}(\theta)) = \theta + \sum_{k=1}^{+\infty} \mathbf{D}_{k,n}(\theta)\mathbf{d}_{0,n}(\theta)^k, \quad (8) $$
where $D_{1,n}(\theta) := 1$ and $D_{q,n}(\theta)$, $q = 2, \dots, N$, are given by the following recursive formula:

$$D_{q,n}(\theta) := \sum_{m=2}^{q} d_{m,n}(\theta) \sum_{k_1+\dots+k_m=q} D_{k_1,n}(\theta) \dots D_{k_m,n}(\theta). \quad (9)$$

If $\theta = 0$, another expansion holds

$$\theta_n = \Phi^{[\infty]}(a_{\cdot,n}, n^{-1/4}) = d_{0,n}(0) + \sum_{\substack{k \ge 3 \\ k \text{ odd}}} \frac{a_{k,n}}{n^{k/4}}, \quad (10)$$

where $a_{k,n}$ for $k = 0, \dots, N$ is $a_{2q,n} := 0$ if $q \ge 0$, and

$$a_{1,n} := n^{1/4}d_{0,n}(0) = \delta(0)P_n(0)/T^{1/4},$$

$$a_{2q+1,n} := \sum_{\substack{m=2 \\ m \text{ even}}}^{2q} n^{1/4}d_{m,n}(0) \sum_{k_1+\dots+k_m=2q} a_{k_1,n} \dots a_{k_m,n} \\ + \sum_{\substack{m=3 \\ m \text{ odd}}}^{2q+1} d_{m,n}(0) \sum_{k_1+\dots+k_m=2q+1} a_{k_1,n} \dots a_{k_m,n}.$$

*Remark 6.* From [20], it can be seen that $\Phi(d_{\cdot,n}(\theta), \Phi(D_{\cdot,n}(\theta), x)) = -x$. For the first coefficients, it holds

$$D_{2,n} = d_{2,n}, \quad D_{3,n} = d_{3,n} + 2d_{2,n}^2,$$

$$D_{4,n} = d_{4,n} + 5d_{2,n}d_{3,n} + 5d_{2,n}^3$$

$$\text{and } D_{5,n} = d_{5,n} + 6d_{2,n}d_{4,n} + 21d_{2,n}^2d_{3,n} + 14d_{2,n}^4 + 3d_{3,n}^2.$$

In the case of the sequence in (10),

$$a_{3,n} = d_{3,n}(0)P_n(0)^3\delta(0)^3/T^{3/4} + n^{1/4}d_{2,n}(0)P_n(0)^2\delta(0)^2/T^{1/2}$$

$$\text{and } a_{5,n} = (d_{5,n}(0) + 3(d_{3,n}(0))^2)P_n(0)^5\delta(0)^5/T^{5/4}$$

$$+ (5n^{1/4}d_{2,n}(0)d_{3,n}(0) + n^{1/4}d_{4,n}(0))P_n(0)^4\delta(0)^4/T$$

$$+ 2n^{1/2}(d_{2,n}(0))^2P_n(0)^3\delta(0)^3/T^{3/4}.$$

*Remark 7.* As $d_{k,n}$ and $P_n$ depend on $\theta$, they cannot be computed under the true parameter in the context of estimation (actually, $P_n(\theta_n) = 0 = d_{0,n}(\theta_n)$ for the MLE $\theta_n$) but can be used for statistical hypothesis testing.

The result shows a sort of “phase transition” between $\theta = 0$ and $\theta \neq 0$. This is clearly related to the dichotomy in the convergence of $d_{2m,n}(0)$ (vanishing) and $n^{1/4}d_{2m,n}(0)$
for $m \ge 1$: Observe, for instance, that the second term in equation (10) is of order $1/n^{3/4}$ and in equation (8) is of order $1/n^{1/2}$. This latter term, in the case $\theta = 0$ goes to 0 as $n$ increases.

Theorem 2 and Proposition 3 prove Theorem 1. Moreover they imply, for all $m \in \mathbb{N}$, approximations of the MLE given by the formal power series

$$ \theta_n^{[m]} := \theta + \Phi^{[m]}(D_{\cdot,n}(\theta), d_{0,n}(\theta)) \quad \text{and} \quad \tilde{\theta}_n^{[m]} := \theta + \Phi^{[m]}(D(\theta), d_{0,n}(\theta)), \quad (11) $$

where $D_k(\theta)$ is defined similarly to $D_{k,n}(\theta)$ in (9) with $d_{k,n}$ replaced by $d_n$ in (6). More precisely, $D_1(\theta) := 1$ and for $q = 2, \dots, N$, $D_q$ is given by the recursive formula:

$$ D_q(\theta) := \sum_{m=2}^{q} d_m(\theta) \sum_{k_1+\dots+k_m=q} D_{k_1}(\theta) \dots D_{k_m}(\theta). \quad (12) $$

Then, for instance,

$$ \tilde{\theta}_n^{[3]} = \theta + d_{0,n} - \frac{1}{2} \frac{\xi_2(\theta)}{\xi_1(\theta)} (d_{0,n})^2 + \frac{1}{2} \left( \left( \frac{\xi_2(\theta)}{\xi_1(\theta)} \right)^2 - \frac{\xi_3(\theta)}{3\xi_1(\theta)} \right) (d_{0,n})^3. $$

The power series $\theta_n^{[m]}$ is the $m$-th order truncation of (8). It has random coefficients $D_{\cdot,n}(\theta)$ and random argument $d_{0,n}(\theta)$ involving the score and its derivatives. The proxy $\tilde{\theta}_n^{[m]}$ is a power series with deterministic coefficients $D(\theta)$ (limit of $D_{\cdot,n}(\theta)$) and the same random argument $d_{0,n}(\theta)$ related to the score and its first derivative. Both $\theta_n^{[m]}$ and $\tilde{\theta}_n^{[m]}$ are non-linear expressions of $d_{0,n}(\theta)$.

The following result enlightens the coefficients' behavior as $\theta$ varies in $(-1, 1)$. It is proved in Section 3.

**Lemma 1.** Let $\{\xi_m\}_{m \in \mathbb{N}}$ defined in (5). Then for all $m \in \mathbb{N}$:

(i) $\xi_0 = 0$;

(ii) $\theta \mapsto \xi_{2m+1}(\theta)$ is even;

(iii) There exists $-\infty < c_{m,1} < c_{m,2} < 0$ such that $-c_{m,1} \le (1-\theta^2)^{2m+1}\xi_{2m+1}(\theta) \le c_{m,2}$ for all $\theta \in (-1, 1)$;

(iv) $\xi_{2m}(\theta)$ is odd, in particular $\xi_{2m}(0) = 0$;

(v) There exists $c_{2m} \in (0, \infty)$ such that $-c_{2m} \le (1-\theta^2)^{2m}\xi_{2m}(\theta) < 0$ for all $\theta \in (0, 1)$.

The following result is proved by induction on the recursive formula (9).

**Corollary 1** (of Proposition 3 and Lemma 1). It holds that

$$ D_{k,n}(\theta) \xrightarrow{n \to \infty} D_k(\theta). $$
<figure><img src="image_5.png" /><figcaption>Figure 2: Functions $\theta \mapsto D_k(\theta)\delta(\theta)^k(1 - \theta^2)^{k/2-1}$ with $D_k$ defined by (12) for $\theta$ restricted to $[0, 1]$ and $m = 1, \dots, 5$.</figcaption></figure>

Besides, for any $k \ge 1$,

$$|D_k(\theta)\delta(\theta)^k| \le \frac{C_k}{(1 - \theta^2)^{k/2-1}}$$

for a constant $C_k \ge 0$.

We plot the coefficients $\theta \mapsto D_k(\theta)\delta(\theta)^k$ in Figure 2, up to this scaling factor $\theta \mapsto (1 - \theta^2)^{1-k/2}$.

*Remark 8.* The coefficients $D_{k,n}(\theta)$ converge towards their limits at rate of $n^{1/4}$. This is a consequence of definition (9) and Proposition 3.

In Figure 3 we show the empirical densities of the MLE $\theta_n$, $\theta_n^{[1]} = \theta + d_{0,n}(\theta)$, $\theta_n^{[5]}$, and of the proxy $\tilde{\theta}_n^{[5]}$ (see (11) for the definitions of $\theta_n^{[5]}$ and $\tilde{\theta}_n^{[5]}$). For $\theta = 0$, $d_{0,n}(0)$ replicates already very well the MLE behavior, therefore there $\theta_n^{[5]}$ and $\tilde{\theta}_n^{[5]}$ are not plotted. For $\theta \ne 0$, we observe that, while the lower order expansions replicate worse the MLE behavior, $\theta_n^{[5]}$ and $\tilde{\theta}_n^{[5]}$ do it quite well, and so do the higher order expansions. The expansions $\theta_n^{[5]}$ and $\tilde{\theta}_n^{[5]}$ are close one to another. Indeed they are respectively random and deterministic polynomials of $d_{0,n}(\theta)$ such that, by Corollary 1, the random coefficients converge to the deterministic ones. To introduce the next paragraph, observe that the density of $\theta_n^{[1]} = \theta + d_{0,n}(\theta)$ is skewed to the left when $\theta$ is close to 1. Of course, since their density is close to the one of the MLE's, $\theta_n^{[5]}$ and $\tilde{\theta}_n^{[5]}$ exhibit a skewed empirical distribution.

*Remark 9.* Let us remind that as for $\theta_n - \theta$, $n^{1/4}d_{0,n}(\theta)$ has asymptotic symmetric distribution: Proposition 3 ensures that $n^{1/4}d_{0,n}(\theta) = P_n(\theta)\delta(\theta)T^{-1/4}$ is asymptotically $\mathcal{L}$-mixed normal distributed. Therefore the alternative proxys
<figure><img src="image_6.png" /><figcaption>Figure 3: Density of MLE, $\theta_n^{[1]} = \theta + d_{0,n}(\theta)$, and the fifth order expansions $\theta_n^{[5]}$ and $\tilde{\theta}_n^{[5]}$ (11), for several values of $\theta$ using 10 000 samples of the SBM using $\Delta t = 10^{-3}$ for $T = 1$ (thus, $n = 1000$).</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 4: Functions $\theta \mapsto f_m(\theta)$ defined by (13) for $\theta$ restricted to $[0, 1]$ and $m = 1, \dots, 5$.</figcaption></figure>

$\Phi^{[m]}(D(\theta), P_{\infty}(\theta)\delta(\theta)(nT)^{-1/4})$ where $P_{\infty}$ follows a $\mathcal{L}$-mixed normal distribution, would not replicate skewness and they would need higher order expansions than the proxys $\tilde{\theta}_n^{[m]}$ in order to have empirical distribution which gets close to the MLE's one.

**Skewness.** We have observed skewness for the empirical distribution $d_{0,n}(\theta)$ and of the MLE. This is due to the fact that the empirical density of the score shows that the more $\theta$ is close to 1 the more the score is skewed to the left. However, the correlation between MLE and score is non-trivial: non-linear dependence. The approximations of the MLE in (11) are non-linear functions of the score and its derivative, actually of $d_{0,n}$ which is basically their ratio. Since $D_2(\theta) = d_2(\theta) = -\frac{1}{2}\frac{\xi_2(\theta)}{\xi_1(\theta)}$, Lemma 1 establishes that $\text{sgn}(D_2(\theta)) = -\text{sgn}(\theta)$. This favors a skewness to the left (resp. right) when $\theta > 0$ (resp. $\theta < 0$) of the distribution of the proxys (11). In other words, the approximations of the MLE proposed above capture the skewness of the MLE.

**Boundary layer effect.** Let us set

$$f_m(\theta) := \frac{(-1)^m (1 - \theta^2)^m}{m!} \xi_m(\theta). \qquad (13)$$

Then $f_{2k+1}$ is positive and bounded (from above and below by a positive constant),
$f_{2k}(0) = 0$ and $f_{2k}(\theta)$ is bounded. Then, (11) rewrites:

$$ \tilde{\theta}_n^{[3]} = \theta + \frac{\sqrt{1-\theta^2}}{\sqrt{f_1(\theta)}} \frac{P_n(\theta)}{(Tn)^{1/4}} + \frac{f_2(\theta)}{f_1(\theta)^2} \left( \frac{P_n(\theta)}{(Tn)^{1/4}} \right)^2 + \frac{1}{\sqrt{1-\theta^2}} \frac{2(f_2(\theta))^2 - f_3(\theta)f_1(\theta)}{f_1(\theta)^{7/2}} \left( \frac{P_n(\theta)}{(Tn)^{1/4}} \right)^3, $$

that is

$$ D_1(\theta)\delta(\theta) = \frac{\sqrt{1-\theta^2}}{\sqrt{f_1(\theta)}}, \quad D_2(\theta)\delta(\theta)^2 = \frac{f_2(\theta)}{f_1(\theta)^2} $$

$$ \text{and } D_3(\theta)\delta(\theta)^3 = \frac{1}{\sqrt{1-\theta^2}} \frac{2f_2(\theta)^2 - f_3(\theta)f_1(\theta)}{f_1(\theta)^{7/2}}. $$

The term of order 2 vanishes for $\theta = 0$ as $f_2(0) = 0$. With this expansion, one sees a "boundary layer" effect for $|\theta|$ close to 1 in the explosion of the third order coefficient. Indeed if we push the expansion up to order $k$ (here we stopped at $k = 3$), the corresponding coefficient explodes as $(1 - \theta^2)^{(k-2)/2}$ when $|\theta|$ is close to 1. We have performed numerical simulations that suggest that no cancellation effect occurs so that the approximation by a polynomial expansion is no longer suitable.

In [20] we discussed a similar boundary layer phenomenon in the case of the Binomial family. One could expect the same to happen for the skewness parameter $\theta$ because of the pathwise construction of SBM done associating independent Bernoulli random variables with parameter $(1 + \theta)/2$ to the excursions from 0 of a reflected Brownian motion and flipping each excursion based on the result of the Bernoulli random variable.

**Change of variable and change of coordinates.** Combining the Faà di Bruno formula with (8), one may give some explicit expansion of $\varphi(\theta_n)$ for any analytic function $\varphi$. Similarly, one may also consider $\theta_n$ in another system of coordinates, as discussed in [20]. On that point, two changes of variables appear to be natural: $\varphi_1(\theta) = \sqrt{1 - \theta^2}$ and $\varphi_1(\theta) = 1/\sqrt{1 - \theta^2}$. The latter one stabilizes the variance. However, no change of coordinate impacts the asymptotic behavior of the score. Besides, we found through numerical experiments that using $\varphi_1$ or $\varphi_2$ does not improve Wald confidence intervals.
## 2 Asymptotic behavior of the score and its derivatives

In this section, we provide results which are necessary to prove Proposition 3. These results are an application of the results in [14] for the Brownian motion, and on the ones from [22] (convergence) and [27] (Central Limit Theorem) for the (true) SBM.

Let us recall that $S_m$ in (4) is nothing else that a rescaled derivative of the score and $\xi_m(\theta)$ is a constant defined in (5). By Lemma 1, it holds that for all $\theta \in (-1, 1)$, $\xi_0(\theta) = 0$ and $\xi_1(\theta)$ is negative, and that for all integer $m$, $\xi_{2m}(0) = 0$.

**Proposition 4.** For any $\theta \in (-1, 1)$, the family $\{\xi_k(\theta)\}_{k \ge 0}$ defined by (5) is such that:

(i) Under $\mathbb{P}_\theta$, there exists a standard (with mean 0 and variance 1) Gaussian random variable $H$ independent from $\mathcal{F}_T$ (the probability space $(\Omega, \mathcal{F}, \mathbb{P}_0)$ has been extended to carry $H$) such that for any $m \ge 0$,

$$ n^{1/4}S_0(n, \theta) \xrightarrow[n \to \infty]{\text{stably}} \sqrt{-\xi_1(\theta)}\sqrt{L_T}H, \\ (S_1(n, \theta), \dots, S_m(n, \theta)) \xrightarrow[n \to \infty]{\text{prob.}} (\xi_1(\theta)L_T, \dots, \xi_m(\theta)L_T). $$

Thanks to the property of the stable convergence, we have joint stable convergence of $(n^{1/4}S_0(n, \theta), S_1(n, \theta), \dots, S_m(n, \theta))$.

(ii) For $\theta = 0$ (the SBM is actually a Brownian motion), under $\mathbb{P}_0$, there exists a Gaussian family $H = \{H_k\}_{k=0,\dots,m}$ independent from $\mathcal{F}_T$ with mean 0 and covariance $\text{Cov}(H_j, H_\ell) = \Psi_{2j,2\ell}$ for $0 \le j, \ell \le m$ described in Section 3.5 (the probability space $(\Omega, \mathcal{F}, \mathbb{P}_0)$ has been extended to carry $H$) such that for any $m \ge 0$,

$$ \{n^{1/4}S_{2k}(n, 0)\}_{k=0,\dots,m} \xrightarrow[n \to \infty]{\text{stably}} \{\sqrt{L_T}H_k\}_{k=0,\dots,m}, \\ \{S_{2k+1}(n, 0)\}_{k=0,\dots,m} \xrightarrow[n \to \infty]{\text{prob.}} \{\xi_{2k+1}(0)L_T\}_{k=0,\dots,m}. $$

Thanks to the property of the stable convergence, we have joint stable convergence of $(n^{1/4}S_0(n, \theta), S_1(n, \theta), n^{1/4}S_2(n, \theta), \dots)$.

**Remark 10.** Following [14, 27], a stronger result holds: on an extended probability space, there exist a $m$-dimensional Brownian motion $\{W_k\}_{k \ge 0}$ independent from $L$
and a $m \times m$-matrix $\Xi(\theta)$ such that

$$ \left\{ \frac{1}{n^{1/4}} \left( \frac{1}{\sqrt{n}} \sum_{i=1}^{\lfloor nt \rfloor} \partial^\ell k_\theta(X_{i-1}\sqrt{n}, X_i\sqrt{n}) - \xi_\ell(\theta)L_t \right) \right\}_{\substack{\ell=0,\dots,m, \ t \in [0,T] \\ \xrightarrow{n \to \infty} \Xi(\theta)W(L_t) \}} \quad (14) $$

A closed-form — yet cumbersome — expression exists for the matrix $\Xi(\theta) = (\Xi_{i,j}(\theta))_{i,j=0,\dots,m}$. Here, we focus on the particular cases $m = 0$ and $\theta = 0$ and we study $\{\xi_k(\theta)\}_{k=1,\dots,\max(m,1)}$ and $\Xi$ for these cases.

Remark 11. When $\xi_m(\theta) = 0$, (14) rewrites

$$ n^{1/4}S_m(n, \theta) \xrightarrow[n \to \infty]{\text{stably}} \Xi_{m,m}(\theta)\sqrt{L_T}G $$

for some constant $\Xi_{m,m}(\theta)$, and $G \sim \mathcal{N}(0, 1)$ independent of $L_T$. In particular,

(i) For any $\theta \in (-1, 1)$ and $m = 0$,

$$ \Xi_{0,0}(\theta) = \sqrt{-\xi_1(\theta)} = \frac{1}{s(\theta)}. $$

This follows from [27] and Lemma 4.

(ii) For $\theta = 0$, $\Xi(0)\Xi(0)^T = \Psi$ for $\Psi$ appearing in Proposition 4.(ii) and discussed in Section 3.5.

**Statistical implications.** The convergence results of Proposition 4 are the key for the convergence of the MLE in Theorem 1. Besides, we could also use these results to construct

* estimators of the local time using $S_1(n, \theta_n)/\xi_1(\theta_n)$;

* Wald confidence interval using $s(\theta_n)$ as a substitute for $s(\theta)$;

* hypothesis testing on the true value of $\theta$ using either a confidence interval around $\theta$, or $S_0(n, \theta)^2/S_1(n, \theta)$ which behaves asymptotically as a $\chi_1^2$ distribution.

## 2.1 Rate of convergence

In this section, we estimate numerically the rate of convergence towards the limit distribution of the score. More precisely, we study empirically the rate convergence of $S_k(n, \theta)$ in (4) towards its limit. In particular, we show that the speeds deteriorates as $|\theta|$ becomes close to 1.

We plot in Figure 5 the Kolmogorov-Smirnov distance $\Delta_{KS}(n)$ between the empirical distributions
<figure><img src="image_8.png" /><figcaption>Figure 5: Behavior of $n \mapsto \Delta_{\text{KS}}(n)$ using 10 000 samples of the SBM with $T = 1$, and comparison with $C/\sqrt{n}$ (dashed line), in log-log scale, for statistics (a) (MLE), (b) (0 for any $\theta$ and 2, 4 if $\theta = 0$) and (c) in the other cases.</figcaption></figure>

(a) of $\delta(\theta)^{-1}n^{1/4}(\theta_n - \theta)$ and $H/\sqrt{L_T}$ with $H \sim \mathcal{N}(0, 1)$, independent from the local time.

(b) of $S_k(n, \theta)/\Xi_{k,k}(\theta)$ and $\sqrt{L_T}H$ with $H \sim \mathcal{N}(0, 1)$, independent from the local time, for $k = 0$ or for $\theta = 0$ and $k$ even.

(c) of $S_k(n, \theta)/\xi_k(\theta)$ and $L_T$ in other cases.

Even for small moderate values of the sample size $n$, the asymptotic regime is reached most of the case. The 5%-quantile of the Kolmogorov-Smirnov statistics is 0.52 while the 1%-quantile is 0.44.

In each case, the statistics $\Delta_{\text{KS}}(n)$ behaves like $C/n^\eta$, with $\eta$ close to 1/2.

For the score of order 0 or when $\theta$ is close to 1, the rate $\eta$ is lower than 1/2. The constant $C$ is greater for the MLE than for the other Kolmogorov-Smirnov
statistics. In each case, the constant increases with $\theta$, as well as the ratio between the constant $C$ for the MLE and the one for the score of order 0.

Note that we are close to the Berry-Esseen rate even though our sample is not composed of independent identically distributed random variables.

In Figure 5, when the Kolmogorov-Smirnov distances are small and $n$ is big, the plotted distance is affected by the Monte-Carlo error. For $\theta = 0$, the terms of order 2 and 4 in the expansion of the MLE (8) vanish. The asymptotic law of these coefficients is the one given in (b). Note that, in the alternative expansion of the MLE (10), the terms of odd order of the latter expansion involve the statistics in (b) and (c). For $\theta \neq 0$, the effect of the terms of orders $k \ge 2$ in the expansion (8), which behave asymptotically as $(1 - \theta^2)^{1-k/2}$, have some increasing effect on the MLE.

# 3 The limiting coefficients

Our aim is to study the limiting coefficients $s(\theta)$ in (3), $\xi_m(\theta)$ in (5) and $\Psi$ appearing in Proposition 4. We also prove Proposition 2 and Lemma 1. We actually establish additional properties useful in numerical studies and applications. For instance, in Section 3.4, we provide an expansion of $\theta \mapsto \xi_m(\theta)$ around 0. In Section 3.5, we study the matrix $\Psi$ in Proposition 4.(ii).

## 3.1 The coefficients $\xi_k$

We study up to Section 3.4 the limiting coefficients $\xi_k$.

We take profit from the explicit expression of the density, hence of $k_{\theta}$. We conveniently rewrite

$$k_{\theta}(x, y) := \begin{cases} \frac{1}{\theta + \exp(2xy)} & \text{if } x \ge 0, y \ge 0, \\ \frac{1}{\exp(2xy) - \theta} & \text{if } x \le 0, y \le 0, \\ \frac{1}{\theta - 1} & \text{if } x \ge 0, y \le 0, \\ \frac{1}{1 + \theta} & \text{if } x \le 0, y \ge 0. \end{cases}$$

**Lemma 2.** For any $m \ge 0$,

$$\partial_{\theta}^{m} k_{\theta}(x, y) = m! (-1)^{m} k_{\theta}^{m+1}(x, y) \quad (15)$$

for any $x, y \in \mathbb{R}$.
Proof. Let us note first that

$$ 
\begin{aligned}
\partial_{\theta} k_{\theta}(x, y) &= \frac{-\operatorname{sgn}(y)^2}{(\operatorname{sgn}(y)\theta + \exp(2(xy)^+))^2} = \frac{-1}{(\operatorname{sgn}(y)\theta + \exp(2(xy)^+))^2} \\
&= -k_{\theta}(x, y)^2. \quad (16)
\end{aligned}
$$

Thus, (15) is true for $m = 0, 1$. Let us assume that (15) is true for some $m \ge 1$. With (16),

$$ 
\begin{aligned}
\partial_{\theta}^{m+1} k_{\theta}(x, y) &= (-1)^m m! \partial_{\theta}(k_{\theta}(x, y)^{m+1}) \\
&= (-1)^m (m+1)! \partial_{\theta} k_{\theta}(x, y) \cdot k_{\theta}(x, y)^m = (-1)^{m+1} (m+1)! k_{\theta}(x, y)^{m+2}.
\end{aligned}
$$

This proves that (15) is true for $m + 1$, and thus for any $m$. $\square$

**Notation 1.** We define for some integer $m \ge 1$,

$$ \chi_m^{\pm\pm}(\theta) := \iint_{\mathbb{R}^2} \mathbf{1}_{\pm x \ge 0, \pm y \ge 0} k_{\theta}(x, y)^m p_{\theta}(1, x, y) \mu_{\theta}(x) \, dx \, dy, \quad (17) $$

$$ \text{and } \chi_m(\theta) := \chi_m^{++}(\theta) + \chi_m^{--}(\theta) + \chi_m^{+-}(\theta) + \chi_m^{-+}(\theta). \quad (18) $$

With Lemma 2, the coefficients $\xi_m$'s given by (5) are related to the $\chi_m$'s by

$$ \xi_m(\theta) = m!(-1)^m \chi_{m+1}(\theta). \quad (19) $$

We now study $\chi_m(\theta)$ for $\theta \in (-1, 1)$ and various values of $m$.

We note the relations

$$ p_{\theta}(1, -x, -y) = p_{-\theta}(1, x, y), \quad (20) $$

$$ k_{\theta}(-x, -y) = -k_{-\theta}(x, y) \quad (21) $$

$$ \text{and } \mu_{\theta}(-x) = \mu_{-\theta}(x). \quad (22) $$

A direct consequence of the symmetry relations (20)-(22) is the following lemma

**Lemma 3 (Symmetry relation).** For any $m \ge 1$ and any $\theta \in (-1, 1)$,

$$ \chi_m^{--}(\theta) = (-1)^m \chi_m^{++}(-\theta) \quad \text{and} \quad \chi_m^{-+}(\theta) = (-1)^m \chi_m^{+-}(-\theta). \quad (23) $$

**Corollary 2.** For any $\theta \in (-1, 1)$ and any $m \ge 1$,

(i) $\theta \mapsto \chi_{2m}(\theta)$ is even and $\chi_{2m}(\theta) > 0$.

(ii) $\theta \mapsto \chi_{2m+1}(\theta)$ is odd; in particular $\chi_{2m+1}(0) = 0$ (we will see below in Lemma 7 that $\chi_{2m+1}(\theta) < 0$ for any $\theta \in (0, 1)$).
*Proof.* Consequence of Lemma 3 and for proving $\chi_{2m}(\theta) > 0$ one just observes the positivity of the integrands of equation (17). $\square$

**Lemma 4** (Values for $m = 1$). For any $\theta \in (-1, 1)$, $\chi_1(\theta) = 0$.

*Proof.* This stems from the very definition of $k_\theta$ and the fact that $\int_{-\infty}^{+\infty} p_\theta(1, x, y) \, dy = 1$ so that $\partial_\theta \int_{-\infty}^{+\infty} p_\theta(1, x, y) \, dy = 0$. In fact,

$$ \chi_1(\theta) = \int_{\mathbb{R}} \mathbb{E}_{\theta}[k_{\theta}(X_0, X_1) | X_0 = x] \mu_{\theta}(x) \, dx $$

and for all $x \in \mathbb{R}$

$$ \mathbb{E}_{\theta}[k_{\theta}(X_0, X_1) | X_0 = x] = \int_{\mathbb{R}} \frac{\partial_{\theta} p_{\theta}(1, x, y)}{p_{\theta}(1, x, y)} p_{\theta}(1, x, y) \, dy = 0. \quad (24) $$

Hence $\chi_1(\theta) = 0$. $\square$

**Computation for** $m \ge 2$. The case $m \ge 2$ and $\theta > 0$ is a bit more cumbersome.

**Lemma 5.** For $m \ge 2$, $\theta \in (-1, 1)$,

$$ \chi_m^{+-}(\theta) = \frac{(-1)^m}{\sqrt{2\pi}} \frac{1+\theta}{(1-\theta)^{m-1}} \quad \text{and} \quad \chi_m^{-+}(\theta) = \frac{1}{\sqrt{2\pi}} \frac{1-\theta}{(1+\theta)^{m-1}}. \quad (25) $$

*Proof.* First,

$$ \frac{1}{\sqrt{2\pi}} \iint_{\mathbb{R}^2} \mathbf{1}_{x \ge 0, y \ge 0} \exp \left( \frac{-(x+y)^2}{2} \right) \, dx \, dy = \int_0^{+\infty} \bar{\Phi}(z) \, dz = \frac{1}{\sqrt{2\pi}} $$

with $\bar{\Phi}(x) := \frac{1}{\sqrt{2\pi}} \int_x^{+\infty} e^{-y^2/2} \, dy$, the complementary distribution function of the unit, centered Gaussian random variable.

Hence,

$$ \chi_m^{+-}(\theta) = \iint_{\mathbb{R}^2} \mathbf{1}_{x \ge 0, y \le 0} k_{\theta}(x, y)^m p_{\theta}(1, x, y) \mu_{\theta}(x) \, dx \, dy = \frac{-(1+\theta)}{(1-\theta)^{m-1}} \int_0^{+\infty} \bar{\Phi}(x) \, dx. $$

The value of $\chi_m^{-+}$ follows from (23). $\square$
We express

$$
\begin{align}
\chi_m^{++}(\theta) &= \iint_{\mathbb{R}^2} \mathbf{1}_{x \ge 0, y \ge 0} k_\theta(x, y)^m p_\theta(1, x, y) \mu_\theta(x) \, dx \, dy \nonumber \\
&= \iint_{\mathbb{R}^2} \mathbf{1}_{x \ge 0, y \ge 0} \frac{1}{\sqrt{2\pi}} (1+\theta) \frac{\exp\left(\frac{-(x-y)^2}{2}\right)}{(\theta + \exp(2xy))^m} \, dx \, dy \nonumber \\
&\quad + \iint_{\mathbb{R}^2} \mathbf{1}_{x \ge 0, y \ge 0} \frac{1}{\sqrt{2\pi}} \theta(1+\theta) \frac{\exp\left(\frac{-(x+y)^2}{2}\right)}{(\theta + \exp(2xy))^m} \, dx \, dy \nonumber \\
&= \iint_{\mathbb{R}^2} \mathbf{1}_{x \ge 0, y \ge 0} \frac{1+\theta}{\sqrt{2\pi}} \exp\left(\frac{-(x+y)^2}{2}\right) \frac{1}{(\theta + \exp(2xy))^{m-1}} \, dx \, dy. \tag{26}
\end{align}
$$

We could also use the following form

$$
\chi_m^{++}(\theta) = \iint_{\mathbb{R}^2} \mathbf{1}_{x \ge 0, y \ge 0} \frac{1+\theta}{\sqrt{2\pi}} \exp\left(\frac{-x^2-y^2}{2}\right) \frac{\exp(xy) + \theta \exp(-xy)}{(\theta + \exp(2xy))^m} \, dx \, dy \quad (27)
$$

which is suitable for Monte Carlo simulation as

$$
\chi_m^{++}(\theta) = \frac{\sqrt{2\pi}}{4}(1+\theta)\mathbb{E}\left[\frac{\exp(|G \cdot G'|) + \theta \exp(-|G \cdot G'|)}{(\theta + \exp(2|G \cdot G'|))^m}\right] \quad (28)
$$

for two independent Gaussian, unit, centered random variables G and G'.

*Remark 12.* The considerations in the remainder of this section or in Lemma 1 (which we are proving in this section), suggest to consider the equivalent expression

$$
\chi_m^{++}(\theta) = \frac{\sqrt{2\pi}}{4(1-\theta)^m} (1+\theta) \mathbb{E} \left[ (1-\theta)^m \frac{\exp(|G \cdot G'|) + \theta \exp(-|G \cdot G'|)}{(\theta + \exp(2|G \cdot G'|))^m} \right]
$$

which reduces the variance of the empirical mean without burdening the computa-
tional cost.

We consider first a close-form formula for the integral involved in $\chi_m^{++}(0)$.

**Lemma 6** (Some values at $\theta = 0$). Let $m \ge 2$, then

$$
\sqrt{2\pi}\chi_m^{++}(0) = \frac{1}{2\sqrt{m(m-1)}} \log\left(2m - 1 + 2\sqrt{m(m-1)}\right) \in (0, 1).
$$

In particular, when $m$ is large, $\sqrt{2\pi}\chi_m^{++}(0) \sim \log(4m)/2m$.

*Remark 13.* We could also write $\chi_m^{++}(0) = \frac{1}{2\sqrt{m(m-1)}} \operatorname{atanh} \frac{2\sqrt{m(m-1)}}{2m-1}$.
Proof of Lemma 6. From (26) we see that $\sqrt{2\pi}\chi_m^{++}(0) = \alpha_{2m-1}$ where

$$ \alpha_n := \int_0^{+\infty} \int_0^{+\infty} \exp\left(\frac{-x^2 - y^2}{2} - nxy\right) dx dy $$

for $n \ge 1$. Let $n \ge 2$ be fixed. Using the change of variable $(x, y) = (\sqrt{v \log(u)}, \sqrt{\log(u)/v})$ for $u \ge 1, v \ge 0$, we obtain that

$$ \alpha_n = \int_0^{+\infty} \int_1^{+\infty} \frac{1}{2vu} \frac{1}{u^{n+v/2+1/2v}} du dv = \int_0^{+\infty} \frac{1}{v^2 + 2nv + 1} dv. $$

Since $n \ge 2$, from (3.3.17) in [1], this leads

$$ \alpha_n = \frac{1}{2\sqrt{n^2-1}} \log \frac{n + \sqrt{n^2-1}}{n - \sqrt{n^2-1}} = \frac{1}{\sqrt{n^2-1}} \log(n + \sqrt{n^2-1}) \in (0, 1) $$

which can be also written as $\frac{1}{\sqrt{n^2-1}} \operatorname{atanh} \frac{\sqrt{n^2-1}}{n}$. $\square$

The following inequalities are consequences of (26). For all $\theta \in [0, 1)$, $m \ge 2$:

$$ \pm \frac{1}{\sqrt{2\pi}(1 \pm \theta)^{m-2}} \sqrt{2\pi} \chi_m^{++}(0) \le \pm \chi_m^{++}(\pm\theta) \le \pm \frac{1 \pm \theta}{\sqrt{2\pi}} \sqrt{2\pi} \chi_m^{++}(0). \quad (29) $$

**Lemma 7.** For any $m \ge 1$ and $\theta \in (-1, 1)$ with $\theta \ne 0$, then $\xi_m(\theta) \ne 0$.

*Proof.* Equality (19) and Corollary 2 ensure that it suffices to prove that $(0, 1) \ni \theta \mapsto \chi_{2m+1}(\theta) < 0$. Let $m \ge 1$ and $\theta \in (0, 1)$ be fixed.

Equations (23), the explicit expression for $\chi_{2m+1}^{-}(\theta)$ in (25), the upper inequality in (29) and the fact that $\sqrt{2\pi}\chi_{2m+1}^{++}(0) < 1$ (see Lemma 6) yield

$$ \begin{aligned} \chi_{2m+1}(\theta) &= \chi_{2m+1}^{++}(\theta) - \chi_{2m+1}^{++}(-\theta) + \chi_{2m+1}^{+-}(\theta) - \chi_{2m+1}^{+-}(-\theta) \\ &\le \frac{1}{\sqrt{2\pi}} \left( 2\theta\sqrt{2\pi}\chi_{2m+1}^{++}(0) - \frac{1+\theta}{(1-\theta)^{2m}} + \frac{1-\theta}{(1+\theta)^{2m}} \right) \\ &< \frac{1}{\sqrt{2\pi}} \left( 2\theta - \frac{1+\theta}{(1-\theta)^{2m}} + \frac{1-\theta}{(1-\theta)^{2m}} \right) = \frac{2\theta}{\sqrt{2\pi}} \left( 1 - \frac{1}{(1-\theta)^{2m}} \right) < 0. \end{aligned} $$

This proves that $\xi_m(\theta)$ does not vanish for $\theta \ne 0$. $\square$

We are now ready to provide the proof of Lemma 1. The proof of Proposition 2 is also provided and requires studying $\chi_2(\theta) = -\xi_1(\theta) = s(\theta)^{-2}$.
## 3.2 Proof of Lemma 1

As (19) relates $\xi$ and $\chi$, Item (i) in Lemma 1 follows from Lemma 4. Item (ii) and (iv) are Corollary 2. We now prove Items (iii) and (v). Actually we provide a finer result.

For every $\theta \in (-1, 1)$ and $k \ge 1$, by (23) and (25), it holds

$$ \begin{aligned} \frac{\xi_{2m}(\theta)}{(2m)!} &= \chi_{2m+1}(\theta) = \chi_{2m+1}^{++}(\theta) - \chi_{2m+1}^{++}(-\theta) + \chi_{2m+1}^{+-}(\theta) - \chi_{2m+1}^{+-}(-\theta) \\ &= \chi_{2m+1}^{++}(\theta) - \chi_{2m+1}^{++}(-\theta) + \frac{1}{\sqrt{2\pi}} \frac{(1+\theta)^{2m+1} - (1-\theta)^{2m+1}}{(1-\theta^2)^{2m}} \end{aligned} $$

and

$$ \begin{aligned} -\frac{\xi_{2m-1}(\theta)}{(2m-1)!} &= \chi_{2m}(\theta) = \chi_{2m}^{++}(\theta) + \chi_{2m}^{++}(-\theta) + \chi_{2m}^{+-}(\theta) + \chi_{2m}^{+-}(-\theta) \\ &= \chi_{2m}^{++}(\theta) + \chi_{2m}^{++}(-\theta) + \frac{1}{\sqrt{2\pi}} \frac{1+\theta}{(1-\theta)^{2m}} + \frac{1}{\sqrt{2\pi}} \frac{1-\theta}{(1+\theta)^{2m}}. \end{aligned} $$

In the proof of Lemma 7 we have already obtained a bound. We find finer ones. Let us recall the bound here: for every $\theta \in (0, 1)$

$$ \frac{\xi_{2m}(\theta)}{(2m)!} = \chi_{2m+1}(\theta) < -\frac{2\theta}{\sqrt{2\pi}} \left( \frac{1}{(1-\theta)^{2m}} - 1 \right). $$

Similarly as in the proof of Lemma 7 we can obtain the desired bounds, by combining inequalities (29) and Lemma 6. For every $\theta \in (0, 1)$ we get

$$ \begin{aligned} \frac{\xi_{2m}(\theta)}{(2m)!} &\ge -\frac{1 + \sqrt{2\pi} \chi_{2m+1}^{++}(0)}{\sqrt{2\pi}} \left( \frac{1}{(1-\theta)^{2m}} - \frac{1}{(1+\theta)^{2m}} \right) \\ &\quad - \theta \frac{1 - \sqrt{2\pi} \chi_{2m+1}^{++}(0)}{\sqrt{2\pi}} \left( \frac{1}{(1-\theta)^{2m}} + \frac{1}{(1+\theta)^{2m}} \right) \end{aligned} \quad (30) $$

$$ \begin{aligned} &\ge -\frac{2}{\sqrt{2\pi}} \left( \frac{1}{(1-\theta)^{2m}} + \frac{1}{(1+\theta)^{2m}} \right) > -\frac{4}{\sqrt{2\pi}} \frac{1}{(1-\theta)^{2m}}; \\ -\frac{\xi_{2m-1}(\theta)}{(2m-1)!} &= \chi_{2m}(\theta) \le \chi_{2m}^{++}(0) \left( 1 + \theta + \frac{1}{(1-\theta)^{2(m-1)}} \right) \\ &\quad + \frac{1}{\sqrt{2\pi}} \frac{1+\theta}{(1-\theta)^{2m-1}} + \frac{1}{\sqrt{2\pi}} \frac{1-\theta}{(1+\theta)^{2m-1}}, \end{aligned} \quad (31) $$

and

$$ \begin{aligned} -\frac{\xi_{2m-1}(\theta)}{(2m-1)!} &\ge \chi_{2m}^{++}(0) \left( 1 - \theta + \frac{1}{(1+\theta)^{2(m-1)}} \right) \\ &\quad + \frac{1}{\sqrt{2\pi}} \frac{1+\theta}{(1-\theta)^{2m-1}} + \frac{1}{\sqrt{2\pi}} \frac{1-\theta}{(1+\theta)^{2m-1}}. \end{aligned} \quad (32) $$
We conclude that for all $\theta \in (-1, 1)$,

$$-\frac{4}{\sqrt{2\pi}} \frac{1}{(1 - |\theta|)^{2m-1}} < \frac{\xi_{2m-1}(\theta)}{(2m-1)!} < -\frac{1 + |\theta|}{\sqrt{2\pi}} \frac{1}{(1 - |\theta|)^{2m-1}}$$

as well as

$$-\frac{4}{\sqrt{2\pi}} \frac{1}{(1 - |\theta|)^{2m}} < \frac{\xi_{2m}(|\theta|)}{(2m)!} < -\frac{2|\theta|}{\sqrt{2\pi}} \left( \frac{1}{(1 - |\theta|)^{2m}} - 1 \right) \le 0,$$

$$0 \le \frac{2|\theta|}{\sqrt{2\pi}} \left( \frac{1}{(1 - |\theta|)^{2m}} - 1 \right) < \frac{\xi_{2m}(-|\theta|)}{(2m)!} < \frac{4}{\sqrt{2\pi}} \frac{1}{(1 - |\theta|)^{2m}}.$$

These inequalities imply Items (iii) and (v) of Lemma 1 and complete its proof.

### 3.3 Proof of Proposition 2

We prove Proposition 2 which completes the analysis of the asymptotic mixed normality in Theorem 1.

The proof of Lemma 1, in particular (31)-(32), show that for all $\theta$,

$$-\frac{2}{\sqrt{2\pi}} \frac{1 + \theta^2}{1 - \theta^2} - \chi_2^{++}(0)(2 + |\theta|) < \xi_1(\theta) < -\frac{2}{\sqrt{2\pi}} \frac{1 + \theta^2}{1 - \theta^2} - \chi_2^{++}(0)(2 - |\theta|),$$

where $\chi_2^{++}(0) \approx 0.62$ by Lemma 6. This would already prove Proposition 2, but the next lemma allow us to establish better bounds.

**Lemma 8.** Let

$$\mathcal{J}(\theta) := \int_{0}^{\infty} \int_{0}^{\infty} \frac{(1 - e^{-2xy})}{(1 - \theta^2 e^{-4xy})} e^{-4xy - \frac{(x+y)^2}{2}} \, dx \, dy.$$

Then $\mathcal{J}$ is increasing on $[0, 1]$ and $\mathcal{J}(\theta) \in [\mathcal{J}(0), \mathcal{J}(1)]$ where

$$\mathcal{J}(0) = \sqrt{2\pi}(\chi_3^{++}(0) - \chi_4^{++}(0)) \approx 0.088,$$

$$\mathcal{J}(1) = \int_{0}^{\infty} \int_{0}^{\infty} \frac{1}{(1 + e^{-2xy})} e^{-4xy - \frac{(x+y)^2}{2}} \, dx \, dy \in [\sqrt{2\pi}\chi_3^{++}(0)/2, \sqrt{2\pi}\chi_4^{++}(0)],$$

and

$$s(\theta)^{-2} = -\xi_1(\theta) = \chi_2(\theta) = \frac{2}{\sqrt{2\pi}} \left( \frac{1 + \theta^2}{1 - \theta^2} + \sqrt{2\pi} \chi_2^{++}(0) - \theta^2 \mathcal{J}(\theta) \right). \quad (33)$$

Recall that Lemma 6 shows that $\sqrt{2\pi}\chi_2^{++}(0) \approx 0.62$, $\sqrt{2\pi}\chi_3^{++}(0)/2 \approx 0.24$, and $\sqrt{2\pi}\chi_4^{++}(0) \approx 0.38$.
*Proof.* First note that $[-1, 1] \ni \theta \mapsto \mathcal{J}(\theta)$ is an even function and increasing on $[0, 1]$, therefore the first part of the statement follows from (26) and Lemma 6. We have already obtained the expressions (25) and (27) hence

$$ \chi_2(\theta) = \frac{1}{\sqrt{2\pi}} \frac{1+\theta^2}{1-\theta^2} + \int_0^\infty \int_0^\infty \frac{(1+\theta)(1-\theta e^{-2xy}) + (1-\theta)(1+\theta e^{-2xy})}{(1-\theta^2 e^{-4xy})} e^{-2xy} e^{-\frac{(x+y)^2}{2}} dy dx. $$

We rewrite the integrand in a suitable way

$$ \begin{aligned} \chi_2(\theta) &= \frac{1}{\sqrt{2\pi}} \frac{1+\theta^2}{1-\theta^2} + \frac{2}{\sqrt{2\pi}} \int_0^\infty \int_0^\infty \frac{(1-\theta^2 e^{-2xy})}{(1-\theta^2 e^{-4xy})} e^{-2xy} e^{-\frac{(x+y)^2}{2}} dy dx \\ &= \frac{1}{\sqrt{2\pi}} \frac{1+\theta^2}{1-\theta^2} + \frac{2}{\sqrt{2\pi}} \int_0^\infty \int_0^\infty \left( 1 - \frac{\theta^2 e^{-2xy}(1-e^{-2xy})}{(1-\theta^2 e^{-4xy})} \right) e^{-2xy} e^{-\frac{(x+y)^2}{2}} dy dx, \end{aligned} $$

where we recognize the right-hand-side of (33).

<figure><img src="image_9.png" /></figure>

*Remark 14.* We have obtained the following bounds for $\chi_2(\theta) = -\xi_1(\theta) = s(\theta)^{-2}$: By (33),

$$ \frac{2}{\sqrt{2\pi}} \frac{1+\theta^2}{1-\theta^2} + 2\chi_2^{++}(0) - \frac{2}{\sqrt{2\pi}} \theta^2 \mathcal{J}(0) \le \chi_2(\theta) \le \frac{2}{\sqrt{2\pi}} \frac{1+\theta^2}{1-\theta^2} + 2\chi_2^{++}(0) - \frac{2}{\sqrt{2\pi}} \theta^2 \mathcal{J}(1). $$

Since $\mathcal{J}(0) = \sqrt{2\pi}(\chi_3^{++}(0) - \chi_4^{++}(0)) \approx 0.088$ and $\mathcal{J}(1) \le \sqrt{2\pi}\chi_4^{++}(0)$, by studying the bounds above which are functions of $\theta$, we can show that

$$ \frac{s(\theta)}{\sqrt{1-\theta^2}} \in [0.79, 0.88]. $$

Numerically, we found that $s(\theta)$ is pretty close to its upper bound.

## 3.4 Expansion around 0

We provide a polynomial expansion around 0 of $\theta \mapsto \xi_m(\theta)$. By (19), it suffices to provide it for $\theta \mapsto \chi_m(\theta)$. Since we could not compute explicitly $\xi_m(\theta)$, this expansion allows us to compute it using their value at $\theta = 0$ given in Lemma 6. Let us also recall that we could rewrite $\chi_m(\theta)$ in a form suitable for Monte Carlo simulation, see (28).

Let us define

$$ \zeta_m(\theta) := \chi_m^{++}(\theta) + \chi_m^{+-}(\theta) - \chi_m^{--}(\theta) - \chi_m^{-+}(\theta). $$
Using the symmetry relations (20)-(22), $\theta \mapsto \zeta_{2m}(\theta)$ is odd, while $\theta \mapsto \zeta_{2m+1}(\theta)$ is even. In particular, $\zeta_{2m}(0) = 0$.

Besides,

$$ \chi_{2m}(0) = \beta_{2m} \text{ and } \zeta_{2m+1}(0) = \beta_{2m+1} $$

with

$$ \beta_m := \sqrt{\frac{2}{\pi}} (\sqrt{2\pi}\chi_m^{++}(0) + (-1)^m) $$

and $\chi_m^{++}(0)$ provided in Lemma 6.

**Lemma 9.** For any $m \ge 1$,

$$ \partial_{\theta}\chi_{m}(\theta) = -(m-1)\chi_{m+1}(\theta) + \zeta_{m}(\theta), \quad (34) $$

$$ \partial_{\theta}\zeta_{m}(\theta) = -(m-1)\zeta_{m+1}(\theta) + \chi_{m}(\theta). \quad (35) $$

*Proof.* Differentiating with respect to $\theta$,

$$ \begin{aligned} \partial_{\theta}\chi_{m}^{\pm\pm}(\theta) = & m \iint_{\mathbb{R}^2} \mathbf{1}_{\pm x \ge 0, \pm y \ge 0} \partial_{\theta}k_{\theta}(x, y)k_{\theta}(x, y)^{m-1}p_{\theta}(1, x, y)\mu_{\theta}(x) \, dx \, dy \\ & + \iint_{\mathbb{R}^2} \mathbf{1}_{\pm x \ge 0, \pm y \ge 0} k_{\theta}(x, y)^m \frac{\partial_{\theta}p_{\theta}(1, x, y)}{p_{\theta}(1, x, y)} p_{\theta}(1, x, y)\mu_{\theta}(x) \, dx \, dy \\ & + \iint_{\mathbb{R}^2} \mathbf{1}_{\pm x \ge 0, \pm y \ge 0} \operatorname{sgn}(x)k_{\theta}(x, y)^m p_{\theta}(1, x, y) \, dx \, dy. \end{aligned} $$

With (15) and since $k_{\theta}(x, y) = \partial_{\theta}p_{\theta}(1, x, y)/p_{\theta}(1, x, y)$,

$$ \begin{aligned} \partial_{\theta}\chi_{m}^{\pm\pm}(\theta) = & -m \iint_{\mathbb{R}^2} \mathbf{1}_{\pm x \ge 0, \pm y \ge 0} \mu_{\theta}(x)k_{\theta}(x, y)^{m+1}p_{\theta}(1, x, y) \, dx \, dy \\ & + \iint_{\mathbb{R}^2} \mathbf{1}_{\pm x \ge 0, \pm y \ge 0} \mu_{\theta}(x)k_{\theta}(x, y)^{m+1}p_{\theta}(1, x, y) \, dx \, dy \\ & + \iint_{\mathbb{R}^2} \mathbf{1}_{\pm x \ge 0, \pm y \ge 0} \operatorname{sgn}(x)k_{\theta}(x, y)^m p_{\theta}(1, x, y) \, dx \, dy. \end{aligned} $$

We then obtain

$$ \partial_{\theta}\chi_{m}^{++}(\theta) = -(m-1)\chi_{m+1}^{++}(\theta) + \chi_{m}^{++}(\theta), $$

$$ \partial_{\theta}\chi_{m}^{+-}(\theta) = -(m-1)\chi_{m+1}^{+-}(\theta) + \chi_{m}^{+-}(\theta), $$

$$ \partial_{\theta}\chi_{m}^{-+}(\theta) = -(m-1)\chi_{m+1}^{-+}(\theta) - \chi_{m}^{-+}(\theta), $$

$$ \partial_{\theta}\chi_{m}^{--}(\theta) = -(m-1)\chi_{m+1}^{--}(\theta) - \chi_{m}^{--}(\theta). $$

With these four expressions, we obtain easily (34) and (35). $\square$
**Corollary 3.** For $k \ge 1$ and $m \ge 1$,

$$ \partial_{\theta}^{k} \chi_{m}(\theta) = \sum_{\ell=0}^{k} (-1)^{\ell} b_{\ell} \binom{k}{\ell} \left( \chi_{m+\ell}(\theta) a_{\ell} + \zeta_{m+\ell}(\theta) (1-a_{\ell}) \right) \quad (36) $$

with

$$ a_{\ell} = \begin{cases} 1 & \text{if } k - \ell \text{ is even,} \\ 0 & \text{if } k - \ell \text{ is odd} \end{cases} \quad \text{and} \quad b_{\ell} = \begin{cases} 1 & \text{if } \ell = 0, \\ (m-1)\cdots(m-\ell-2) & \text{if } \ell \ge 1. \end{cases} $$

*Proof.* We form the formal power series in $\mathbb{R}[[\lambda, \mu]]$ for unknown $\lambda$ and $\mu$ by

$$ X(\theta) = \sum_{m \ge 1} \lambda^m \chi_m(\theta) + \sum_{m \ge 1} \mu^m \zeta_m(\theta). $$

We then define two linear operators $S$ and $J$ on $\mathbb{R}[[\lambda, \mu]]$ defined by $S(\nu^{m+1}) = -(m-1)\nu^m$ for $\nu = \lambda, \nu$ and $J(\lambda^m) = \mu^m, J(\mu^m) = \lambda^m$ so that $J$ exchanges $\mu$ and $\lambda$ and $J^2$ is idempotent. The operators $S$ and $J$ commute. With our choice of $S$,

$$ S(X(\theta)) = - \sum_{m \ge 1} (m-1) \lambda^m \chi_{m+1}(\theta) - \sum_{m \ge 1} \mu^m \zeta_{m+1}(\theta), $$

so that we write (34)-(35) as

$$ \partial_{\theta} X(\theta) = (S + J)X(\theta). $$

Therefore,

$$ \partial_{\theta}^{k} X(\theta) = \sum_{\ell=0}^{k} \binom{k}{\ell} S^{\ell} J^{k-\ell} X(\theta). $$

We then deduce (36). $\square$

From (36),

$$ \partial_{\theta}^{2} \chi_{m}(\theta) = (m-1)m\chi_{m+2}(\theta) - 2(m-1)\zeta_{m+1}(\theta) + \chi_{m}(\theta), $$

$$ \partial_{\theta}^{3} \chi_{m}(\theta) = -(m-1)m(m+1)\chi_{m+3}(\theta) + 3(m-1)m\zeta_{m+2}(\theta) - 3(m-1)\chi_{m+1} + \zeta_{m}(\theta). $$

The Taylor expansion around 0 leads to the following result.

**Corollary 4.** For any $\theta$ with $|\theta| < 1$, $\chi_m(\theta) = \chi_m(0) + \sum_{k \ge 1} c_{m,k} \theta^k$ where $c_{2m,1} = 0$, $c_{2m+1,1} = -2m\chi_{2m+2}(0)$ and each of the $c_{m,k}$ is a linear superposition of values of $\beta_{\ell}$ for $\ell = m, \dots, \ell + k$.

In particular,

$$ \chi_2(\theta) \approx 1.29 + 2.17 \cdot \theta^2 + \dots \text{ and } \chi_3(\theta) \approx -0.42 \cdot \theta - 2.93 \cdot \theta^2 + \dots. $$
## 3.5 The limiting covariance

We now show how to compute, at least numerically, the covariance matrix $\Psi$ appearing in Proposition 4.(ii). Recall that $\xi_{2i}(0) = 0$ for any $i \ge 0$.

According to [14, Theorem 1.2, p. 511] or [27],

$$ \Psi_{2i,2j} := \int_{\mathbb{R}} \mathbb{E}[(\partial_{\theta}^{2i} k_0(B_0, B_1) \partial_{\theta}^{2j} k_0(B_0, B_1)) | B_0 = x] dx \\ + 2 \int_{\mathbb{R}} \mathbb{E}\left[\sum_{\ell \ge 1} \partial_{\theta}^{2i} k_0(B_0, B_1) \partial_{\theta}^{2j} k_0(B_{\ell}, B_{\ell+1}) \mid B_0 = x\right] dx, $$

for a Brownian motion $B$. Owing to (15) in Lemma 2,

$$ \Psi_{2i,2j} = (2i)!(2j)! \int_{\mathbb{R}} \mathbb{E} \left[ \left. \sum_{\ell \le 0} k_0^{2i+1}(B_0, B_1) k_0^{2j+1}(B_\ell, B_{\ell+1}) \right| B_0 = x \right] dx. $$

We decompose this expression as

$$ \Psi_{2i,2j} = (2i)!(2j)!(\underline{\Psi}_{2i,2j}^I + \underline{\Psi}_{2i,2j}^{II} + \underline{\Psi}_{2i,2j}^{III}) $$

with

$$ \underline{\Psi}_{2i,2j}^I = \int_{\mathbb{R}} \mathbb{E}[k_0^{2i+1}(B_0, B_1)k_0^{2j+1}(B_0, B_1) | B_0 = x] dx, $$

$$ \underline{\Psi}_{2i,2j}^{II} = \int_{\mathbb{R}} \mathbb{E}[k_0^{2i+1}(B_0, B_1)k_0^{2j+1}(B_1, B_2) | B_0 = x] dx \\ + \int_{\mathbb{R}} \mathbb{E}[k_0^{2j+1}(B_0, B_1)k_0^{2i+1}(B_1, B_2) | B_0 = x] dx, $$

$$ \underline{\Psi}_{2i,2j}^{III} = \int_{\mathbb{R}} \mathbb{E} \left[ \left. \sum_{\ell \ge 2} k_0^{2i+1}(B_0, B_1) k_0^{2j+1}(B_\ell, B_{\ell+1}) \right| B_0 = x \right] dx \\ + \int_{\mathbb{R}} \mathbb{E} \left[ \left. \sum_{\ell \ge 2} k_0^{2j+1}(B_0, B_1) k_0^{2i+1}(B_\ell, B_{\ell+1}) \right| B_0 = x \right] dx. $$

**Lemma 10.** It holds that

$$ \underline{\Psi}_{2i,2j}^I = \chi_{2(i+j)+2}(0) = \frac{2}{\sqrt{2\pi}} (1 + \sqrt{2\pi}\chi_{2(i+j+1)}^{++}(0)) = \frac{-\xi_{2(i+j)+1}(0)}{(2(i+j)+1)!}. \quad (37) $$

$$ \text{In particular, } \Xi_{0,0}(0) = \sqrt{\underline{\Psi}_{0,0}^I} = \sqrt{-\xi_1(0)}. $$
Proof. Eq. (37) follows from the definition of $\xi_m$ in (5) and (19). When $i = j = 0$, by the Markov property and (24), we have $\underline{\Psi}_{0,0}^{II} = \underline{\Psi}_{0,0}^{III} = 0$. $\square$

By (18) and Lemma 5, $\chi_{2(i+j)+2}(0) = \frac{2}{\sqrt{2\pi}}(1 + \sqrt{2\pi}\chi_{2(i+j)+2}^{++}(0))$, and an approximation of $\chi^{++}(0)$ is given in Lemma 6.

For $m \ge 1$, $m$ odd, we set

$$ \begin{aligned} \kappa_m(x, y) &:= \exp(-2m(xy)_+) \operatorname{sgn}(y) = k_0(x, y)^m \\ \text{and } \hat{\kappa}_m(x, y) &:= \exp(-2m(xy)_+) \operatorname{sgn}(x) = \kappa_m(y, x). \end{aligned} $$

We also define

$$ \begin{aligned} K_m(x) &:= \int \kappa_m(x) \exp\left(\frac{-(x-y)^2}{2}\right) \frac{dy}{\sqrt{2\pi}} = \mathbb{E}[\kappa_m(x, B_1) | B_0 = x] \\ \text{and } \hat{K}_m(x) &:= \int \hat{\kappa}_m(x) \exp\left(\frac{-(x-y)^2}{2}\right) \frac{dy}{\sqrt{2\pi}} = \mathbb{E}[\hat{\kappa}_m(x, B_1) | B_0 = x] \end{aligned} $$

for a Brownian motion $B$. With some straightforward computations,

$$ \begin{aligned} K_m(x) &= \operatorname{sgn}(x) \exp(2x^2m(m-1))\Phi(-(2m-1)|x|) - \operatorname{sgn}(x)\Phi(-|x|) \\ \text{and } \hat{K}_m(x) &= \operatorname{sgn}(x) \exp(2x^2m(m-1))\Phi(-(2m-1)|x|) + \operatorname{sgn}(x)\Phi(-|x|). \end{aligned} $$

**Notation 2** ((Scaled) complementary error function). The complementary error function and the scaled complementary error function are

$$ \operatorname{erfc}(x) := \frac{2}{\sqrt{\pi}} \int_{x}^{+\infty} e^{-y^2} dy = 2\bar{\Phi}(x\sqrt{2}) \quad \text{and} \quad \operatorname{erfcx}(x) := \exp(x^2) \operatorname{erfc}(x) $$

for $x \ge 0$.

The functions $x \in [0, 1] \mapsto \operatorname{erfc}(x)$ and $x \in [0, 1] \mapsto \operatorname{erfcx}(x)$ take their values in $[0, 1]$. Thanks to the Mill's ratio, for $x$ large, $\operatorname{erfcx}(x) \sim 1/x\sqrt{\pi}$.

We also introduce

$$ b_m := \frac{2m-1}{\sqrt{2}} \text{ so that } 2m(m-1) = b_m^2 - \frac{1}{2}. $$

As $\Phi(-|x|) = \bar{\Phi}(|x|)$, we rewrite $K_m$ and $\hat{K}_m$ as

$$ \begin{aligned} K_m(x) &= \frac{1}{2} \operatorname{sgn}(x) e^{-x^2/2} \left( \operatorname{erfcx}(b_m|x|) - \operatorname{erfcx}(|x|/\sqrt{2}) \right) \\ \text{and } \hat{K}_m(x) &= \frac{1}{2} \operatorname{sgn}(x) e^{-x^2/2} \left( \operatorname{erfcx}(b_m|x|) + \operatorname{erfcx}(|x|/\sqrt{2}) \right). \end{aligned} $$
**Lemma 11.** For $i, j \ge 0$,

$$ \begin{aligned} \underline{\Psi}_{2i,2j}^{III} = \sum_{\ell \ge 2} & \frac{1}{\sqrt{2\pi}\sqrt{\ell-1}} \int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} \hat{K}_{2i}(y) K_{2j}(z) \exp\left(\frac{-(z-y)^2}{2(\ell-1)}\right) dy dz \\ & + \sum_{\ell \ge 2} \frac{1}{\sqrt{2\pi}\sqrt{\ell-1}} \int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} \hat{K}_{2j}(y) K_{2i}(z) \exp\left(\frac{-(z-y)^2}{2(\ell-1)}\right) dy dz. \end{aligned} $$

*Proof.* We define for $m, n \ge 0$,

$$ A(\ell, m, n) := \int_{-\infty}^{+\infty} \mathbb{E}[\kappa_m(B_0, B_1)\kappa_n(B_\ell, B_{\ell+1}) | B_0 = x] \, dx \quad (38) $$

so that

$$ \underline{\Psi}_{2i,2j}^{III} = \sum_{\ell \ge 2} 2A(\ell, 2i+1, 2j+1). $$

We rewrite $A(\ell, m, n)$ using its integral form:

$$ \begin{aligned} A(\ell, m, n) = & \frac{1}{(2\pi)^{3/2}\sqrt{\ell-1}} \int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} \kappa_m(x, y) \kappa_n(z, u) \\ & \times \exp\left(\frac{-(x-y)^2}{2}\right) \exp\left(\frac{-(z-y)^2}{2(\ell-1)}\right) \exp\left(\frac{-(z-u)^2}{2}\right) dx dy dz du. \end{aligned} \quad (39) $$

We compute first $A(\ell, m, n)$ for $\ell \ge 2$. Since $\hat{\kappa}_m(y, x) = \kappa_m(x, y)$, inverting $x$ and $y$, and using the definitions of $K_m$ and $\hat{K}_m$ leads to

$$ A(\ell, m, n) = \frac{1}{\sqrt{2\pi}\sqrt{\ell-1}} \int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} \hat{K}_m(y) K_n(z) \exp\left(\frac{-(z-y)^2}{2(\ell-1)}\right) dy dz. $$

This gives the result.

<figure><img src="image_10.png" /></figure>

*Remark 15.* With a change of variable $(y, z) \rightarrow (y/\sqrt{j-1}, z/\sqrt{j-1})$, for $A$ introduced in (38),

$$ \begin{aligned} A(\ell, m, n) = & \sqrt{2\pi}\sqrt{\ell-1} \int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} \hat{K}_m(y\sqrt{\ell-1}) K_n(z\sqrt{\ell-1}) \exp(yz) \\ & \times \exp\left(\frac{-z^2}{2}\right) \exp\left(\frac{-y^2}{2}\right) \frac{dy}{\sqrt{2\pi}} \frac{dz}{\sqrt{2\pi}}. \end{aligned} $$

For two independent Gaussian random variables $G, G' \sim \mathcal{N}(0, 1)$, we may write

$$ A(\ell, m, n) = \sqrt{2\pi}\sqrt{\ell-1} \mathbb{E}[K_m(G'\sqrt{\ell-1})\hat{K}_n(G\sqrt{\ell-1})\exp(G \cdot G')], $$

where $G, G' \sim \mathcal{N}(0, 1)$ are independent. However, the variance is high so that Monte Carlo techniques are unstable here.
**Lemma 12.** For $i, j \ge 1$,

$$ \underline{\Psi}_{2i,2j}^{II} = \int_{0}^{+\infty} e^{-x^2} \operatorname{erfcx}(b_{2i+1}x) \operatorname{erfc}(b_{2j+1}x) \, dx - \frac{2}{\sqrt{\pi}} (\sqrt{2} - 1). $$

*Proof.* As above,

$$ \begin{aligned} A(1, m, n) &= \int_{-\infty}^{+\infty} \hat{K}_m(x) K_n(x) \, dx \\ &= \frac{1}{2} \int_{0}^{+\infty} e^{-x^2} \operatorname{erfcx}(b_m x) \operatorname{erfc}(b_n x) \, dx - \frac{1}{2} \int_{0}^{+\infty} e^{-x^2} \operatorname{erfcx}(x/\sqrt{2})^2 \, dx. \end{aligned} $$

The quantity $\underline{\Psi}_{2i,2j}^{III}$ is then

$$ \underline{\Psi}_{2i,2j}^{II} = A(1, 2i + 1, 2j + 1) + A(1, 2j + 1, 2i + 1). $$

According to [29, 4.7.2],

$$ \int_{0}^{+\infty} e^{-x^2} \operatorname{erfcx}\left(\frac{x}{\sqrt{2}}\right)^2 dx = \int_{0}^{+\infty} \operatorname{erfc}\left(\frac{x}{\sqrt{2}}\right)^2 dx = \frac{2}{\sqrt{\pi}}(\sqrt{2}-1) \approx 0.467. $$

No closed-form expression for $\int_{-\infty}^{+\infty} e^{-x^2} \operatorname{erfcx}(b_m x) \operatorname{erfc}(b_n x) dx$ seems to exist. However, this integral is easy to compute numerically as there exists various implementation of erfcx. This could be done through a quadrature method to compute the integral or using a Monte Carlo method as

$$ \int_{0}^{+\infty} e^{-x^2} \operatorname{erfcx}(b_m x) \operatorname{erfc}(b_n x) dx \\ = \frac{\sqrt{\pi}}{\sqrt{2m-1}} \mathbb{E} \left[ \operatorname{erfcx} \left( \frac{\sqrt{2m-1}}{\sqrt{2}} |G| \right) \operatorname{erfcx} \left( \frac{\sqrt{2n-1}}{\sqrt{2}} |G| \right) \right]. $$

This concludes the proof. $\square$

**Numerical computations.** We give some the first values of $\Psi$ obtained by numerical computations:

$$ \{\Psi_{2i,2j}\}_{i,j=0,1,2} \approx \begin{bmatrix} 1.295 & 1.422 & 12.576 \\ 1.422 & 1.891 & 18.135 \\ 12.576 & 18.135 & 181.421 \end{bmatrix}. $$

Let us conclude by an observation: for $m \ge 2$, $m$ odd, we found the empirical rule of thumbs that $\underline{\Psi}_{m,m}/\underline{\Psi}_{m,m}^I$ behaves, for moderate values of $m$, as $C_m/\sqrt{2m}$ for a pre-factor $C_m$ varying slowly. In addition, $\underline{\Psi}_{m,m}^I$ converges to $\frac{2}{\sqrt{2\pi}}$ as $m \to \infty$ (see Lemma 6 and Lemma 10) while $\underline{\Psi}_{m,m}^{II}$ converges to $-2(\sqrt{2}-1)/\sqrt{\pi}$.
# 4 Proofs of the asymptotic results

The goal of this section is to prove Proposition 4. Proposition 3 is proved in Section 2 as a corollary of Proposition 4. Theorem 2 is a direct application of [20, Theorem 3]. The convergence in Theorem 1 follows from combining the expansion in Theorem 2 with Proposition 3.

Proposition 4 is a consequence of the results stated in Remark 10, in particular of the convergence (14) which is a simple application of [27]. The next Lemma 13 is an integrability condition which ensures that the assumptions of the main result in [27] are satisfied. Although (14) concerns derivatives of $k_\theta$, the lemma concerns its powers. This is justified by Lemma 2 in Section 3 where it is shown that $\partial_\theta^m k_\theta(x, y)$ is proportional to $k_\theta(x, y)^m$ with a multiplicative factor depending only on $m \in \mathbb{N}$.

**Lemma 13.** For $m = 1, 2, 3, \dots$, and any $\gamma > 3$, there exists $a > 0$ such that for every $\theta \in (-1, 1)$ there exists a measurable, bounded function $h_\theta: \mathbb{R} \to \mathbb{R}$ which satisfies

$$ \int_{\mathbb{R}} h_{\theta}(x) |x|^{\gamma} dx < +\infty \quad \text{and} \quad |k_{\theta}^{m}(x, y)| \le h_{\theta}(x) e^{a|y-x|}, \quad \text{for any } x, y \in \mathbb{R}. $$

*Proof.* Fix $a > 0$. Observe that

$$ |k_{\theta}(x, y)| \le \frac{1}{(1 - |\theta|)} \mathbf{1}_{\{xy \le 0\}} + \frac{1}{(1 - |\theta|)} e^{-2xy} \mathbf{1}_{\{xy > 0\}} =: h_{\theta}(x, y) \le \frac{1}{(1 - |\theta|)}. $$

Now, $\mathbf{1}_{\{xy \le 0\}} \le e^{-a|x|}e^{a|y-x|}$. Let us find a bound for $e^{-2xy}e^{-a|y-x|}\mathbf{1}_{\{xy > 0\}}$. Let us first assume $xy > 0$ and $|y| \ge |x|$ then

$$ e^{-2xy}e^{-a|y-x|} = e^{-2|x||y|}e^{-a|y|+a|x|} \le e^{-2x^2} $$

and if $xy > 0$ and $|y| < |x|$ then

$$ e^{-2xy}e^{-a|y-x|} = e^{-2|x||y|}e^{-a|x|+a|y|} \le e^{-2x^2} \mathbf{1}_{\{|x| \le a/2\}} + e^{-a|x|} \mathbf{1}_{\{|x| > a/2\}}. $$

We conclude that for any $x, y \in \mathbb{R}$,

$$ |k_{\theta}(x, y)|^m \le (1 - |\theta|)^{-(m-1)} h_{\theta}(x, y) \le (1 - |\theta|)^{-m} (e^{-2x^2} + e^{-a|x|}) e^{a|y-x|} $$

$$ \text{with } \int_{\mathbb{R}} (e^{-2x^2} + e^{-a|x|}) |x|^{\gamma} dx < +\infty. \quad \square $$

<figure><img src="image_11.png" /></figure>
### 4.0.1 Item (i): Case $\theta \neq 0$

Item (i) of Proposition 4 is an immediate consequence of the result stated in Remark 10, Remark 11, and Lemma 1 which ensures that when $\theta \in (-1, 1)$, $\xi_0(\theta) = 0$, and when in addition $\theta \neq 0$, $\xi_m(\theta) \neq 0$ as soon as $m \ge 1$ and $\xi_{2m}(0) = 0$.

Using the fact that given two sequences, one converging in probability and the other converges stably, joint stable convergence holds (see [2, Theorem 1]), we get the joint stable convergence of the vector $(S_k(\theta))_{k=0,...,m}$ for any $m \ge 0$.

### 4.0.2 Item (ii): Case $\theta = 0$

From the above arguments we have only to focus on the convergence of

$$S\sim(n, 0) := (n^{1/4}S_0(n, 0), n^{1/4}S_2(n, 0), \dots, n^{1/4}S_{2m}(n, 0))$$

since $S_{2m+1}(n, 0)$ converges in probability to $\xi_{2m+1}(\theta)L_T$ for each $m \ge 0$.

The result is a direct consequence of [14, Theorem 1.2, p. 511] (which can be applied since Lemma 13 holds).

## References

[1] Milton Abramowitz and Irene A. Stegun. Handbook of mathematical functions with formulas, graphs, and mathematical tables, volume 55 of National Bureau of Standards Applied Mathematics Series. 1964.

[2] D. J. Aldous and G. K. Eagleson. On mixing and stability of limit theorems. Ann. Probab., 6:325-331, 1978.

[3] Héctor Araya, Meryem Slaoui, and Soledad Torres. Bayesian inference for fractional oscillating Brownian motion. Comput. Statist., 37(2):887-907, 2022.

[4] Yizhou Bai, Yongjin Wang, Haoyan Zhang, and Xiaoyang Zhuo. Bayesian estimation of the skew Ornstein-Uhlenbeck process. Computational Economics, 60(2):479-527, July 2021.

[5] Manuel Barahona, Laura Rifo, Maritza Sepúlveda, and Soledad Torres. A simulation-based study on Bayesian estimators for the skew Brownian motion. Entropy, 18(7):Paper No. 241, 14, 2016.

[6] Olivier Bardou and Miguel Martinez. Statistical estimation for reflected skew processes. Stat. Inference Stoch. Process., 13(3):231-248, 2010.
[7] Said Karim Bounebache and Lorenzo Zambotti. A skew stochastic heat equation. *J. Theor. Probab.*, 27(1):168–201, 2014.

[8] R.S. Cantrell and C. Cosner. Diffusion models for population dynamics incorporating individual behavior at boundaries: Applications to refuge design. *Theoretical Population Biology*, 55(2):189–207, 1999.

[9] Marc Decamps, Marc Goovaerts, and Wim Schoutens. Self exciting threshold interest rates models. *Int. J. Theor. Appl. Finance*, 9(7):1093–1122, 2006.

[10] Alexander Gairat and Vadim Shcherbakov. Density of skew Brownian motion and its functionals with application in finance. *Math. Finance*, 27(4):1069–1088, 2017.

[11] J. M. Harrison and L. A. Shepp. On skew Brownian motion. *Ann. Probab.*, 9(2):309–313, 1981.

[12] I.A. Ibragimov and R.Z. Has'minskii. *Statistical Estimation Asymptotic Theory*. Springer, 1981.

[13] Andrey Itkin, Alexander Lipton, and Dmitry Muravey. Multilayer heat equations and their solutions via oscillating integral transforms. *Physica A*, 601:34, 2022. Id/No 127544.

[14] J. Jacod. Rates of convergence to the local time of a diffusion. *Ann. Inst. H. Poincaré Probab. Statist.*, 34(4):505–544, 1998.

[15] Yu.A. Kutoyants. *Parameter Estimation for Stochastic Processes*. Heldermann, 1984.

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

[17] A. Lejay and P. Pigato. Statistical estimation of the oscillating brownian motion. *Bernoulli*, 24(4B):3568–3602, 2018.

[18] Antoine Lejay. Estimation of the bias parameter of the skew random walk and application to the skew Brownian motion. *Stat. Inference Stoch. Process.*, 21(3):539–551, 2018.

[19] Antoine Lejay. Estimation of the bias parameter of the skew random walk and application to the skew Brownian motion. *Stat. Inference Stoch. Process.*, 21(3):539–551, 2018.
[20] Antoine Lejay and Sara Mazzonetto. Beyond the delta method, 2022. Preprint https://hal.inria.fr/hal-03738371v1.

[21] Antoine Lejay, Ernesto Mordecki, and Soledad Torres. Is a Brownian motion skew? *Scand. J. Stat.*, 41(2):346-364, 2014.

[22] Antoine Lejay, Ernesto Mordecki, and Soledad Torres. Two consistent estimators for the skew brownian motion. *ESAIM: Probability and Statistics*, 23:567-583, 2019.

[23] Antoine Lejay and Géraldine Pichot. Simulating diffusion processes in discontinuous media: a numerical scheme with constant time steps. *J. Comput. Phys.*, 231(21):7299-7314, 2012.

[24] Antoine Lejay and Paolo Pigato. Statistical estimation of the oscillating Brownian motion. *Bernoulli*, 24(4B):3568-3602, 2018.

[25] Antoine Lejay and Paolo Pigato. A threshold model for local volatility: evidence of leverage and mean reversion effects on historical data. *Int. J. Theor. Appl. Finance*, 22(4):1950017, 24, 2019.

[26] Antoine Lejay and Paolo Pigato. Maximum likelihood drift estimation for a threshold diffusion. *Scand. J. Stat.*, 47(3):609-637, 2020.

[27] Sara Mazzonetto. Rates of convergence to the local time of oscillating and skew brownian motions, 2019.

[28] Pedro P. Mota and Manuel L. Esquível. On a continuous time stock price model with regime switching, delay, and threshold. *Quant. Finance*, 14(8):1479-1488, 2014.

[29] Edward W. Ng and Murray Geller. A table of integrals of the error functions. *J. Res. Nat. Bur. Standards Sect. B*, 73B:1-20, 1969.

[30] Otso Ovaskainen and Stephen J. Cornell. Biased movement at a boundary and conditional occupancy times for diffusion processes. *J. Appl. Probab.*, 40(3):557-580, 2003.

[31] Jorge M. Ramirez, Enrique A. Thomann, Edward C. Waymire, Roy Haggerty, and Brian Wood. A generalized Taylor-Aris formula and skew diffusion. *Multiscale Model. Simul.*, 5(3):786-801, 2006.

[32] Christian Y. Robert. How large is the jump discontinuity in the diffusion coefficient of a time-homogeneous diffusion? *Econometric Theory*, page 1-33, 2022.
[33] Fei Su and Kung-Sik Chan. Quasi-likelihood estimation of a threshold diffusion process. *J. Econometrics*, 189(2):473–484, 2015.

[34] Fei Su and Kung-Sik Chan. Testing for threshold diffusion. *J. Bus. Econom. Statist.*, 35(2):218–227, 2017.

[35] J. B. Walsh. A diffusion with discontinuous local time. In *Temps locaux*, volume 52-53, pages 37-45. Société Mathématique de France, 1978.

[36] M. Zhang. Calculation of diffusive shock acceleration of charged particles by skew Brownian motion. *The Astrophysical Journal*, 541:428–435, 2000.