# Simulation of reflected Brownian motion on two dimensional wedges

Pierre Bras, Arturo Kohatsu-Higa

► To cite this version:

Pierre Bras, Arturo Kohatsu-Higa. Simulation of reflected Brownian motion on two di-
mensional wedges. Stochastic Processes and their Applications, 2023, 156, pp.349-378.
10.1016/j.spa.2022.11.011. hal-03898480

HAL Id: hal-03898480

https://hal.science/hal-03898480

Submitted on 14 Dec 2022

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Simulation of Reflected Brownian motion on two dimensional wedges

Pierre Bras* and Arturo Kohatsu-Higa†

**Abstract**

We study a correlated Brownian motion in two dimensions, which is reflected, stopped or killed in a wedge represented as the intersection of two half spaces. First, we provide explicit density formulas, hinted by the method of images. These explicit expressions rely on infinite oscillating sums of Bessel functions and may demand computationally costly procedures.

We propose suitable recursive algorithms for the simulation of the laws of reflected and stopped Brownian motion which are based on generalizations of the reflection principle in two dimensions. We study and give bounds for the complexity of the proposed algorithms.

**Key words:** Reflected Brownian motion, Wedge, Hitting times, Reflection principle, Method of images, Monte Carlo simulation.

**MSC 2020:** 60H35, 65C05, 65C30

# 1 Introduction

The objective of the present article is to provide exact formulas and algorithms for the simulation of a normally reflected two-dimensional Brownian motion starting at $x_0 \in \mathbb{R}^2$.

The reflected process is denoted by $X = (X_t)_{t \in [0,T]}$ and the domain of reflection is a wedge $\mathcal{D}$, i.e. a subset of $\mathbb{R}^2$ delimited by two (non parallel) lines, so as to give Monte Carlo methods for the estimation of $E^{x_0}[f(X_T)]$, where $f: \mathcal{D} \to \mathbb{R}$ is a measurable function such that $f(X_T) \in L^1$ and $T > 0$ is the finite time horizon. Our first goal is to obtain an explicit formula for the density of $X_t$ (see Theorem 3). Unfortunately, this formula involves oscillating infinite sums of Bessel functions. Instead of directly computing these sums, we propose an alternative simulation method which uses an extension of the reflection principle in two dimensions for a particular type of wedges with angle $\pi/m$, $m \in \mathbb{N}$. As a first step, we obtain a simulation method for $E^{x_0}[f(W_{T \wedge \tau})]$, where $W = (W_t)_{t \ge 0}$ is a two-dimensional Brownian motion and $\tau$ stands for the first time the process $W$ touches the boundary of the wedge $\mathcal{D}$. Applying the methodology for the stopped process recursively one obtains an algorithm for the reflected process $X$. The algorithms we present here rely on an idea which has links to the rectangle method [DL06]: we include a smaller wedge of angle $\pi/m$ inside the first wedge and simulate the exit from the smaller wedge.

The algorithm for the simulation of the stopped process improves the method proposed in [Met09]. In the reflected case, we prove that although some lower moments are finite, the average number of simulations for this method is infinite. This effect is created due to the large number of reflections that may happen if the process enters the corner of the wedge. We then propose a modified approximation scheme which takes into account the asymptotic behaviour of the density near the corner. With this modification, we show that the expected number of iterations of the algorithm is finite and measure the error of approximation in total variation distance. We then give adaptations of our algorithms to more general Itô processes reflected in a wedge.

The expression of the density of the stopped Brownian motion in a wedge of $\mathbb{R}^2$ has been obtained in [Iye85] using the method of images. However, the author does not provide full arguments for the verification of the initial condition and some exchanges of infinite summation and integrations are not fully explained. For the history of the expression of the density of the stopped Brownian motion on a wedge, see [BnS97].

<footnote>* Sorbonne Université, Laboratoire de Probabilités, Statistique et Modélisation, UMR 8001, case 188, 4 pl. Jussieu, F-75252 Paris Cedex 5, France. E-mail: pierre.bras@sorbonne-universite.fr. Supported in part by École Normale Supérieure, Département de Mathématiques et Applications.</footnote>

<footnote>†Ritsumeikan University. Department of Mathematical Sciences. E-mail: khts00@fc.ritsumei.ac.jp. Supported in part by KAKENHI 20K03666.</footnote>
In [CBM15] and [Met10] the authors correct some mistakes in formulas appearing in [Iye85] and give formulas for a number of other random variables related to the stopped Brownian motion in two dimensions with general correlation coefficients, such as the survival probability and the first-passage time distribution. However, formulas from [Iye85] and [Met10] are not directly applicable to simulation algorithms of stopped processes, as they involve infinite sums and Bessel functions. [Met09, Section 2] gives an approximate simulation algorithm using the semi-analytic expression of the density established in [Iye85]. But this method implies the approximation of an infinite sum of Bessel functions and it leads to bias when simulating the stopping time $\tau$.

More recent papers ([EFW13], [BSCD13]) extend the results of [Iye85] and [Met10] to three dimensions and give algorithms and density formulas for the stopped Brownian motion in three dimensions, using the method of images in $\mathbb{R}^3$. These algorithms which point towards the possibility of larger dimension extensions are however limited to some specific values of correlation. In [KLR18], the authors use the spectral decomposition method which applies in a larger generality but still relies on efficient computational methods for eigenvalue calculations and truncations of infinite sums. These references concentrate on the stopped Brownian motion case and there is no discussion about the simulation of reflected Brownian motion in wedges.

Simulation algorithms for the reflected Brownian motion have been only partly studied, although they share similarities with the stopped case. Simulation algorithms for quantities related to reflected Brownian motion on an orthant in multi-dimensions can be found in [BM18] and [BC15]. They do not use the method of images, but the characterization through the Skorokhod problem and the so-called $\varepsilon$-strong methodology which relies on the regularity of the reflection functional. Still, the method proposed has a theoretically infinite mean running time.

Other approximations methods which are closer to random walk versions of the reflected Brownian motion on tile domains are studied in [Dub04] and [Kag07]. Furthermore, other simulations methods for reflected Brownian motion which do not cover the case of the wedge can be found in [CPS98], [BGT04], [BC08] and [BST10] the references therein.

The reflected Brownian motion and more generally reflected processes have applications in finance (for example, in stochastic models where the process, which can model an interest rates, is constrained to be non-negative [Ha09] or has other constraints like barriers such as in [IKP13]); in queueing models [GNR86], etc. In particular, in the recent years there has been a lot of developments concerning the study of stationary measures of reflected Brownian motion such as [DM09] which is slightly related to the ideas which we use here to construct a simulation method on a wedge. For a review, on this topic and the relationship with queuing theory see e.g. [Die10].

Considering the particular case of a wedge may open the way to new simulation algorithms for reflected processes adapted to non-smooth domains with corner points in two or more dimensions (see Section C).

The paper is organized as follows. In Section 3 we state the framework of the problem, giving a parametrization of the wedge in $\mathbb{R}^2$. We also give the formula for the density of the stopped Brownian motion in a wedge of angle $\alpha \in (0, 2\pi)$. The density of the stopped Brownian motion had already been obtained in [Iye85]. A more general case is studied in [BnS97] which provides a detailed history and a full proof that the density formula is the fundamental solution of the associated partial differential equation.

In Section 4, we provide the density of the reflected Brownian motion on the wedge. First, we use the method of images in a special type of wedges as it serves to understand how one induces the density of the reflected Brownian motion in the general case. In the Appendix we include the proof in the reflected case on more general cones following [BnS97].

In Section 5, we give explicit algorithms for the simulation of the stopped and reflected Brownian motions. To do so, we first derive from [Met10] and the formulas proved in Section 4 expressions of densities of the exit time and exit point from the wedge. Then, we give an algorithm for the simulation of the stopped Brownian motion using the explicit expressions in the case $\alpha = \pi/m$. The simulation algorithm is based on the simulation of a random sequence of stopped Brownian motion processes on domains with angle $\alpha = \pi/m$. This first algorithm is useful for the simulation of reflected Brownian motion, but it also has its own interest. Then we give the simulation algorithm for the reflected Brownian motion with unit diffusion matrix.

In Section 6, we study the mean number of iterations of the algorithm to finish, denoted $\mathbb{E}[N]$. For the stopped Brownian motion, $\mathbb{E}[N]$ is bounded above by a constant. The reflected case is more difficult to deal with. In this case, $N$ is equal to the number of times a Brownian motion goes from one frontier of a wedge to another. We first give theoretical bounds for $\mathbb{E}[N]$ and prove that although the number of iterations is almost surely finite, $\mathbb{E}[N]$ is infinite, which is due to the events when the Brownian motion comes close to the origin. We then propose a modification of the algorithm: we stop the algorithm when the Brownian motion comes too close
to the origin, measured by a parameter $\epsilon$ and use an approximation of the density for the reflected Brownian motion close to the origin, induced by the density formula in Section 4. We then prove that for the modified algorithm, $E[N]$ has an upper bound which grows as $\epsilon^{1-p}$ for all $p \in (1, 2)$, and that the error of approximation is of order $\epsilon$.

In Section 7, we show how the algorithms can directly be applied for the simulation of the reflected and stopped Brownian motions with non-zero drift. Then we give an algorithm for the simulation of a class of Itō diffusion processes, reflected or stopped in a wedge in $\mathbb{R}^2$. This algorithm is in fact a direct application of the precedent algorithms, combined with a Euler-Maruyama scheme.

Finally, in Section 8, we perform Monte Carlo simulations for the estimation of reflected and stopped Brownian motions and Itō processes. We show that the bias in the algorithm from [Met09] can be significant compared with our method which is exact. Simulations prove the need for the approximation algorithm for the reflected Brownian motion, as the exact algorithm takes too much time to be efficient. Clearly, one remaining subject which is not considered here is the reflecting case with non-unit diffusion coefficient. This remains an open problem.

The authors would like to express thanks to the anonymous referee which provided references [BnS97] and [DM09] which improved the results and shortened the proof of Theorem 3.

## 2 Notations

We give a brief list of notations that are used through the text.

The natural logarithm is denoted with log. We frequently use the notation $\pm$ in order to denote the two rays that define a wedge. Sometimes rather than writing two equations we just write one using the symbol $\pm$ or $\mp$ meaning that we are stating two equations, one using the top symbols and another using the bottom symbols that appear throughout the equation.

$I_a$ stands for the modified Bessel function of the first kind given for $a \ge 0$, $x \ge 0$ by:

$$I_a(x) = \sum_{k=0}^{\infty} \frac{x^{a+2k}}{2^{a+2k} k! \Gamma(k+a+1)}. \qquad (2.1)$$

This function satisfies the differential equation

$$x^2 \frac{d^2 I_a(x)}{dx^2} + x \frac{dI_a(x)}{dx} - (x^2 + a^2)I_a(x) = 0, \quad x \ge 0. \qquad (2.2)$$

We consider the space $\mathbb{R}^2$ endowed with the Euclidean norm denoted by $|\cdot|$.

For $\mathcal{D} \subset \mathbb{R}^2$ and $k \in \mathbb{N}$, we denote by $C_b^k(\mathcal{D})$ the set of real-valued functions defined on $\mathcal{D}$ which are bounded and have bounded partial derivatives up to the order $k$.

Now we define the probabilistic setting as follows. Let $(\Omega, \mathcal{F}, \mathbb{P}, (\mathcal{F}_t)_{t \ge 0})$ be a filtered probability space satisfying the usual conditions. Let $T > 0$ be the time horizon and let $W = (W_t)_{t \ge 0}$ be a two-dimensional correlated Brownian motion on $(\Omega, \mathcal{F}, \mathbb{P}, (\mathcal{F}_t)_{t \ge 0})$. We denote by $\Sigma$ the covariance matrix of $W$ and we assume that $\Sigma$ is a non-singular matrix. We define a family of probabilities on $(\Omega, \mathcal{F}, (\mathcal{F}_t)_{t \ge 0})$ by:

$$\forall x_0 \in \mathbb{R}^2, \mathbb{P}^{x_0} := \mathbb{P}_{|W_0=x_0}.$$

We will denote by $E^{x_0}$ the expectation under $\mathbb{P}^{x_0}$.

If $W = (W_t)_{t \ge 0}$ is a Brownian motion we will say “stopped Brownian motion” to designate the process $(W_{t \wedge \tau})_{t \ge 0}$, and “killed Brownian motion” to designate the process $(W_t \mathbb{1}_{\tau > T})_{t \in [0, T]}$, where $\tau$ is some stopping time defined later.

We will frequently use in the proofs the Bessel process $(R_t)_{t \ge 0}$ in dimension $\delta \ge 2$, i.e. a real-valued process satisfying the stochastic differential equation:

$$R_t = r_0 + B_t + \int_0^t \frac{\delta - 1}{2} \frac{ds}{R_s},$$

where $(B_t)_{t \ge 0}$ is a standard one-dimensional Brownian motion starting from zero. The index of the Bessel process is defined as $\nu = \frac{\delta}{2} - 1$. For a general reference to these matters we refer to [JYC09, Section 6].
<figure><img src="image_2.png" /><figcaption>Figure 3.1: Example of a wedge of angle $\alpha$.</figcaption></figure>

# 3 Setting of the problem

Let $a \in \mathbb{R} \cup \{\pm\}$ and consider the subset $\mathcal{D} \subset \mathbb{R}^2$ defined by one of the four following cartesian equations:

$$ \mathcal{D} = \left\{ \begin{array}{ll} \{(x, y) \in \mathbb{R}^2, y \ge 0 \text{ and } y \le ax\} & \text{with } a > 0, \\ \{(x, y) \in \mathbb{R}^2, y \ge 0 \text{ and } y \ge ax\} & \text{with } a < 0, \\ \{(x, y) \in \mathbb{R}^2, y \ge 0 \text{ or } y \ge ax\} & \text{with } a > 0, \\ \{(x, y) \in \mathbb{R}^2, y \ge 0 \text{ or } y \le ax\} & \text{with } a < 0. \end{array} \right. \qquad (3.1) $$

The set $\mathcal{D}$ is called a wedge<fnref n="1" />. The angle of the wedge is denoted by $\alpha$ and chosen such that $\alpha \in (0, 2\pi)$. An example is given in Figure 3.1. We write the wedge $\mathcal{D}$ in polar coordinates as follows:

$$ \mathcal{D} = \{(r \cos(\theta), r \sin(\theta)), r \in \mathbb{R}^{+}, \theta \in [0, \alpha]\}. $$

We define its boundary as $\partial\mathcal{D} = \partial\mathcal{D}^- \cup \partial\mathcal{D}^+$, where

$$ \partial\mathcal{D}^- := \{(x, y) \in \mathcal{D}, y = 0\} = \{(r \cos(\theta), r \sin(\theta)), r \in \mathbb{R}^{+}, \theta = 0\}, $$

$$ \partial\mathcal{D}^+ := \{(x, y) \in \mathcal{D}, y = ax\} = \{(r \cos(\theta), r \sin(\theta)), r \in \mathbb{R}^{+}, \theta = \alpha\}. $$

In general, we use the notation $\mathcal{D} = \langle\alpha\rangle$ to denote a wedge which boundary is determined by the rays $\partial\mathcal{D}^\pm$. In a similar way, we also generalize this notation to a general wedge which boundaries are determined by the rays at the angles $0 \le \alpha < \beta < 2\pi$ as $\mathcal{V} = \langle\alpha, \beta\rangle$.

We use the following stopping times:

$$ \tau = \inf \{t > 0, W_t \in \partial\mathcal{D}\}. $$

$$ \tau^{\pm} = \inf \{t > 0, W_t \in \partial\mathcal{D}^{\pm}\}. $$

As in [Met10], we parametrize the covariance matrix of $W$ via its Cholesky decomposition and write $\Sigma = \sigma\sigma^T$ with

$$ \sigma = \begin{pmatrix} \sigma_1 \sqrt{1 - \rho^2} & \sigma_1 \rho \\ 0 & \sigma_2 \end{pmatrix}, $$

where $\rho \in (-1, 1)$ and $\sigma_1, \sigma_2 \ge 0$. Then assuming that $\sigma$ is invertible, we consider the process $W' := \sigma^{-1}W$, which is a standard two-dimensional Brownian motion with independent components. Furthermore, using the explicit formula for the inverse matrix, we deduce that

$$ W^2 = 0 \iff W'^2 = 0, $$

$$ W^2 = aW^1 \iff \begin{cases} W'^2 = \frac{a\sigma_1\sqrt{1-\rho^2}}{\sigma_2-a\sigma_1\rho}W'^1 & \text{if } \sigma_2 - a\sigma_1\rho \neq 0, \\ W'^1 = 0 & \text{if } \sigma_2 - a\sigma_1\rho = 0. \end{cases} $$

<footnote n="1">¹The cases where $a = 0$ are not considered here because they lead to cases where the problem can be simplified to a one dimensional problem.</footnote>
<table><caption>Figure 3.2: Table summing up the four different cases of wedge $\mathcal{D}$</caption><thead><tr><th>Cartesian Equation</th><th>{y &ge; 0 and y &le; ax}, a &gt; 0</th><th>{y &ge; 0 and y &ge; ax}, a &lt; 0</th><th>{y &ge; 0 or y &ge; ax}, a &gt; 0</th><th>{y &ge; 0 or y &le; ax}, a &lt; 0</th></tr></thead><tbody><tr><th>Value of &alpha;</th><td>arctan(a)</td><td>&pi; + arctan(a)</td><td>&pi; + arctan(a)</td><td>2&pi; + arctan(a)</td></tr><tr><th>Wedge for the original problem</th><td></td><td></td><td></td><td></td></tr><tr><th>Equation of W' if &sigma;<sub>2</sub> &minus; a&sigma;<sub>1</sub>&rho; &gt; 0</th><td>{y &ge; 0 and y &le; a'<sub>1</sub>x}, a' &gt; 0</td><td>{y &ge; 0 and y &ge; a'<sub>1</sub>x}, a' &lt; 0</td><td>{y &ge; 0 or y &ge; a'<sub>1</sub>x}, a' &gt; 0</td><td>{y &ge; 0 or y &le; a'<sub>1</sub>x}, a' &lt; 0</td></tr><tr><th>Equation of W' if &sigma;<sub>2</sub> &minus; a&sigma;<sub>1</sub>&rho; &lt; 0</th><td>{y &ge; 0 and y &ge; a'<sub>1</sub>x}, a' &lt; 0</td><td>{y &ge; 0 and y &le; a'<sub>1</sub>x}, a' &gt; 0</td><td>{y &ge; 0 or y &le; a'<sub>1</sub>x}, a' &lt; 0</td><td>{y &ge; 0 or y &ge; a'<sub>1</sub>x}, a' &gt; 0</td></tr></tbody></table>

In the case that $\sigma_2 - a\sigma_1\rho = 0$, then the problem reduces to two independent one-dimensional Brownian motions being reflected or stopped in the first quadrant, and therefore the problem can be reduced to using two independent reflected or stopped Brownian motions in one dimension. For explicit formulas, see (4.2) and (4.1). Thus in the following we assume that $\sigma_2 - a\sigma_1\rho \neq 0$. So $\tau^-$ is the first passage time of $W'$ through the horizontal axis and $\tau^+$ is the time of the first passage of $W'$ through the line $y = a'x$, where

$$a' = \frac{a\sigma_1\sqrt{1 - \rho^2}}{\sigma_2 - a\sigma_1\rho}, \qquad (3.2)$$

so that $\tau = \tau^- \wedge \tau^+$ is the exit time of $W'$ from the wedge $\mathcal{D}' = \langle \alpha' \rangle$, where $\alpha'$ is defined as $\arctan(a')$ or $\pi + \arctan(a')$ or $2\pi + \arctan(a)$ depending on the cases. The process $W'$ starts at $x'_0 = \sigma^{-1} \cdot x_0$. The possible cases are summed up in Table 3.2.

After the above transformation the original wedge is transformed into a new one according to this table. Therefore mutatis-mutandis, we shall omit the primes in the notation and without loss of generality we directly assume that the process $W$ has the correlation matrix $\Sigma = I_2$ and the wedge is given as in (3.1) with angle $\alpha$.

In some results, we assume that the angle of the wedge satisfies the condition:

$$\alpha = \frac{\pi}{m},$$

for some $m \in \mathbb{N}$. We define adjacent wedges as follows. For $k = 0, 1, \dots, 2m-1$, let

$$\mathcal{D}_k := \{(r \cos(\theta), r \sin(\theta)) : r \in [0, +\infty), \theta \in [k\alpha, (k+1)\alpha]\}.$$

We remark that $\mathcal{D}_0 = \mathcal{D}$ and that $\mathbb{R}^2 = \bigcup_{k=0}^{2m-1} \mathcal{D}_k$. We also define the transformation $T_k$:

$$T_k : \mathcal{D}_0 \to \mathcal{D}_k \\ T_k((r \cos \theta, r \sin \theta)) := (r \cos(\vartheta_k), r \sin(\vartheta_k))$$

$$\vartheta_k := \begin{cases} (k+1)\alpha - \theta; & k \text{ odd,} \\ k\alpha + \theta; & k \text{ even.} \end{cases}$$

From the definition, it follows that $T_k$ is an isometric bijection between $\mathcal{D}_0$ and $\mathcal{D}_k$. This setting is represented in Figure 3.3.

The density of the killed Brownian motion and some related random quantities appear in [Met10]. We first quote the explicit density result for killed Brownian motion which is fully proven in [BnS97, Lemma 1] for generalized multidimensional cones.
<figure><img src="image_3.png" /><figcaption>Figure 3.3: Partition of $\mathbb{R}^2$ using $\{\mathcal{D}_k\}_k$</figcaption></figure>

**Theorem 1 (Killed case).** Assume that $\mathcal{D} = \langle \pi/m \rangle$ for some $m \in \mathbb{N}$. Then we have the following formula for the density of killed Brownian motion on the wedge $\mathcal{D}$:

$$ \forall t > 0, x, y \in \mathcal{D}, \mathbb{P}^x(W_t \in dy, \tau > t) = \frac{1}{2\pi t} \sum_{k=0}^{2m-1} (-1)^k e^{-\frac{|x-T_k y|^2}{2t}} dy. \quad (3.3) $$

In general, let $\alpha \in (0, 2\pi)$, $\mathcal{D} = \langle \alpha \rangle$ and $x = (r_0 \cos(\theta_0), r_0 \sin(\theta_0)) \in \mathcal{D}$. Then for any $y = (r \cos(\theta), r \sin(\theta)) \in \mathcal{D}$, the density of the killed Brownian motion on the wedge $\mathcal{D}$ is given by

$$ \mathbb{P}^x(W_t \in dy, \tau > t) = \frac{2r}{t\alpha} e^{-\frac{r^2+r_0^2}{2t}} \sum_{n=1}^{\infty} I_{n\pi/\alpha} \left( \frac{rr_0}{t} \right) \sin \left( \frac{n\pi\theta}{\alpha} \right) \sin \left( \frac{n\pi\theta_0}{\alpha} \right) drd\theta. \quad (3.4) $$

# 4 Analytic formulas for the density of the reflected process

In this section, we give explicit expressions for the densities of the reflected Brownian motion. Before doing this, we review the one dimensional case which will also help explain the motivation for the simulation methods to be introduced later. The one-dimensional case can be treated using the classic reflection principle, which directly gives the density of the killed and reflected Brownian motions. In the one-dimensional case, we can assume that $\mathcal{D} = [0, +\infty)$ and that $\Sigma = 1$. Then for all non-negative measurable functions $f: \mathcal{D} \to \mathbb{R}$ and $x_0 \ge 0$,

$$ \mathbb{E}^{x_0} [f(W_T) \mathbf{1}_{\tau>T}] = \mathbb{E}^{x_0} \left[ f(W_T) \mathbf{1}_{W_T>0} \left( 1 - e^{-\frac{x_0(x_0+W_T)}{T}} \right) \right], \quad (4.1) $$

$$ \mathbb{E}^{x_0} [f(X_T)] = \mathbb{E}^{x_0} \left[ f(W_T) \mathbf{1}_{W_T>0} \left( 1 + e^{-\frac{x_0(x_0+W_T)}{T}} \right) \right]. \quad (4.2) $$

The above formulas show that the simulation of the killed and reflected Brownian motions in the one-dimensional case can be easily performed using changes of measures on the original unrestricted Brownian motion.

The proof in the reflected case in $\mathbb{R}^2$ shares common arguments with the killed case. Throughout these arguments, one uses results from classical theory of partial differential equations on wedges. For a general reference on this topic, we refer the reader to [LSU68, chapter IV] and [KZ16] where one can find existence and uniqueness results for the partial differential equations that are related to the problems in this article.

**Definition 2.** We define the normal reflection process of a continuous adapted stochastic process $\xi$ on the subset $\mathcal{D}$, as the unique solution $(X_t, L_t)_{t \ge 0}$ of the following problem (see for example [Pil14]):

1. $X_t = \xi_t + L_t$, where $L = (L_t)_{t \ge 0}$ is some adapted process.

2. $\forall t \ge 0$, $X_t \in \mathcal{D}$.

3. $L$ is a continuous process of bounded variation and $L_0 = 0$.
4. $\forall t \ge 0$, $L_t = \int_0^t n(X_s)1_{X_s \in \partial D} d|L|_s$, which means that the process $L$ increases only when $X_t \in \partial D$ and that the reflection is normal to the boundary.

Therefore the normally reflected Brownian motion is the process that is obtained from the above equation by taking $\xi = W$. In the one dimensional case, $L$ corresponds to the local time of Brownian motion.

The existence and uniqueness of a Brownian motion reflected in a wedge, which is a non-smooth domain, has been proved in the decorrelated case in [VW85]. In this article we will only consider normal reflections in the case of unitary diffusion matrix. Note that if we use the decorrelation step as in (3.2) for the reflected Brownian motion, this transformation changes the normal reflected process into an obliquely reflected process. We do not know how to obtain the density of such a process.

**Theorem 3 (Reflected case).** Assume that $D = \langle \pi/m \rangle$ for some $m \in \mathbb{N}$. Then we have the following formulas for the density of the reflected Brownian motion on the wedge $D = \langle \pi/m \rangle$:

$$ \forall t > 0, x, y \in D, \mathbb{P}^x(X_t \in dy) = \frac{1}{2\pi t} \sum_{k=0}^{2m-1} e^{-\frac{|x-T_k y|^2}{2t}} dy. \quad (4.3) $$

In the general case for $\alpha \in (0, 2\pi)$, $D = \langle \alpha \rangle$ and $x = (r_0 \cos(\theta_0), r_0 \sin(\theta_0)) \in D$ we have that for any $y = (r \cos(\theta), r \sin(\theta)) \in D$, the density of the reflected Brownian motion on the wedge $D = \langle \alpha \rangle$ is given by

$$ \mathbb{P}^x(X_t \in dy) = \frac{2r}{t\alpha} e^{-\frac{r^2+r_0^2}{2t}} \left( \frac{1}{2} I_0\left(\frac{rr_0}{t}\right) + \sum_{n=1}^{\infty} I_{n\pi/\alpha}\left(\frac{rr_0}{t}\right) \cos\left(\frac{n\pi\theta}{\alpha}\right) \cos\left(\frac{n\pi\theta_0}{\alpha}\right) \right) drd\theta. \quad (4.4) $$

The proof of the above theorem is given in B. The idea of using the simple case $D = \langle \pi/m \rangle$ to infer a general result appears in the literature of the image method. It has been used in a non-trivial way in order to deduce stationary measures for reflected processes (see e.g. [DM09]). A formula similar to (4.4) appears in [CJ59, page 379 (8)], which is proved using Laplace transform inversion methods. We have decided to include a proof which uses the image method because it gives some intuition on the simulations methods to be introduced later.

# 5 Exact simulation algorithms

Theorems 1 and 3 give formulas that can be directly used to simulate the final values of the reflected and stopped processes in the case that the wedge angle is $\alpha = \pi/m$ (see Sections 5.2 and 5.3). However, direct algorithms that arise due to these results demand the use of approximations and larger computational time in general. Instead, we investigate an alternative simulation method using the tractable case $\alpha = \pi/m$ and give algorithms for the exact simulation of the reflected and stopped processes in any wedge.

## 5.1 Formulas for the simulation of $(\tau, W_\tau)$ in the case $\alpha = \pi/m$

Before providing the simulation method, we give a formula that will be the key to provide a simplified simulation method which avoids the calculation of Bessel functions.

**Theorem 4.** Let $V = \langle \alpha^{-}, \alpha^{+} \rangle$ and assume that $\alpha^{+} - \alpha^{-} = \pi/m =: \alpha$ for some $m \in \mathbb{N}$ and let $\tau$ be the hitting time on the wedge $V$ for $W$ which starts at $x_0 = (r_0 \cos(\theta_0), r_0 \sin(\theta_0)) \in V$, then

$$ \mathbb{P}^{x_0}(\tau \in dt, W_{\tau} \in dy^{\pm}) = \frac{r_0}{2\pi t^2} e^{-\frac{r^2+r_0^2}{2t}} \sum_{k=0}^{m-1} \sin(\gamma_k^{\pm}) e^{\frac{rr_0}{t}} \cos(\gamma_k^{\pm}) dr dt, \quad (5.1) $$

where $y^{\pm} = (r \cos(\alpha^{\pm}), r \sin(\alpha^{\pm})) \in \partial V^{\pm}$. We use the notation

$$ \gamma_k^{\pm} := \pm \alpha^{\pm} \pm \frac{2k\pi}{m} \mp \theta_0. \quad (5.2) $$
Proof. We can assume that $\alpha^{-} = 0$ without loss of generality. We use the formulas [Met10, (1.5), (1.6)]. If $y \in \partial D^{-}$, writing $y = (r, 0)$ for some $r > 0$, we have:

$$ \mathbb{P}^{x_0}(\tau \in dt, W_{\tau} \in dy) = \frac{\pi}{\alpha^2 tr} e^{-\frac{r^2+r_0^2}{2t}} \sum_{n=1}^{\infty} n \sin\left(\frac{n\pi\theta_0}{\alpha}\right) I_{n\pi/\alpha}\left(\frac{rr_0}{t}\right) drdt. \quad (5.3) $$

Using (3.4) and switching the order of derivation and integration, which can be justified (see the proof of Theorem 3 with Corollary 16 in the Appendix), we obtain for $y \in \partial D^{-}$:

$$ \mathbb{P}^{x_0}(\tau \in dt, W_{\tau} \in dy) = \frac{1}{2r^2} \left. \frac{\partial}{\partial \theta} \right|_{\theta=0} \mathbb{P}^{x_0}(\tau > t, W_t \in d(r \cos(\theta), r \sin(\theta))). $$

So that in the case $\alpha = \pi/m$, using (3.3) in polar coordinates and the definition of $\vartheta_k$:

$$ \begin{aligned} \mathbb{P}^{x_0}(\tau \in dt, W_{\tau} \in dy) &= \frac{1}{2r^2} \left. \frac{\partial}{\partial \theta} \right|_{\theta=0} \frac{r}{2\pi t} e^{-\frac{r^2+r_0^2}{2t}} \sum_{k=0}^{2m-1} (-1)^k e^{\frac{1}{t}rr_0 \cos(\theta_0 - \vartheta_k)} dr dt \\ &= \frac{r_0}{2\pi t^2} e^{-\frac{r^2+r_0^2}{2t}} \sum_{k=0}^{m-1} \sin(\theta_0 - 2k\alpha) e^{\frac{rr_0}{t} \cos(\theta_0 - 2k\alpha)} dr dt. \end{aligned} \quad (5.4) $$

Likewise, we obtain a similar formula for the case $y \in \partial D^{+}$.

From here, we present a method to simulate $\tau$ and $W_{\tau}$ in the case of the wedge $\mathcal{V} = \langle \alpha^{-}, \alpha^{+} \rangle$ with $\alpha^{+} - \alpha^{-} = \pi/m$. For the simulation of $W_{\tau}$, we use the same method which is proposed in [Met09, Section 2.1] but we restrict to the case $\alpha = \pi/m$. On the other hand, we propose an unbiased simulation method for $\tau$ in comparison with the method proposed in [Met09] which has a bias difficult to estimate (see [Met09, Proposition 2.1.8] and the remark that follows).

**Simulation of $W_{\tau}$:** First, we simulate the selection of the boundary line on which $W_{\tau}$ arrives, with a Bernoulli distribution. These two boundaries are denoted by $\partial \mathcal{V}^{\pm} = \{(r \cos(\alpha^{\pm}), r \sin(\alpha^{\pm})) : r \ge 0\}$. Following [Met09, Corollary 2.1.5] we have:

$$ \mathbb{P}^{x_0}(W_{\tau} \in \mathcal{V}^{+}) = \frac{\theta_0 - \alpha^{-}}{\alpha^{+} - \alpha^{-}}, \quad (5.5) $$

We simulate on which frontier $W_{\tau}$ arrives first and then we directly simulate $W_{\tau}$, without simulating $\tau$ yet. Following [Met09, Proposition 2.1.6] one can compute explicit formulas for the distribution function of the radius of exit given that it exits at $\mathcal{V}^{+}$ or $\mathcal{V}^{-}$. Also its inverse can be computed and therefore the inverse transformation method for simulation can be applied. Using these formulas, we can simulate the radius $r_{\tau}$ as:

$$ r_{\tau} = \begin{cases} r_0 \left( \cos\left(\frac{\pi\theta_0}{\alpha}\right) - \frac{\sin\left(\frac{\pi\theta_0}{\alpha}\right)}{\tan\left((\pi - \pi\theta_0/\alpha)(U-1)\right)} \right)^{\alpha/\pi} & \text{if } W_{\tau} \in \mathcal{V}^{-}, \\ r_0 \left( -\cos\left(\frac{\pi\theta_0}{\alpha}\right) - \frac{\sin\left(\frac{\pi\theta_0}{\alpha}\right)}{\tan\left((\pi\theta_0/\alpha)(U-1)\right)} \right)^{\alpha/\pi} & \text{if } W_{\tau} \in \mathcal{V}^{+}, \end{cases} \quad (5.6) $$

with $U \sim \mathcal{U}([0, 1])$.

**Simulation of $\tau$:** Knowing $W_{\tau}$, we simulate $\tau$ according to the conditional formula (see (5.4) and recall that $dy^{\pm} = (dr \cos(\alpha^{\pm}), dr \sin(\alpha^{\pm}))$):

$$ \mathbb{P}^{x_0}(\tau \in dt | W_{\tau} = y^{\pm}) = ((\mathbb{P}^{x_0}(W_{\tau} \in dy^{\pm}))^{-1} dr) \frac{r_0}{2\pi t^2} \sum_{k=0}^{m-1} \sin(\gamma_k^{\pm}) \exp \left( - \frac{(r - r_0 \cos(\gamma_k^{\pm}))^2 + r_0^2 \sin^2(\gamma_k^{\pm})}{2t} \right) dt. \quad (5.7) $$

Since some of the coefficients $\sin(\gamma_k^{\pm})$ are negative, we cannot directly simulate according to this distribution as a mixture. Instead we select the indexes $k_0^{\pm}, \dots, k_{p^{\pm}}^{\pm}$ such that for all $i$, $\sin(\gamma_{k_i^{\pm}}^{\pm}) \ge 0$ and write:

$$ \mathbb{P}^{x_0}(\tau \in dt | W_{\tau} = y^{\pm}) \le ((\mathbb{P}^{x_0}(W_{\tau} \in dy^{\pm}))^{-1} dr) \frac{r_0}{2\pi t^2} \sum_{i=0}^{p^{\pm}} \sin(\gamma_{k_i^{\pm}}^{\pm}) \exp \left( - \frac{(r - r_0 \cos(\gamma_{k_i^{\pm}}^{\pm}))^2 + r_0^2 \sin^2(\gamma_{k_i^{\pm}}^{\pm})}{2t} \right) dt. \quad (5.8) $$
This is the distribution of a discrete mixture of random variables $\xi_k$ for $0 \le k \le p^{\pm}$, where $\xi_k$ is the inverse of an exponential random variable of parameter $((r - r_0 \cos(\gamma_{k_i^{\pm}}))^{2} + r_0^{2} \sin^{2}(\gamma_{k_i^{\pm}}))/2$. With this information, we build an acceptance-rejection sampling, based on the density on the right hand side of (5.8) and accept the sample value $t$ under the condition

$$ U < \left( \sum_{k=0}^{m-1} \sin(\gamma_k^{\pm}) \exp \left( - \frac{(r - r_0 \cos(\gamma_k^{\pm}))^2 + r_0^2 \sin^2(\gamma_k^{\pm})}{2t} \right) \right) \times \left( \sum_{i=0}^{p^{\pm}} \sin(\gamma_{k_i^{\pm}}^{\pm}) \exp \left( - \frac{(r - r_0 \cos(\gamma_{k_i^{\pm}}^{\pm}))^2 + r_0^2 \sin^2(\gamma_{k_i^{\pm}}^{\pm})}{2t} \right) \right)^{-1}, $$

with $U \sim U([0, 1])$.

**Simulation in the case W does not exit the wedge:** So far, we have provided a methodology in order to simulate the exit location and time for the simpler wedge $V = \langle \alpha^{-}, \alpha^{+} \rangle$ with $\alpha^{+} - \alpha^{-} = \pi/m$. Now we discuss the case when W does not leave this wedge in the time interval $[0, T]$. To simulate the final value $W_T$ of the Brownian motion starting from $x_0$ conditionally to the fact that it does not exit the wedge $V$ on $[0, T]$, we need to simulate according to the density (3.3), which is a sum where some terms are negative. We propose the following method.

A slight generalization of (3.3) to the wedge $V = \langle \alpha^{-}, \alpha^{+} \rangle$ with

$$ \tilde{\vartheta}_k := \begin{cases} \vartheta_k + 2\alpha^{-} & \text{if } k \text{ odd} \\ \vartheta_k & \text{if } k \text{ even,} \end{cases} $$

gives, for $y \in V$:

$$ \mathbb{P}^{x_0}(W_T \in dy | \tau > T) = \frac{\mathbb{P}^{x_0}(\tau > T)^{-1}}{2\pi T} \sum_{k=0}^{2m-1} (-1)^k e^{-\frac{|x_0 - \tilde{T}_k y|^2}{2T}} dy \le \frac{\mathbb{P}^{x_0}(\tau > T)^{-1}}{2\pi T} \sum_{k=0}^{m-1} e^{-\frac{|x_0 - \tilde{T}_{2k} y|^2}{2T}} dy, \quad (5.9) $$

where $\tilde{T}_k(r \cos(\theta), r \sin(\theta)) = (r \cos(\tilde{\vartheta}_k), r \sin(\tilde{\vartheta}_k))$. We then apply acceptance-rejection sampling with the reference density proportional to $\frac{1}{2\pi T} \sum_{k=0}^{m-1} e^{-\frac{|x_0 - \tilde{T}_{2k} y|^2}{2T}}$ and accept the sample $y$ under the condition

$$ U < \left( \sum_{k=0}^{2m-1} (-1)^k e^{-\frac{|x_0 - \tilde{T}_k y|^2}{2T}} \mathbb{1}_{y \in V} dy \right) \left( \sum_{k=0}^{m-1} e^{-\frac{|x_0 - \tilde{T}_{2k} y|^2}{2T}} \right)^{-1}, \quad (5.10) $$

with $U \sim U([0, 1])$.

## 5.2 Algorithm for the simulation of the stopped Brownian motion: General case

In this subsection, we give a recursive algorithm to simulate the final value $W_{T \wedge \tau}$ of the stopped process on the wedge $D$, for any angle $\alpha \in (0, 2\pi)$. This algorithm will be used for the simulation algorithm of the reflected Brownian motion, but it also has its own interest as it can be used as itself for the simulation of the stopped Brownian motion. In what follows, we denote $F_n$ the sigma algebra generated by the algorithm until the step $n$, i.e. $F_n = \sigma(\tau_1, y_1, \dots, \tau_n, y_n)$, where $\tau_1, \dots, \tau_n$ and $y_1, \dots, y_n$ are random variables generated by the algorithm. We also denote by $P_n$ and $E_n$ the conditional probability and expectation with respect to $F_n$.

First, we propose an algorithm for the case where the angle is $\pi/m$. At the step $n$ of the algorithm we will define a wedge $V_n \subset D$ of angle $\pi/m$ for some $m \in \mathbb{N}$ and then simulate the process stopped on this wedge. Define $\Theta := \max(\{\pi/m, m \in \mathbb{N}\} \cap (0, \alpha])$. Note that since $\alpha \in (0, 2\pi)$, we have $\alpha/2 < \Theta \le \alpha$.

**Algorithm I:** Start at the point $x_0 = (r_0 \cos(\theta_0), r_0 \sin(\theta_0)) =: y_0$ at time $T_0 = 0$. The algorithm follows these steps for $n \in \mathbb{N}$:

1. Given $y_n = (r_n \cos(\theta_n), r_n \sin(\theta_n)) \in V_n$, define the angles:

$$ \beta_{n+1}^{-} := \begin{cases} 0 & \text{if } \theta_n \in [0, \frac{\Theta}{2}] \\ \theta_n - \frac{\Theta}{2} & \text{if } \theta_n \in [\frac{\Theta}{2}, \alpha - \frac{\Theta}{2}] \\ \alpha - \Theta & \text{if } \theta_n \in [\alpha - \frac{\Theta}{2}, \alpha] \end{cases} $$
<figure><img src="image_4.png" /><figcaption>Figure 5.1: Example of the domains of simulation at step $n + 1$.</figcaption></figure>

$$ \beta_{n+1}^{+} := \beta_{n+1}^{-} + \Theta. $$

By the definition of $\Theta$, we have $\Theta \le \alpha$ and therefore $0 \le \Theta/2 \le \alpha - \Theta/2 \le \alpha$.

2. Consider the wedge $V_{n+1} = \langle \beta_{n+1}^{-}, \beta_{n+1}^{+} \rangle$ of angle $\Theta$, which satisfies $V_{n+1} \subset \mathcal{D}$ and $y_n = (r_n \cos(\theta_n), r_n \sin(\theta_n)) \in V_{n+1}$.

3. Simulate on which one of the two sets $\partial V_{n+1}^{+}$ or $\partial V_{n+1}^{-}$, the process $(W_t)_{t \ge T_n}$ starting from $y_n$, at time $T_n$, exits the wedge $V_{n+1}$ and define $\theta_{n+1} := \beta_{n+1}^{\pm}$ according to the simulation result. This is done using a Bernoulli variable and the formula (5.5).

4. Simulate $y_{n+1} \in \partial V_{n+1}$, which is the value of $W$ when it reaches $\partial V_{n+1}$ for the first time after starting at $y_n$ at time $T_n$, and knowing which one of the two events $y_{n+1} \in \partial V_{n+1}^{\pm}$ has occurred, using formula (5.6). Denote the exit point as

$$ y_{n+1} = (r_{n+1} \cos(\theta_{n+1}), r_{n+1} \sin(\theta_{n+1})). $$

5. Simulate $\tau_{n+1}$, the time for the process $W$ starting at $y_n = (r_n \cos(\theta_n), r_n \sin(\theta_n))$ at time $T_n$ to reach $\partial V_{n+1}$, knowing the event $W_{T_n+\tau_{n+1}} = y_{n+1}$. This can be done using the formula (5.7) and the acceptance-rejection procedure described after it. Define

$$ T_{n+1} := T_n + \tau_{n+1}, $$

so that $W_{T_{n+1}} = y_{n+1}$.

If $T_n + \tau_{n+1} < T$ then the algorithm iterates. This algorithm stops under one of these two conditions:

* Condition 1: $T_n + \tau_{n+1} \ge T$. Then simulate the final value $W_T$ knowing that $W_{T_n} = y_n$ and that for $t \in [T_n, T]$, $W_t \in V_n$. We do this simulation conditionally to the event $T_n + \tau_{n+1} > T$ and "forget" the exact simulated value of $\tau_{n+1}$ as well as the value of $W_{T_n+\tau_{n+1}}$. This is justified in Proposition 11, using $X := (W_{T_n+\tau_{n+1}}, T_{n+1} \mathbf{1}_{\{T_{n+1} < T\}})_{\mathcal{F}_n}$, $Y := (W_{T_n+\tau_{n+1}}, T_{n+1})_{\mathcal{F}_n}$ and $A = \mathbb{R}^2 \times [0, T]$. We simulate the value of $W_T$ using (5.9) and the acceptance-rejection procedure described after it.

* Condition 2: $\theta_{n+1} = 0$ or $\theta_{n+1} = \alpha$. Then we obtain $\tau := T_{n+1} < T$, i.e. the Brownian motion reaches the wedge $\mathcal{D}$ before the time $T$ and the result of the simulation is

$$ W_{T_n+\tau} = (r_{n+1} \cos(\theta_{n+1}), r_{n+1} \sin(\theta_{n+1})). $$

An illustration of this algorithm is given in Figure 5.1. Next, we prove that this algorithm can be carried out in a finite number of steps.

**Proposition 5.** Algorithm I ends in finite time, i.e. the number of required iterations to finish the algorithm, denoted by $N$, is such that $N < \infty$ almost surely. More precisely, for all $K \in \mathbb{N}$, $\mathbb{P}^{x_0}(N \ge K) \le 2^{-\lfloor K/2 \rfloor}$.
Proof. Fix $K \in \mathbb{N}$. For a wedge $\mathcal{V}_n$, $n \le N$ in Algorithm I, we say that "$\mathcal{V}_n$ intersects the boundary of $\mathcal{D}$" if $\beta_n^+ = \alpha$ or $\beta_n^- = 0$. For $n \le N$, we have two possible cases:

* $\mathcal{V}_n$ does not intersect the boundary of $\mathcal{D}$. Then, using the fact that $\Theta > \frac{\alpha}{2}$ and the definition of $\beta_{n+1}^{\pm}$, $\mathcal{V}_{n+1}$ intersects the boundary of $\mathcal{D}$, and $\theta_n$ is closer to $\beta_{n+1}^{\pm}$ than $\beta_{n+1}^{\mp}$, where here the sign $\pm$ is such that $\partial\mathcal{V}_{n+1}^{\pm} \cap \partial\mathcal{D} \ne \{0\}$. So using (5.5), $\mathbb{P}_n(y_{n+1} \notin \partial\mathcal{D}) \le 1/2$.

* $\mathcal{V}_n$ intersects the boundary of $\mathcal{D}$. Then if $\theta_{n+1} \notin \{0, \alpha\}$, $\mathcal{V}_{n+1}$ does not intersect the boundary of $\mathcal{D}$ and by the same reasoning as above $\mathbb{P}_n(y_{n+2} \notin \partial\mathcal{D}) \le 1/2$.

Then, by tower property of conditional expectations, we have:

$$ \mathbb{P}^{y_0}(N \ge K) \le \mathbb{P}^{y_0}(y_1 \notin \partial\mathcal{D}, \dots, y_K \notin \partial\mathcal{D}) \le 2^{-\lfloor \frac{K}{2} \rfloor}. $$

So that:

$$ \mathbb{P}^{y_0}(N = \infty) = \lim_{K \to \infty} \mathbb{P}^{y_0}(N \ge K) = 0. $$

□

## 5.3 Algorithm for the simulation of the reflected Brownian motion

We now state an algorithm to simulate the final value of the reflected Brownian motion $X_T$ in a wedge $\mathcal{D} = \langle \alpha \rangle$ of any angle $\alpha \in (0, 2\pi)$.

**Algorithm II:** Choose $\Theta := \max(\{\pi/m, m \in \mathbb{N}\} \cap \{0, \alpha\})$ as before. Start at the point $x_0 = (r_0 \cos(\theta_0), r_0 \sin(\theta_0)) =: y_0$ at time $T_0 = 0$. In general, suppose that for $n \in \mathbb{N}$ the point $y_n = (r_n \cos(\theta_n), r_n \sin(\theta_n))$ has already been simulated.

1. Define the angles:

$$ \beta_{n+1}^{-} := \theta_{n} - \frac{\Theta}{2}, $$

$$ \beta_{n+1}^{+} := \theta_{n} + \frac{\Theta}{2} = \beta_{n+1}^{-} + \Theta. $$

2. Consider the wedge $\mathcal{V}_{n+1} := \langle \beta_{n+1}^{-}, \beta_{n+1}^{+} \rangle$. Note that although $y_n \in \mathcal{V}_{n+1}$, we do not necessarily have that $\mathcal{V}_{n+1} \subset \mathcal{D}$.

3. Simulate the random variables $\tau_{n+1}$ and $z_{n+1}$, which are respectively the time and the final point of a Brownian motion $Z = (Z_t)_{t \ge 0}$ starting at $y_n$ at time 0 and stopped on the wedge $\mathcal{V}_{n+1}$ using the steps 3-5 of Algorithm I in Section 5.2. Note that, since $\theta_n$ is in the middle of the wedge $\mathcal{V}_{n+1}$, we have by symmetry $\mathbb{P}^{z_n}(Z_{\tau_{n+1}} \in \partial\mathcal{V}_{n+1}^{\pm}) = 1/2$.

4. If $T_n + \tau_{n+1} < T$, then we define $r_{n+1}$ and $\tilde{\theta}_{n+1}$ to be respectively the radius and the angle of $z_{n+1}$. Since it is possible that $z_{n+1} \notin \mathcal{D}$, we define $\theta_{n+1}$ as follows:

$$ \theta_{n+1} = \begin{cases} \tilde{\theta}_{n+1} & \text{if } \tilde{\theta}_{n+1} \in [0, \alpha], \\ -\tilde{\theta}_{n+1} & \text{if } \tilde{\theta}_{n+1} = \beta_{n+1}^{-} \text{ and } \beta_{n+1}^{-} < 0, \\ 2\alpha - \tilde{\theta}_{n+1} & \text{if } \tilde{\theta}_{n+1} = \beta_{n+1}^{+} \text{ and } \beta_{n+1}^{+} > \alpha. \end{cases} $$

And then we define $y_{n+1} := (r_{n+1} \cos(\theta_{n+1}), r_{n+1} \sin(\theta_{n+1}))$, so that $y_{n+1} \in \mathcal{D}$. Note that if $z_{n+1} \notin \mathcal{D}$, then $y_{n+1}$ is the reflection of $z_{n+1}$ with respect to the line $\{\theta = 0\}$ in the second case or $\{\theta = \alpha\}$ in the third case in the definition of $\theta_{n+1}$. Then define the time:

$$ T_{n+1} = T_n + \tau_{n+1} $$

and the algorithm iterates.
<figure><img src="image_5.png" /><figcaption>Figure 5.2: Example of domains of simulation at the step $n + 1$ in the reflected case.</figcaption></figure>

5. If $T_n + \tau_{n+1} > T$, we simulate $z_{n+1} \in \mathcal{V}_{n+1}$ as the value at time $T - T_n$ of a Brownian motion starting at $y_n$ and conditionally to the fact that it stays in the wedge $\mathcal{V}_{n+1}$ in the time interval $[0, T - T_n]$. This can be done using the acceptance-rejection method given in (5.9) and what follows based on Proposition 11 as it was done in the termination Condition 1 in Algorithm I in Section 5.2. Then we define $y_{n+1}$ in function of $z_{n+1}$ as in step 4. Finally, we stop the algorithm, and the resulting value of the simulation is $y_{n+1}$.

An illustration of this algorithm is given in Figure 5.2.

**Proposition 6.** Define $N := \inf\{n \in \mathbb{N}, \tau_1 + \dots + \tau_n > T\}$. Then $N < \infty$ almost surely. That is, the algorithm terminates in finite time.

*Proof.* The stopping times $(\tau_i)_i$ follow the same law as $(\tilde{\tau}_i)$, defined by $\tilde{\tau}_0 = 0$ and:

$$ \tilde{\tau}_{i+1} := \inf \left\{ t > 0 : \left| \theta(B_{t+\tilde{\tau}_i}) - \theta(B_{\tilde{\tau}_i}) \right| \ge \frac{\pi}{2m} \right\}, \quad (5.11) $$

where $m = \pi/\Theta \in \mathbb{N}$, $B = (B_t)_{t \ge 0}$ is a standard two-dimensional Brownian motion and $\theta(B_t)$ denotes the angle of the process<fnref n="2" /> B. Then:

$$ \mathbb{P}^{x_0} \left( \sum_{i=1}^{\infty} \tau_i < \infty \right) = \mathbb{P}^{x_0} \left( \sum_{i=1}^{\infty} \tilde{\tau}_i < \infty \right). $$

But, if $\mathcal{T} := \sum_{i=1}^{\infty} \tilde{\tau}_i < \infty$, then there exist two random sequences, $(t_n^1)$ and $(t_n^2)$, increasing and converging to $\mathcal{T}$, such that $\theta(B_{t_n^1}) = k_n^i \pi/(2m) + \theta_0$ for all $n \in \mathbb{N}$ and $i \in \{1, 2\}$, and where $k_n^i \in \mathbb{Z}$ has the same parity as $i$, so that the difference between $\theta(B_{t_n^1})$ and $\theta(B_{t_n^2})$ is at least $\pi/(2m)$ for all $n, n' \in \mathbb{N}$. Taking $n \to \infty$, we have necessarily that $B_{\mathcal{T}} = 0$, which occurs with probability 0 on any closed time interval. $\square$

# 6 Folding number in a wedge and complexity

In this section we study the complexity of Algorithms I and II in Sections 5.2 and 5.3 in separate cases. Recall that $N$ denotes the number of iterations that an algorithm requires to finish one simulation, i.e.

$$ N := \inf\{n \in \mathbb{N} : \tau_1 + \dots + \tau_n > T\}. $$

For Algorithm I, $N$ is the number of times we touch the boundaries of the corresponding sets $\mathcal{V}_n$, $n \in \mathbb{N}$, before touching one of the boundaries of $\mathcal{D}$ or reaching time $T$. For Algorithm II, $N$ is the number of times that a Brownian motion reflected in a wedge $\mathcal{V}$ of angle $\pi/(2m)$ goes from one boundary $\partial\mathcal{V}^{\pm}$ to another $\partial\mathcal{V}^{\mp}$ until time $T$. For this reason, we may call $N$ the number of folds for the simulation of the reflected Brownian motion.

In this section, we give theoretical properties and bounds for $N$ for each algorithm separately.

<footnote n="2">²We choose the angle of the process $(\theta(B_t))_{t \ge 0}$ so that it is continuous with respect to $t$.</footnote>
## 6.1 Majoration of N for the simulation algorithm of the stopped Brownian motion

For Algorithm I, by Proposition 5, we have for all $K \in \mathbb{N}$, $\mathbb{P}^{x_0}(N \ge K) \le 2^{-\lfloor K/2 \rfloor}$, so

$$ \mathbb{E}^{x_0}[N] = \sum_{K=1}^{\infty} \mathbb{P}^{x_0}(N \ge K) \le \sum_{K=1}^{\infty} 2^{-\lfloor \frac{K}{2} \rfloor} = 3. $$

So Algorithm I ends in finite time almost surely and its complexity is finite in expectation.

## 6.2 Majoration of the number of folds of the Brownian motion in a wedge

We now investigate the complexity for Algorithm II. We decompose the standard two-dimensional Brownian motion $W = (W_t)_{t \ge 0}$ in polar coordinates and introduce the notations that will be used in this section and the next one as follows. The successive stopping times $(\tau_i)_i$ which appear in Algorithm II, have the same law as $(\tilde{\tau}_i)_i$, defined as in the proof of Proposition 6. Then we have

$$ N \stackrel{\mathcal{L}}{=} \inf\{n \in \mathbb{N} : \tilde{\tau}_1 + \dots + \tilde{\tau}_n > T\}. $$

Moreover, using the skew-product representation of the Brownian motion (see [RY99, page 194]):

$$ B_t = R_t U_{F(t)}, \qquad (6.1) $$

where $(R_t)_{t \ge 0} = (|B_t|)_{t \ge 0}$ is a Bessel process of dimension 2, $(U_t)_{t \ge 0}$ is a Brownian motion on $\mathbb{S}^1 \subset \mathbb{R}^2$ and $F$ is defined as

$$ F(t) = \int_{0}^{t} \frac{ds}{R_{s}^{2}}, \qquad (6.2) $$

which is a strictly increasing process almost surely. Moreover, the processes $(R_t)_{t \ge 0}$ and $(U_t)_{t \ge 0}$ are independent. If we let $(\theta(U_t))_{t \ge 0}$ be the angle of the process $(U_t)_{t \ge 0}$ then $(\theta(U_t))_{t \ge 0}$ is a standard one-dimensional Brownian motion.

Define the stopping times $(s_i)_i$ by $s_0 = 0$ and

$$ s_{i+1} = \inf\{t > 0 : |\theta(U_{t+s_i}) - \theta(U_{s_i})| \ge \pi/(2m)\}. $$

Then, the random variables $s_i$ are i.i.d. and have the law of the time for one-dimensional Brownian motion starting at 0 to reach the double barrier $\pm \pi/(2m)$. Then using (5.11) we have that for all $K \in \mathbb{N}$,

$$ \sum_{i=1}^{K} s_i \stackrel{\mathcal{L}}{=} F \left( \sum_{i=1}^{K} \tilde{\tau}_i \right). $$

So the estimation of N can be simplified as

$$ \mathbb{P}^{x_0}(N \ge K) = \mathbb{P}^{x_0}\left(\sum_{i=1}^{K} \tilde{\tau}_i \le T\right) = \mathbb{P}^{x_0}\left(\sum_{i=1}^{K} s_i \le F(T)\right). \qquad (6.3) $$

Note that $F$ is a random change of time, independent of $(s_i)_i$. Furthermore, the Laplace transform of $s_i$ is explicitly given in [JYC09, Section 3.5, Proposition 3.5.1.3] as

$$ \mathbb{E}[e^{-\lambda s_1}] = \frac{2}{e^{\frac{\sqrt{2\pi}\sqrt{\lambda}}{2m}} + e^{-\frac{\sqrt{2\pi}\sqrt{\lambda}}{2m}}}. $$

In particular, we have:

$$ (\mathbb{E}[s_i], \mathrm{Var}(s_i)) = \left( \frac{\pi^2}{4m^2}, \frac{2}{3} \left( \frac{\pi}{2m} \right)^4 \right) =: (\mu, \sigma^2). $$

**Proposition 7.** Let $0 < a < \frac{1}{2}$. Then we have $\mathbb{E}^{x_0}[N^a] < \infty$. On the other hand, if $a \ge 1/2$ then $\mathbb{E}^{x_0}[N^a] = \infty$.
Equation (6.3) together with the estimate in (A.3) give the result which is a consequence of the fact that the tails of $F(T)$ are large because the probability that the Bessel process $R$ is close to zero is large and therefore due to (6.1) the simulation may spend a large time close to the origin.

*Proof.* We have that for all $K \in \mathbb{N}$ and $\delta > 1$:

$$ \mathbb{P}^{x_0}(N^a \geq K) = \int_{0}^{\frac{\mu K^{1/a}}{\delta}} \mathbb{P}^{x_0}\left(\sum_{i=1}^{\lfloor K^{1/a} \rfloor} s_i \leq y\right) \mathbb{P}^{x_0}(F(T) \in dy) + \int_{\frac{\mu K^{1/a}}{\delta}}^{\infty} \mathbb{P}^{x_0}\left(\sum_{i=1}^{\lfloor K^{1/a} \rfloor} s_i \leq y\right) \mathbb{P}^{x_0}(F(T) \in dy) \\ =: \mathcal{I}_1(K) + \mathcal{I}_2(K). $$

For $\mathcal{I}_1(K)$, we will use the fact that

$$ 2 \exp \left( \frac{\pi^2}{4\delta m^2} \right) < \exp \left( \frac{\sqrt{2\pi}}{2m} \right) + \exp \left( -\frac{\sqrt{2\pi}}{2m} \right) $$

for $\delta$ large enough. Then, we obtain that $\mathcal{I}_1(K)$ decreases exponentially fast to zero. In fact, using Markov inequality

$$ I_1(K) \le \mathbb{P}^{x_0}\left(\sum_{i=1}^{\lfloor K^{1/a} \rfloor} s_i \le \frac{\mu K^{1/a}}{\delta}\right) \le e^{\frac{\mu K^{1/a}}{\delta}} \left(\mathbb{E}^{x_0}[e^{-s_1}]\right)^{\lfloor K^{1/a} \rfloor} \le \left(\frac{2e^{\frac{\pi^2}{4\delta m^2}}}{e^{\frac{\sqrt{2\pi}}{2m}} + e^{-\frac{\sqrt{2\pi}}{2m}}}\right)^{\lfloor K^{1/a} \rfloor} . \quad (6.4) $$

Using (A.3), we have:

$$ I_2(K) \le \mathbb{P}^{x_0}\left(F(T) \ge \frac{\mu K^{1/a}}{\delta}\right) \underset{K \to \infty}{\sim} CK^{-\frac{1}{2a}}, \quad C = \frac{1}{\sqrt{2\Gamma(1/2)}} \left(\int_{\frac{r_0^2}{2T}}^{\infty} \frac{e^{-u}}{u} du\right) \left(\frac{\mu}{\delta}\right)^{-\frac{1}{2}}. \quad (6.5) $$

Since $0 < a < 1/2$, the result follows.

Similarly, in the case $a \ge 1/2$, one has that

$$ \mathbb{P}^{x_0}(N^a \geq K) \geq \int_{2\mu K^{1/a}}^{\infty} \mathbb{P}^{x_0}\left(\sum_{i=1}^{\lfloor K^{1/a} \rfloor} s_i \leq y\right) \mathbb{P}^{x_0}(F(T) \in dy). $$

For $K$ large enough and $y \ge 2\mu K^{1/a}$, we have $\mathbb{P}^{x_0}(\sum_{i=1}^{\lfloor K^{1/a} \rfloor} s_i \le y) \ge 1/2$ by the Central Limit Theorem, and $\mathbb{P}^{x_0}(F(T) \ge 2\mu K^{1/a}) \sim CK^{-1/(2a)}$ as $K \to \infty$ for some constant $C > 0$ which implies the result in this case. $\square$

From this result, we conclude that the average number of iterations of the algorithm in the reflected case is infinite. A modification of the above algorithm with finite number of iterations in expectation is provided below.

## 6.3 Proposition of modification of Algorithm II

When the simulated radii $r_1, \dots, r_n$ in Algorithm II of Section 5.3 become small, the process $(X_t)_{t \in [0,T]}$ goes from one boundary $\partial V^{\pm}$ to the other boundary $\partial V^{\mp}$ many times and the number of iterations becomes high. In Section C is hinted that this problem is in fact specific to the dimension 2. To remedy this problem, we study the asymptotic behavior of $\mathbb{P}^{x_0}(X_t \in dy)$ when $r_0$ is close to zero. In fact, if the ratio $rr_0/t$ is small, then the density of the reflected Brownian motion expressed in (4.4) can be approximated by the distribution obtained by taking only $n = 0$:

$$ \mathbb{P}^{x_0}(X_t \in dy) \approx \frac{r}{t\alpha} e^{-\frac{r^2+r_0^2}{2t}} I_0\left(\frac{rr_0}{t}\right) drd\theta. \quad (6.6) $$

The function above can be renormalized so that it becomes a density which can be simulated using inequality (A.1) and the acceptance-rejection method with the reference density

$$ \propto re^{-\frac{(r-r_0)^2}{2t}} drd\theta = \left( (r-r_0)e^{-\frac{(r-r_0)^2}{2t}} dr + r_0 e^{-\frac{(r-r_0)^2}{2t}} dr \right) d\theta. \quad (6.7) $$

Then we modify Algorithm II as follows. We choose some small $\varepsilon \in (0, r_0)$ and:
* Simulate the sequence of $(r_n)_n$ and $(\theta_n)_n$ as in Algorithm II in Section 5.3.

* If after some iteration $n$, we have

$$ \frac{r_n^2}{T - T_n} < \varepsilon \qquad (6.8) $$

(note that $T - T_n > 0$ is the remaining time of the simulated process after step $n$), then we directly simulate the final point $\bar{X}_T$ according to the approximation (6.6):

$$ \mathbb{P}^{x_0}(\bar{X}_T \in dy | X_{T_n} = y_n) = C(\varepsilon) \frac{r}{(T - T_n)\alpha} e^{-\frac{r^2 + r_n^2}{2(T - T_n)}} I_0\left(\frac{rr_n}{T - T_n}\right) drd\theta, \qquad (6.9) $$

$$ C(\varepsilon) := \left( \int_0^\alpha \int_0^\infty \frac{r}{(T - T_n)\alpha} e^{-\frac{r^2 + r_n^2}{2(T - T_n)}} I_0\left(\frac{rr_n}{T - T_n}\right) dr d\theta \right)^{-1}. $$

Using that $I_0(x) \ge 1$ for all $x \ge 0$, one obtains that $C(\varepsilon) \le e^{r_n^2/(2(T-T_n))}$. This upper bound for $C(\varepsilon)$ is enough in order to implement the acceptance-rejection sampling method. In this case, the algorithm directly ends after one additional iteration.

**Proposition 8.** Denote by $\bar{N}$ the number of iterations of the modified algorithm, (6.9), then for any $1 \le a < p < 2$, there exists a quantity $C(a, p, x_0, T, m) > 0$ such that for $\varepsilon$ small enough:

$$ \mathbb{E}^{x_0}[\bar{N}^a] \le \frac{C(a, p, x_0, T, m)}{\varepsilon^{p-1}}. \qquad (6.10) $$

*Proof.* We use a similar argument as in the proof of Proposition 7. Let $\zeta_\varepsilon := \inf\{t \in [0, T] : R_t^2/(T-t) \le \varepsilon\}$ and let $K \in \mathbb{N}$ fixed. Using Markov's inequality, we have:

$$ \begin{aligned} \mathbb{P}^{x_0}(\bar{N} - 1 \ge K^{1/a}) &\le \mathbb{P}^{x_0}\left(\sum_{i=1}^{\lfloor K^{1/a} \rfloor} s_i \le F(T \wedge \zeta_\varepsilon), F(T \wedge \zeta_\varepsilon) \le \frac{\mu K^{1/a}}{\delta}\right) + \mathbb{P}^{x_0}\left(F(T \wedge \zeta_\varepsilon) > \frac{\mu K^{1/a}}{\delta}\right) \\ &\le e^{\frac{\mu K^{1/a}}{\delta}} (\mathbb{E}^{x_0}[e^{-s_1}])^{\lfloor K^{1/a} \rfloor} + \frac{\delta^p}{(\mu K^{1/a})^p} \mathbb{E}^{x_0}[F(T \wedge \zeta_\varepsilon)^p]. \end{aligned} $$

From here the result follows by using Lemma 15. $\square$

**Proposition 9.** Denote by $\bar{X}_T$, the result obtained from the above approximation algorithm. Then, for $\alpha \in (0, 2\pi)$ and $\varepsilon$ small enough, the error made by the approximation in total variation distance satisfies:

$$ d_{\text{TV}}(\bar{X}_T, X_T) \le C_m \varepsilon^{\min(1, \frac{\pi}{2\alpha})}, $$

where $C_m$ is a constant which only depends on $m$.

*Proof.* If the approximation is not used, then $\bar{X}_T = X_T$. Else, we denote by $n$ the step when the approximation is used, so that we simulate the last step starting from $(r_n \cos(\theta_n), r_n \sin(\theta_n))$ with remaining time $T - T_n =: t'$. Since the $\varepsilon$-condition is reached, we have $r_n^2/t' < \varepsilon$. Then if we denote by $d_{\text{TV}_n}(\bar{X}_T, X_T)$ the total variation distance conditioned to the filtration $\mathbb{F}_n$ up to time $T_n$, and $d_n := 1/2$ for $n = 0$ and $d_n := 1$ otherwise, we have:

$$ \begin{aligned} d_{\text{TV}_n}(\bar{X}_T, X_T) &\le \int_{\mathcal{D}} |(C(\varepsilon) - 1 + 1) \frac{r}{t'\alpha} e^{-\frac{r^2+r_n^2}{2t'}} I_0\left(\frac{rr_n}{t'}\right) \\ &\quad - \frac{2r}{t'\alpha} e^{-\frac{r^2+r_n^2}{2t'}} \sum_{k=0}^{\infty} d_n I_{k\pi/\alpha}\left(\frac{rr_n}{t'}\right) \cos\left(\frac{k\pi\theta}{\alpha}\right) \cos\left(\frac{k\pi\theta_n}{\alpha}\right) |dr d\theta \\ &\le \frac{|C(\varepsilon) - 1|}{C(\varepsilon)} + \int_0^\infty \frac{2r}{t'} e^{-\frac{r^2+r_n^2}{2t'}} \sum_{k=1}^{\infty} I_{k\pi/\alpha}\left(\frac{rr_n}{t'}\right) dr. \end{aligned} $$
Using that for all $x \ge 0$, $I_0(x) \ge 1$, we have $C(\varepsilon) \le e^{r_n^2/(2t')}$. On the other hand, we have

$$C(\varepsilon)^{-1} \le \frac{1}{t'} \int_0^\infty re^{-r^2/(2t')} I_0\left(\frac{rr_n}{t'}\right) dr = \int_0^\infty ue^{-u^2/2} I_0\left(\frac{r_n u}{\sqrt{t'}}\right) du.$$

Using (2.1) one immediately obtains that $|I'_0(x)|+|I''_0(x)| \le 2e^x$. Therefore the function $A(\delta) := \int_0^\infty ue^{-u^2/2} I_0(\delta u) du$, $\delta \ge 0$ is twice differentiable. As $I_0(0) = 1$ and $I'_0(0) = 0$, we obtain $A(0) = 1$, $A'(0) = 0$ and then $A(\delta) = 1+O(\delta^2)$ as $\delta \to 0$ which gives $C(\varepsilon) \ge 1-C\varepsilon$ for small enough $\varepsilon > 0$. These bounds on $C(\varepsilon)$ imply that $|C(\varepsilon)-1|/(C(\varepsilon)\varepsilon)$ is bounded. For the second term, using the expression (2.1) and a multiplicity property of the number of terms of the type $k\pi/\alpha$ in the interval $[0, k]$, we have

$$\sum_{k=1}^{\infty} I_{k\pi/\alpha} \left( \frac{rr_n}{t'} \right) = \left( \frac{rr_n}{2t'} \right)^{\frac{\pi}{\alpha}} \sum_{k=0}^{\infty} \sum_{m=0}^{\infty} \frac{1}{m! \Gamma \left( m + \frac{(k+1)\pi}{\alpha} + 1 \right)} \left( \frac{rr_n}{2t'} \right)^{2m + \frac{k\pi}{\alpha}} \\ \le \left( \frac{rr_n}{2t'} \right)^{\frac{\pi}{\alpha}} \sum_{m=0}^{\infty} \frac{1}{m!} \left( \frac{rr_n}{2t'} \right)^{2m} \sum_{k=0}^{\infty} \frac{1}{\Gamma \left( \frac{(k+1)\pi}{\alpha} + 1 \right)} \left( \frac{rr_n}{2t'} \right)^{\frac{k\pi}{\alpha}} \le \left( \frac{rr_n}{2t'} \right)^{\frac{\pi}{\alpha}} e^{\left( \frac{rr_n}{2t'} \right)^2} \frac{\alpha}{\pi} e^{\frac{rr_n}{2t'}}.$$

With this inequality and the change of variable $u = rr_n/(2t')$, we obtain

$$d_{\text{TV}_n}(\bar{X}_T, X_T) \le C\varepsilon + \frac{8\alpha t'}{\pi r_n^2} e^{-\frac{r^2}{2t'}} \int_0^\infty u^{1+\frac{\pi}{\alpha}} e^{u^2+u} e^{-\frac{2t'}{r_n^2}u^2} du =: C\varepsilon + I(\varepsilon).$$

The integral $I(\varepsilon)$ converges if $\gamma \equiv \gamma(\varepsilon) := 2t'/r_n^2 > 1$, which is satisfied as soon as $\varepsilon < 2$. Then, denoting $\beta := 1 + \pi/\alpha$,

$$I(\varepsilon) \le \frac{8\alpha t'}{\pi r_n^2} \int_0^\infty u^\beta e^{-(\gamma-1)u^2+u} du = \frac{8\alpha t'}{\pi r_n^2} e^{\frac{1}{4(\gamma-1)}} \int_0^\infty u^\beta e^{-(\gamma-1)(u-\frac{1}{2(\gamma-1)})^2} du \\ \le \frac{8\alpha t'}{\pi r_n^2} e^{\frac{1}{4(\gamma-1)}} \left( \int_{\frac{1}{2(\gamma-1)}}^\infty u^\beta e^{-(\gamma-1)u^2} du + \left( \frac{1}{2(\gamma-1)} \right)^\beta \int_{-\frac{1}{2(\gamma-1)}}^0 e^{-(\gamma-1)u^2} du \right) \\ \le \frac{4\alpha\gamma}{\pi} e^{\frac{1}{4(\gamma-1)}} \left( \frac{1}{\sqrt{\gamma-1}1+\beta} \int_0^\infty v^\beta e^{-v^2} dv + \frac{1}{2^\beta(\gamma-1)^{\beta+1/2}} \int_{-\frac{1}{2\sqrt{\gamma-1}}}^0 e^{-v^2} dv \right).$$

The first term converges to zero as $O(\varepsilon^{\pi/(2\alpha)})$ while the second term converges to zero as $O(\varepsilon^{1+\pi/\alpha})$ when $\varepsilon \to 0$. Taking all the elements into account we obtain that

$$d_{\text{TV}}(\bar{X}_T, X_T) \underset{\varepsilon \to 0}{=} O\left(\varepsilon^{\min(1, \frac{\pi}{2\alpha})}\right).$$

□

# 7 Adaptation of the algorithms to general processes

## 7.1 Stopped Brownian motion with constant drift

Let $b \in \mathbb{R}^2$ and consider the Brownian motion with drift: $\widetilde{W}_t := W_t + bt$. Then by Girsanov's theorem, $\widetilde{W}$ is a Brownian motion under the probability

$$Q_{|\mathcal{F}_t}^{x_0} := e^{-b \cdot W_t - b^2 t/2} P_{|\mathcal{F}_t}^{x_0}, \qquad (7.1)$$

so that for any bounded test function $f: \mathcal{D} \times [0, T] \to \mathbb{R}$:

$$\mathbb{E}^{x_0} \left[ e^{-b \cdot \widetilde{W}_{T \wedge \tilde{\tau}}} + \frac{b^2 T \wedge \tilde{\tau}}{2} f(\widetilde{W}_{T \wedge \tilde{\tau}}, T \wedge \tilde{\tau}) \right] = \mathbb{E}^{x_0} [f(W_{T \wedge \tau}, T \wedge \tau)].$$

So we have:

$$\mathbb{E}^{x_0} \left[ f(\widetilde{W}_{T \wedge \tilde{\tau}}, T \wedge \tilde{\tau}) \right] = \mathbb{E}^{x_0} \left[ e^{b \cdot W_{T \wedge \tau} - \frac{b^2 T \wedge \tau}{2}} f(W_{T \wedge \tau}, T \wedge \tau) \right].$$

To simulate the Brownian motion with drift stopped in the wedge $\mathcal{D}$, we proceed as follows:
1. Simulate $W_{T \wedge \tau}$ and $T \wedge \tau$ for the stopped process without drift according to Algorithm I in Section 5.2.

2. The result of the simulation becomes

$$e^{-\frac{b^2 T \wedge \tau}{2}} e^{b \cdot W_{T \wedge \tau}} f(W_{T \wedge \tau}) \qquad (7.2)$$

We can also deduce an explicit formula for the density of the stopped Brownian motion with drift using (3.4). In fact, for $y \in \mathcal{D}$, we have

$$\mathbb{P}^{x_0}(\tilde{W}_t \in dy, \tilde{\tau} > t) = \frac{2r}{t\alpha} e^{-\frac{b^2 t}{2}} e^{b \cdot y} e^{-(r^2+r_0^2)/2t} \sum_{n=1}^{\infty} I_{n\pi/\alpha} \left(\frac{rr_0}{t}\right) \sin\left(\frac{n\pi\theta}{\alpha}\right) \sin\left(\frac{n\pi\theta_0}{\alpha}\right) drd\theta, \quad (7.3)$$

with the same notations as before: $\tilde{\tau} := \inf\{s : W_s + bs \in \partial\mathcal{D}\}$, $x = (r_0 \cos(\theta_0), r_0 \sin(\theta_0))$ and $y = (r \cos(\theta), r \sin(\theta))$.

## 7.2 Adaptation of the simulation algorithms for Itō processes

In this section, we present some extensions of the above simulation methods applied to the approximation of stopped and reflected diffusions at some fixed time $T > 0$. Let us start, first with the case of killed diffusion.

Consider the following Itō process in $\mathbb{R}^2$:

$$\begin{cases} dY_t = b(Y_t)dt + \sigma(Y_t)dW_t \\ Y_0 = x_0, \end{cases} \qquad (7.4)$$

where $(W_t)_{t \ge 0}$ is a standard two-dimensional Brownian motion and $x_0 \in \mathcal{D}$. We define $\tau$ to be the exit time of $(Y_t)_{t \ge 0}$ from $\mathcal{D}$.

**Simulation of $Y_{T \wedge \tau}$:** Choose $t_i := Ti/n$, $i = 0, ..., n$, a uniform partition of $[0, T]$ and consider the following Euler-Maruyama scheme. Given $\bar{Y}_{t_k} \in \text{int}(\mathcal{D})$ (with $\bar{Y}_{t_0} = x_0$) and supposing that $\bar{\tau} := \inf\{s; \bar{Y}_s \in \partial\mathcal{D}\} > t_k$, simulate $\bar{Y}_{t_{k+1}}$ as the final value of the following Brownian motion with drift, stopped on the wedge $\mathcal{D}$, $t \in [t_k, t_{k+1}]$:

$$\bar{Y}_{t \wedge \bar{\tau}} := \bar{Y}_{t_k} + b(\bar{Y}_{t_k})(t - t_k) + \sigma(\bar{Y}_{t_k}) \cdot (W_{t \wedge \bar{\tau}} - W_{t_k}), \quad (7.5)$$

This is done using the algorithm in Section 7.1. More precisely, we simulate $\bar{Y}_{T \wedge \bar{\tau}}$ taking $b \equiv 0$ while keeping track of the Brownian increments; then instead of $f(\bar{Y}_{T \wedge \bar{\tau}})$, the result of the simulation becomes

$$\exp\left(-\frac{1}{2} \int_{0}^{T \wedge \bar{\tau}} |\sigma^{-1}(\bar{Y}_s) \cdot b(\bar{Y}_s)|^2 ds + \int_{0}^{T \wedge \bar{\tau}} (\sigma^{-1}(\bar{Y}_s) \cdot b(\bar{Y}_s)) \cdot dW_s\right) f(\bar{Y}_{T \wedge \bar{\tau}}). \quad (7.6)$$

Note that we have to take into account the decorrelation step and the change of angle as described in Section 3 at every step. If the process (7.5) exits the wedge $\mathcal{D}$ in the interval $[t_k, t_{k+1}]$ then the algorithm directly stops. This scheme has weak order one. That is, for any $f, b, \sigma \in C_b^5(\bar{\mathcal{D}})$ one has that there exists a constant $C_f > 0$

$$|\mathbb{E}^{x_0}[f(Y_{T \wedge \tau})] - \mathbb{E}^{x_0}[f(\bar{Y}_{T \wedge \bar{\tau}})]| \le C_f n^{-1}.$$

The proof is straightforward if one follows the same line of proof as in [Gob01]. In particular, see Section 2.2.1 and note that our case is simpler as the proposed scheme does not have the possibility of touching the boundary before it is stopped or reaches $T$. For a discussion about the associated partial differential equation with Dirichlet conditions, see the discussion right after (B.1).

Since the Girsanov change of measure uses an exponential function, it may not be suitable to processes with large drifts as it may increase the variance. For this reason, we also propose a two-step Euler-Maruyama scheme for the stopped process.

**Simulation of $Y_{T \wedge \tau}$:** Choose $t_i := Ti/n$, $i = 0, ..., n$, a uniform partition of $[0, T]$ and consider the following two-step Euler-Maruyama scheme. Given $\bar{Y}_{t_k} \in \text{int}(\mathcal{D})$ (with $\bar{Y}_{t_0} = x_0$) and supposing that $\bar{\tau} := \inf\{s; \bar{Y}_s \in \partial\mathcal{D}\} > t_k$ and that $\tilde{\tau} := \inf\{s; \tilde{Y}_s \in \partial\mathcal{D}\} > t_k$, simulate $\bar{Y}_{t_{k+1} \wedge \bar{\tau}}$ in two steps as follows for $t \in [t_k, t_{k+1}]$:

$$\tilde{Y}_{t \wedge \bar{\tau}} := \bar{Y}_{t_k} + b(\bar{Y}_{t_k})(t \wedge \bar{\tau} - t_k) \qquad (7.7)$$
$$ \bar{Y}_{t \wedge \bar{\tau}} := \bar{Y}_{t_{k+1} \wedge \bar{\tau}} + (W_{t \wedge \bar{\tau}} - W_{t_k}). \qquad (7.8) $$

The second step is performed using the algorithm in Section 5.2. Note that we have to take into account the decorrelation step and the change of angle as described in Section 3. If the process exits the wedge $\mathcal{D}$ in the interval $[t_k, t_{k+1}]$ then the algorithm directly stops. This scheme has weak order one. That is, for any $f, b, \sigma \in C_b^5(\bar{D})$, there exists a constant $C_f > 0$ such that

$$ |\mathbb{E}^{x_0}[f(Y_{T \wedge \bar{\tau}})] - \mathbb{E}^{x_0}[f(\bar{Y}_{T \wedge \bar{\tau}})]| \leq C_f n^{-1}. \qquad (7.9) $$

The proof is similar to the proof for the convergence of the two-step Euler scheme for the reflected processes in Proposition 10.

Now let $(Z_t)_{t \ge 0}$ to be the solution of a stochastic differential equation with drift coefficient $b$ which is normally reflected on the boundary of the wedge $\mathcal{D}$. That is,

$$ dZ_t = b(Z_t)dt + dW_t + dL_t. $$

For details on existence and uniqueness for this equation we refer to [Sai87].

In this case, we do not know how to extend the simulation method using the Girsanov change of measure; instead we propose a two-step scheme.

**Simulation of $Z_T$:** Choose $t_i := Ti/n$, $i = 0, ..., n$ a uniform partition of $[0, T]$ and consider the following two-step scheme. Given $\bar{Z}_{t_k}$ (with $\bar{Z}_{t_0} = x_0$), simulate $\bar{Z}_{t_{k+1}}$ as the final value of the reflection on the wedge $\mathcal{D}$ of the Brownian motion and drift in two steps as follows for $t \in [0, t_{k+1} - t_k]$

$$ \bar{Z}_{t+t_k} := \bar{Z}_{t_k} + b(\bar{Z}_{t_k})(t-t_k) + L_t^1 - L_{t_k}^1 \qquad (7.10) $$

$$ \bar{Z}_{t+t_k} := \tilde{Z}_{t_{k+1}} + (W_t - W_{t_k}) + L_t^2 - L_{t_k}^2. \qquad (7.11) $$

The terms $L^i$, $i = 1, 2$, denote the respective local time terms for each of the two steps. The simulation of the second step of this algorithm uses the argument described in either Sections 5.3 or 6.3. The weak rate of convergence for these methods is as follows:

**Proposition 10.** Assume that $f, b \in C_b^5(\bar{D})$ then

$$ |\mathbb{E}^{x_0}[f(Z_T)] - \mathbb{E}^{x_0}[f(\bar{Z}_T)]| \leq C_f n^{-1}. \qquad (7.12) $$

Denoting by $\hat{Z}_T$, the result of the algorithm using the approximation described in Section 6.3, we have

$$ |\mathbb{E}^{x_0}[f(Z_T)] - \mathbb{E}^{x_0}[f(\hat{Z}_T)]| \leq C_f \left(n^{-1} + \varepsilon^{\min(1, \frac{\pi}{2\sigma})}\right). \qquad (7.13) $$

*Proof.* The proof of (7.13) is a small modification of (7.12). Therefore, we prove the rate of convergence for the algorithm provided in Section 5.3.

The argument follows as in the classical diffusion case which can be found in the proof of the main result in [TT90]. We use this argument and refer the reader to [TT90] for more details. First, consider $u$ to be the solution of the following backward partial differential equation with Neumann conditions:

$$ \begin{aligned} \partial_t u(t, x) + \mathcal{L}u(t, x) &= 0, && (t, x) \in [0, T] \times \mathcal{D} \\ u(T, x) &= f(x), && x \in \mathcal{D} \\ \nabla u(t, x) \cdot n(x) &= 0, && t > 0, x \in \partial \mathcal{D}. \end{aligned} \qquad (7.14) $$

Here, $\mathcal{L}u(t, x) := b(x)\nabla u(t, x) + \frac{1}{2}\Delta u(t, x)$. Following the same discussion as in (B.1), one obtains that $u \in C^{7/2,7}([0, T] \times \bar{D})$. This property will be used when doing Taylor's expansions for $u$. Next, note that

$$ \mathbb{E}^{x_0}[f(Z_T)] - \mathbb{E}^{x_0}[f(\bar{Z}_T)] = \sum_{i=0}^{n-1} \left( \mathbb{E}^{x_0}[u(t_i, \bar{Z}_{t_i}) - u(t_{i+1}, \bar{Z}_{t_{i+1}})] \right) =: - \sum_{i=0}^{n-1} A_i. \qquad (7.15) $$
We will now consider each term within the above sum. Using Itō-Tanaka formula (here one uses (7.14) to cancel all local time terms), we have

$$
\begin{aligned}
A_i &= \mathbb{E}^{x_0} \left[ u(t_i, \tilde{Z}_{t_{i+1}}) - u(t_i, \bar{Z}_{t_i}) + \int_{t_i}^{t_{i+1}} \left( \partial_t + \frac{1}{2}\Delta \right) u(s, \bar{Z}_s) ds \right] \\
&= \int_{t_i}^{t_{i+1}} \mathbb{E}^{x_0} \left[ b(\bar{Z}_{t_i}) \nabla u(t_i, \tilde{Z}_s) + \left( \partial_t + \frac{1}{2}\Delta \right) u(s, \bar{Z}_s) \right] ds.
\end{aligned}
$$

The following step relies on using the Taylor expansion for each of the terms involving $\nabla u(t_i, \tilde{Z}_s)$ and $(\partial_t + (1/2)\Delta) u(s, \bar{Z}_s)$ at $(t_i, \bar{Z}_{t_i})$ together with the fact that $u$ solves (7.14). Without going into much detail, let us consider the expansion of the term $\partial_t u(s, \bar{Z}_s)$:

$$
\begin{aligned}
\partial_t u(s, \bar{Z}_s) ={}& \partial_t u(t_i, \bar{Z}_{t_i}) + \int_0^1 \nabla \partial_t u(t_i, \alpha \bar{Z}_s + (1-\alpha) \bar{Z}_{t_i}) \cdot (\bar{Z}_s - \bar{Z}_{t_i}) d\alpha \\
& + \int_0^1 \partial_t \partial_t u(\alpha s + (1-\alpha)t_i, \bar{Z}_s) \cdot (s - t_i) d\alpha
\end{aligned}
$$

Each of the above derivatives of $u$ is bounded. In the case of the increments of $\bar{Z}_s - \bar{Z}_{t_i} = (W_s - W_{t_i}) + L_s^2 - L_{t_i}^2 + \tilde{Z}_{t_{i+1}} - \bar{Z}_{t_i}$ one has to take one further step in the Taylor expansion to be able to use the fact that the expectation of the Brownian increment is zero. Namely,

$$
\begin{aligned}
& \int_0^1 \nabla \partial_t u(t_i, \alpha \bar{Z}_s + (1-\alpha) \bar{Z}_{t_i}) \cdot (\bar{Z}_s - \bar{Z}_{t_i}) d\alpha \\
& = \nabla \partial_t u(t_i, \bar{Z}_{t_i}) \cdot (\bar{Z}_s - \bar{Z}_{t_i}) + \int_0^1 (1-\alpha)(\bar{Z}_s - \bar{Z}_{t_i})^\top \cdot \nabla^2 \partial_t u(t_i, \alpha \bar{Z}_s + (1-\alpha) \bar{Z}_{t_i}) \cdot (\bar{Z}_s - \bar{Z}_{t_i}) d\alpha.
\end{aligned}
$$

After this and canceling the first terms of the expansion using (7.14) and taking expectations we see that the sum of most terms is of order $O(n^{-2})$. The remaining terms are bounded by $C \int_{t_i}^{t_{i+1}} (L_s^j - L_{t_i}^j) ds$ for $j = 1, 2$. Considering these terms within (7.15), one obtains for $L_T^j = \sum_{i=0}^{n-1} (L_{t_{i+1}}^j - L_{t_i}^j)$ that these sums are bounded by

$$
\sum_{i=0}^{n-1} \int_{t_i}^{t_{i+1}} \mathbb{E}^{x_0} [L_s^j - L_{t_i}^j] ds \le n^{-1} \mathbb{E}^{x_0} [L_T^j] \le C n^{-1},
$$

where we used that $t \mapsto L_t^j$ is increasing. The last inequality follows from an estimate for the Skorohod problem in convex domains (see [Sai87, Theorem 4.2]). In the case of the algorithm described in Section 6.3 one may modify the above arguments adding at each time $t_i$, the possibility of ending the simulation in the partition interval using the criteria provided in (6.8). Using Proposition 9 one obtains the result. $\square$

# 8 Simulations

We implement the algorithms³ described in Sections 5.2 and 5.3 as well as the approximated version in Section 6.3. We also give simulation results obtained with Metzler's algorithm [Met09, Section 2.1]. This last algorithm simulates $W_\tau$ (instead of $W_{\tau \wedge T}$) directly for a wedge of general angle $\alpha$, by:

1. Simulate the radius $|W_\tau|$ using (5.6).

2. Conditionally simulate $\tau$ by acceptance-rejection, approximating the infinite sum by a partial sum, and using a Cauchy distribution as reference density.

However, Metzler's algorithm cannot be adapted to simulate $W_{\tau \wedge T}$ and is biased for the simulation of $\tau$ (see [Met09, Proposition 2.1.8]).

In Tables 1 and 2, we consider the Monte Carlo estimation of various expectations. We choose positive test functions in order to avoid cancellations and $r_0$ so that $\mathbb{E}^{x_0}[\tau]$ is of the same order as $T$, so that the stopped and reflected processes are significantly different than the standard Brownian motion.
<table><caption>Table 1: Simulations with $\alpha = 0.9$, $r_0 = 1.5$, $\theta_0 = 0.3$, $f(x, y) = x^2 + y^2$, $T = 1$.</caption><thead><tr><th></th><th>E<sup>x0</sup></th><th>95 % interval</th><th>Time (s)</th><th>MC iterations</th><th>E<sup>x0</sup>[N]</th></tr></thead><tbody><tr><td>E<sup>x0</sup>[f(W<sub>τ</sub>)] with [Met09]</td><td>3.473</td><td>± 0.082</td><td>0.67</td><td>50000</td><td>1</td></tr><tr><td>E<sup>x0</sup>[f(W<sub>τ</sub>)]</td><td>3.489</td><td>± 0.085</td><td>1.04</td><td>50000</td><td>1.45</td></tr><tr><td>E<sup>x0</sup>[τ] with [Met09]</td><td>0.742</td><td>± 0.030</td><td>7.06</td><td>20000</td><td>-</td></tr><tr><td>E<sup>x0</sup>[τ]</td><td>0.590</td><td>± 0.024</td><td>4.92</td><td>20000</td><td>-</td></tr><tr><td>E<sup>x0</sup>[f(W<sub>τ∧T</sub>)]</td><td>2.980</td><td>± 0.049</td><td>5.79</td><td>10000</td><td>1.37</td></tr><tr><td>E<sup>x0</sup>[(W<sub>τ∧T</sub>) · e<sub>1</sub>]</td><td>1.441</td><td>± 0.012</td><td>6.26</td><td>10000</td><td>1.37</td></tr><tr><td>E<sup>x0</sup>[f(X<sub>T</sub>)], ε = 0</td><td>4.138</td><td>± 0.224</td><td>13.2</td><td>1000</td><td>75.6</td></tr><tr><td>E<sup>x0</sup>[f(X<sub>T</sub>)], ε = 0.03</td><td>4.313</td><td>± 0.072</td><td>28.2</td><td>10000</td><td>5.11</td></tr></tbody></table>

<table><caption>Table 2: Simulations with $\alpha = 0.58$, $r_0 = 3$, $\theta_0 = 0.4$, $f(r \cos(\theta), r \sin(\theta)) = \sin^2(\theta)$, $T = 1$.</caption><thead><tr><th></th><th>E<sup>x0</sup></th><th>95 % interval</th><th>Time (s)</th><th>MC iterations</th><th>E<sup>x0</sup>[N]</th></tr></thead><tbody><tr><td>E<sup>x0</sup>[f(W<sub>τ∧T</sub>)]</td><td>0.195</td><td>± 0.003</td><td>8.85</td><td>10000</td><td>1.28</td></tr><tr><td>E<sup>x0</sup>[f(X<sub>T</sub>)], ε = 0.03</td><td>0.117</td><td>± 0.003</td><td>9.26</td><td>5000</td><td>2.73</td></tr></tbody></table>

We note that for the estimation of $E[f(W_τ)]$, our simulation method finds a similar value as the algorithm proposed by Metzler, as both are exact simulation methods. However, for the estimation of $E[τ]$, the bias in the method proposed by Metzler seems to be significant.

As a second way to check our algorithm, we also simulated the projection on the first axis of $W_{τ∧T}$. Note that the simulation gives a result close to the initial point $1.5 \cdot \cos(0.3) \approx 1.4330$ which is reported on the sixth line of Table 1. This is because $t \mapsto W_{τ∧t}$ is a martingale.

For the estimation of $E[f(X_T)]$ (reflected Brownian motion), the exact algorithm takes too much time, as hinted in Proposition 7, hence the need to use the approximation version from Section 6.3. We only give a simulation example with 1000 samples, as we could not get a more proper estimation. Indeed, if we increase the number of samples, there appear trajectories where $N$ becomes large and where the algorithm does not stop in reasonable time.

<figure><img src="image_6.png" /><figcaption>Figure 8.1: (left) Histogram of $N$ for the exact reflected algorithm with 50000 Monte Carlo iterations. The horizontal axis represents the values of $N$ while the vertical axis gives the number of trajectories such that $N$ is in the interval noted in the horizontal axis. The last bar counts all the values greater than 150. On the right, average number of iterations in function of $\varepsilon$. The parameters are the same as in Table 1 and for each value of $\varepsilon$, 50000 Monte Carlo iterations are used.</figcaption></figure>

To have a better idea of the behavior of the algorithm, we provide on the left side of Figure 8.1, the histogram of $N$ for the exact Algorithm II from Section 5.3; if, for a trajectory, $N$ exceeds 150, then we stop the algorithm. We observe that for most of the trajectories, $N$ is not too large (does not exceed 5). However, for a few iterations (around 0.2%) which corresponds to the case where the radius $r$ becomes small, $N$ is large, which slows down the

<footnote n="3">³The programs with a demonstration notebook are available at https://github.com/Bras-P/simulation_reflected_brownian_motion_wedge</footnote>
<table><caption>Table 3: Simulation of the process $dY_t = -\mu(Y_t - \kappa) + \sigma dW_t$, with $\alpha = 0.9$, $r_0 = 1.5$, $\theta_0 = 0.3$, $\mu = (0.1, 0.2)$, $\kappa = (0.7, 0.5)$, $\sigma = I_2$, $f(x, y) = x^2 + y^2$, $T = 1$. Z corresponds to the reflected process in the same domain. The first stopped process is simulated using the Girsanov method while the second one using the two-step Euler-Maruyama scheme.</caption><thead><tr><th></th><th>E<sup>x0</sup></th><th>95 % interval</th><th>Time (s)</th><th>MC iterations</th><th>Discretization step</th></tr></thead><tbody><tr><td>E<sup>x0</sup>[f(Y<sub>τ∧T</sub>)] [Girsanov]</td><td>2.83</td><td>± 0.06</td><td>456</td><td>5000</td><td>T/500</td></tr><tr><td>E<sup>x0</sup>[f(Y<sub>τ∧T</sub>)] [Two-step]</td><td>2.84</td><td>± 0.06</td><td>457</td><td>5000</td><td>T/500</td></tr><tr><td>E<sup>x0</sup>[f(Z<sub>T</sub>)], ε = 0.01</td><td>3.77</td><td>± 0.09</td><td>1140</td><td>5000</td><td>T/500</td></tr></tbody></table>

algorithm. On the right side of Figure 8.1 we trace $E^{x_0}[N]$ for the approximated reflected algorithm in function of $\varepsilon$. However, we could not trace the bias in function of $\varepsilon$, as it seems negligible in comparison with the size of the confidence interval. This is why we do not need to lengthen the simulation time of each trajectory by choosing a very small $\varepsilon$; $\varepsilon = 0.02$ or $\varepsilon = 0.05$ with $E^{x_0}[N] \approx 5$ is sufficient.

In Table 3, we implement the algorithms adapted to general Itō processes from Section 7.2.

# A Appendix: Auxiliary Lemmas

We start with a simple lemma about the simulation of a r.v. $X$ using only partial information of a previously simulated r.v. $Y$.

**Proposition 11.** Let $X$ and $Y$ be two random variables. We want to simulate $X$. Let $A$ be a deterministic set. We consider the following algorithm :

1. Simulate $Y$.

2. If $Y \in A$, simulate $X$ conditionally to the value of $Y$ in the previous step, and return $X$.

3. If $Y \notin A$, simulate $X$ conditionally to the event $Y \notin A$, and return $X$.

Let us denote $\tilde{X}$ the value returned by the algorithm. Then:

$$ \tilde{X} \stackrel{\mathcal{L}}{=} X. $$

*Proof.* The random variable $\tilde{X}$ is defined conditionally to $Y$ as:

$$ \mathbb{P}(\tilde{X} \in dx | Y = y) = \begin{cases} \mathbb{P}(X \in dx | Y = y) & \text{if } y \in A \\ \mathbb{P}(X \in dx | Y \notin A) & \text{if } y \notin A \end{cases} $$

Then:

$$ \begin{aligned} \mathbb{P}(\tilde{X} \in dx) &= \int \mathbb{P}(\tilde{X} \in dx | Y = y) \mathbb{P}(Y \in dy) \\ &= \int_{y \in A} \mathbb{P}(X \in dx | Y = y) \mathbb{P}(Y \in dy) + \int_{y \notin A} \mathbb{P}(X \in dx | Y \notin A) \mathbb{P}(Y \in dy) \\ &= \mathbb{P}(X \in dx, Y \in A) + \mathbb{P}(X \in dx | Y \notin A) \mathbb{P}(Y \notin A) = \mathbb{P}(X \in dx). \end{aligned} $$

□

## A.1 Estimates on Bessel processes

**Proposition 12.** For all $x \ge 0$, we have:

$$ 1 \le I_0(x) \le e^x. \qquad (A.1) $$

Furthermore, for $x \ge 0$ and $\nu \ge 0$, we have:

$$ I_{\nu}(x) \le e^{x} + \frac{1}{\pi(\nu + x)}. \qquad (A.2) $$
Proof. The inequality follows from the integral representation [Wat44, p.181]:

$$I_{\nu}(x) = \frac{1}{\pi} \int_{0}^{\pi} e^{x \cos \theta} \cos(\nu \theta) d\theta - \frac{\sin(\nu \pi)}{\pi} \int_{0}^{\infty} e^{-x \cosh(t) - \nu t} dt.$$

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

Now we provide two results on estimates related to the law of $F(T) = \int_{0}^{T} \frac{ds}{R_{s}^{2}}$.

**Proposition 13.** Let $(R_t)_{t \ge 0}$ be a Bessel process of dimension 2, or equivalently, of index 0. Then

$$\mathbb{P}^{r_0} [F(T) \ge x] \underset{x \to \infty}{\sim} \frac{1}{\sqrt{2\Gamma(1/2)}} \left( \int_{\frac{r_0^2}{2T}}^{\infty} \frac{e^{-u}}{u} du \right) x^{-1/2}. \quad (A.3)$$

Proof. Step 1: We use [JYC09], Proposition 6.2.5.1, with $a = 0$ and $\nu = 0$:

$$\mathbb{E}^{r_0} \left[ e^{-\frac{\lambda^2}{2} F(T)} \right] = \frac{r_0^{\lambda}}{\Gamma(\lambda/2)} \int_{0}^{\infty} v^{\frac{\lambda}{2}-1} (1+2vT)^{-(1+\lambda)} e^{-\frac{r_0^2 v}{1+2vT}} dv.$$

Performing the change of variable $w := \frac{Tv}{1+2vT}$, so that $v = \frac{w}{T(1-2w)}$, we have for all $\lambda > 0$:

$$\mathbb{E}^{r_0} \left[ e^{-\frac{\lambda^2}{2} F(T)} \right] = \frac{r_0^{\lambda}}{\Gamma(\lambda/2) T^{\lambda/2}} \int_{0}^{1/2} w^{\frac{\lambda}{2}-1} (1-2w)^{\lambda/2} e^{-r_0^2 w/T} dw. \quad (A.4)$$

Step 2: In fact, we prove the asymptotic expansion:

$$\mathbb{E}^{r_0} \left[ e^{-\frac{\lambda^2}{2} F(T)} \right] \underset{\lambda \to 0^+}{=} 1 - \frac{\lambda}{2} \int_{\frac{r_0^2}{2T}}^{\infty} \frac{e^{-u}}{u} du + O(\lambda^2). \quad (A.5)$$

We have for all $x > 0$, $\Gamma(x+1) = x\Gamma(x)$ and $\Gamma(1) = 1$ so that

$$\Gamma(x) \underset{x \to 0^+}{\sim} \frac{1}{x}.$$

Then, using (A.4) and the change of variable $w = Tu/r_0^2$:

$$\mathbb{E}^{r_0} \left[ e^{-\frac{\lambda^2}{2} F(T)} \right] = \frac{1}{\Gamma\left(\frac{\lambda}{2}\right)} \int_{0}^{\frac{r_0^2}{2T}} u^{\frac{\lambda}{2}-1} \left(1 - \frac{2Tu}{r_0^2}\right)^{\frac{\lambda}{2}} e^{-u} du.$$

Using the definition of the Gamma function and decomposing the integral in two parts, we have

$$1 - \mathbb{E}^{r_0} \left[ e^{-\frac{\lambda^2}{2} F(T)} \right] = \frac{1}{\Gamma\left(\frac{\lambda}{2}\right)} \left[ \int_{0}^{\frac{r_0^2}{2T}} u^{\frac{\lambda}{2}-1} e^{-u} \left( 1 - \left( 1 - \frac{2Tu}{r_0^2} \right)^{\frac{\lambda}{2}} \right) du + \int_{\frac{r_0^2}{2T}}^{\infty} u^{\frac{\lambda}{2}-1} e^{-u} du \right] =: I_1 + I_2.$$

We use the expansion $(1-x)^\alpha = \sum_{k=0}^\infty (\alpha \log(1-x))^k/k!$ for $|x| < 1$ and $\alpha \to 0$ in order to analyze the first integral. This gives the following asymptotic equivalence

$$I_1 = \underset{\lambda \to 0^+}{-\frac{\lambda^2}{4}} \int_0^{\frac{r_0^2}{2T}} u^{-1} e^{-u} \log \left( 1 - \frac{2Tu}{r_0^2} \right) du + o(\lambda^2) = O(\lambda^2).$$

And on the other side:

$$I_2 = \frac{1}{\Gamma\left(\frac{\lambda}{2}\right)} \int_{\frac{r_0^2}{2T}}^{\infty} u^{\frac{\lambda}{2}-1} e^{-u} du = \underset{\lambda \to 0^+}{\frac{\lambda}{2}} \int_{\frac{r_0^2}{2T}}^{\infty} \frac{e^{-u}}{u} du + O(\lambda^2).$$

From here we obtain the asymptotic expansion (A.5).

Step 3: We get (A.3) from the expansion (A.5) and using Karamata Tauberian theorems (see [EKM97], Corollary A3.10):
**Proposition 14.** Let df be a measure on $\mathbb{R}^+$, $F(x) := \int_0^x df$, $\hat{f}(\lambda) := \int_0^\infty e^{-\lambda x} F(x)dx$, $0 \le \gamma < 1$ and let $L > 0$. Then the following are equivalent:

$$1. \ 1 - \hat{f}(\lambda) \sim L\lambda^{\gamma}, \quad \lambda \to 0^{+}$$

$$2. \ 1 - F(x) \sim \frac{L}{\Gamma(1 - \gamma)} x^{-\gamma}, \quad x \to \infty.$$

□

Next we will give an estimate for $\mathbb{E}[F(T \wedge \zeta_{\varepsilon})^p]$ in the case of the modified algorithm for reflected Brownian motion.

**Lemma 15.** We define the stopping time $\zeta_{\varepsilon} := \inf\{t \in [0, T] : R_t^2/(T-t) \le \varepsilon\}$. Then we have the following bound for $p \in (1, 2)$ and $\varepsilon > 0$ small enough:

$$\mathbb{E}^{r_0} [F(T \wedge \zeta_{\varepsilon})^p] \le C \left(1 + \log\left(\frac{2T}{r_0^2}\right)\right) (\varepsilon T)^{1-p}. \quad (A.6)$$

*Proof.* By [JYC09, (6.2.3)], the density of $R_t$ in $r$ is

$$\frac{r}{t} \left(\frac{r}{r_0}\right)^{\nu} e^{-\frac{r^2+r_0^2}{2t}} I_{\nu}\left(\frac{rr_0}{t}\right), \quad (A.7)$$

where $\nu = \frac{d}{2} - 1$ and $d$ is the dimension of the Bessel process. Taking $d = 2$, we have:

$$\begin{aligned} \mathbb{E}^{r_0} [F(T \wedge \zeta_{\varepsilon})^p] &\le \int_0^T \mathbb{E}^{r_0} \left[ \frac{1}{R_t^{2p}} \mathbf{1}_{R_t \ge \sqrt{\varepsilon(T-t)}} \right] dt = \int_0^T \int_{\sqrt{\varepsilon(T-t)}}^{\infty} \frac{1}{r^{2p-1}t} e^{-\frac{r^2+r_0^2}{2t}} I_0\left(\frac{rr_0}{t}\right) dr dt \\ &= \int_{\sqrt{\varepsilon T}}^{\infty} \int_0^T \frac{1}{r^{2p-1}t} e^{-\frac{r^2+r_0^2}{2t}} I_0\left(\frac{rr_0}{t}\right) dt dr + \int_0^{\sqrt{\varepsilon T}} \int_{T-r^2/\varepsilon}^T \frac{1}{r^{2p-1}t} e^{-\frac{r^2+r_0^2}{2t}} I_0\left(\frac{rr_0}{t}\right) dt dr =: I_1 + I_2. \end{aligned}$$

Then, using (A.1) for both $I_1$ and $I_2$:

$$\begin{aligned} I_1 &\le \int_{\sqrt{\varepsilon T}}^{\infty} \frac{1}{r^{2p-1}} \int_0^T \frac{1}{t} e^{-\frac{(r-r_0)^2}{2t}} dt dr = \int_{\sqrt{\varepsilon T}}^{\infty} \frac{1}{r^{2p-1}} \int_{\frac{(r-r_0)^2}{2T}}^{\infty} \frac{e^{-v}}{v} dv dr \\ &= \int_{\sqrt{\varepsilon T}}^{\sqrt{2T}+r_0} \frac{1}{r^{2p-1}} \int_{\frac{(r-r_0)^2}{2T}}^{\infty} \frac{e^{-v}}{v} dv dr + \int_{\sqrt{2T}+r_0}^{\infty} \frac{1}{r^{2p-1}} \int_{\frac{(r-r_0)^2}{2T}}^{\infty} \frac{e^{-v}}{v} dv dr. \end{aligned}$$

But we have, for all $a \ge 0$:

$$\int_a^\infty \frac{e^{-x}}{x} dx \le \begin{cases} e^{-a} & \text{if } a \ge 1, \\ e^{-1} + \log(1/a) & \text{if } a < 1. \end{cases}$$

So that

$$I_1 \le \int_{\sqrt{\varepsilon T}}^{\sqrt{2T}+r_0} \frac{dr}{r^{2p-1}} \left( e^{-1} + \log \left( \frac{2T}{(r-r_0)^2} \right) \right) + \int_{\sqrt{2T}+r_0}^{\infty} \frac{dr}{r^{2p-1}} e^{-\frac{(r-r_0)^2}{2T}} \underset{\varepsilon \to 0}{\sim} \frac{(\varepsilon T)^{1-p}}{2(p-1)} \left( e^{-1} + \log \left( \frac{2T}{r_0^2} \right) \right).$$

For the second integral, we have:

$$\begin{aligned} I_2 &\le \int_0^{\sqrt{\varepsilon T}} \frac{1}{r^{2p-1}} \int_{T-r^2/\varepsilon}^T e^{-\frac{(r-r_0)^2}{2t}} \frac{dt}{t} dr = \int_0^{\sqrt{\varepsilon T}} \frac{1}{r^{2p-1}} \int_{\frac{(r-r_0)^2}{2T}}^{\frac{(r-r_0)^2}{2(T-r^2/\varepsilon)}} \frac{e^{-v}}{v} dv dr \\ &\le \int_0^{\sqrt{\varepsilon T}} \frac{1}{r^{2p-1}} \int_{\frac{(r-r_0)^2}{2T}}^{\frac{(r-r_0)^2}{2(T-r^2/\varepsilon)}} \frac{dv}{v} dr \\ &= - \int_0^{\sqrt{\varepsilon T}} \frac{1}{r^{2p-1}} \log \left(1 - \frac{r^2}{\varepsilon T}\right) dr = - \frac{1}{(\varepsilon T)^{p-1}} \int_0^1 \frac{1}{u^{2p-1}} \log(1-u^2) du, \end{aligned}$$

where the last integral converges as soon as $p < 2$. □
B Proofs in Section 4

Proof of Theorem 3 in the case that $\mathcal{D} = \langle \pi/m \rangle$. We apply the method of images in $\mathbb{R}^2$. Let $f: \mathcal{D} \to [0, +\infty)$ be a non-negative continuous function with compact support. The heat equation with boundary conditions

$$ 
\begin{array}{ll}
\partial_t u(t, x) = \frac{1}{2} \Delta u(t, x), & (t, x) \in [0, +\infty) \times \mathcal{D} \\
u(0, x) = f(x), & x \in \mathcal{D} \\
\nabla u(t, x) \cdot n(x) = 0, & t > 0, x \in \partial \mathcal{D},
\end{array}
\qquad (B.1) $$

where $n(x)$ denotes the inward unitary orthogonal vector on the boundary, can be rewritten as a partial differential equation in polar coordinates $(r, \theta)$ as described in [MNP00, Chapter 1]. Furthermore if we perform the change of variable $z = \log(r)$, we obtain a parabolic problem with non-constant coefficients in a strip with mixed type boundary condition. The existence and uniqueness for this PDE is treated in [LSU68, chapter IV, Theorem 5.3]. From this reference and the Feyman-Kac representation theorem we obtain that

$$ u(t, x) := \mathbb{E}^x [f(X_t)] $$

satisfies the heat equation on $\mathcal{D}$ described above. Now we define the function

$$ 
\begin{array}{ll}
\tilde{f} : \mathbb{R}^2 \to \mathbb{R} \\
y \mapsto f(T_k^{-1} y) \quad \text{for } y \in \mathcal{D}_k,
\end{array}
$$

and for $t \ge 0$ and $x \in \mathbb{R}^2$,

$$ \tilde{u}(t, x) := \mathbb{E}^x [\tilde{f}(W_t)]. $$

Then $\tilde{u}$ satisfies the heat equation on $\mathbb{R}^2$ with $\tilde{f}$ as initial condition, so that

$$ \tilde{u}(t, x) = \frac{1}{2\pi t} \int_{\mathbb{R}^2} \tilde{f}(y) e^{-\frac{|x-y|^2}{2t}} dy = \frac{1}{2\pi t} \sum_{k=0}^{2m-1} \int_{\mathcal{D}_k} \tilde{f}(y) e^{-\frac{|x-y|^2}{2t}} dy = \frac{1}{2\pi t} \int_{\mathcal{D}} f(y) \sum_{k=0}^{2m-1} e^{-\frac{|x-T_k y|^2}{2t}} dy. $$

On the other hand, $u$ and $\tilde{u}$ satisfy the same boundary equation on $\mathcal{D}$. Indeed, for all $x \in \mathcal{D}$, $\tilde{f}(x) = f(x)$ so $u$ and $\tilde{u}$ satisfy the same initial conditions. For $x \in \partial \mathcal{D}^{-}$, $k = 0, \dots, 2m-1$ and $y \in \mathcal{D}$, we have

$$ 
\begin{array}{l}
(x - T_k y) \cdot n(x) = -(x - T_{2m-1-k} y) \cdot n(x), \\
|x - T_k y|^2 = |x - T_{2m-1-k} y|^2,
\end{array}
$$

so that

$$ 
\begin{aligned}
\nabla \tilde{u}(t, x) \cdot n(x) &= - \frac{1}{2\pi t^2} \int_{\mathcal{D}} f(y) \sum_{k=0}^{2m-1} (x - T_k y) \cdot n(x) e^{-\frac{|x-T_k y|^2}{2t}} dy \\
&= - \frac{1}{2\pi t^2} \int_{\mathcal{D}} f(y) \sum_{k=0}^{m-1} (x - T_k y) \cdot n(x) e^{-\frac{|x-T_k y|^2}{2t}} dy \\
&\quad + \frac{1}{2\pi t^2} \int_{\mathcal{D}} f(y) \sum_{k=0}^{m-1} (x - T_{2m-1-k} y) \cdot n(x) e^{-\frac{|x-T_{2m-1-k} y|^2}{2t}} dy = 0.
\end{aligned}
$$

If $x \in \partial \mathcal{D}^{+}$, we get the same result noting that for all $y \in \mathcal{D}$ and $k = 0, \dots, 2m-1$,

$$ 
\begin{array}{l}
(x - T_{(k+1) \bmod 2m} y) \cdot n(x) = -(x - T_{(2m-k) \bmod 2m} y) \cdot n(x), \\
|x - T_{(k+1) \bmod 2m} y|^2 = |x - T_{(2m-k) \bmod 2m} y|^2,
\end{array}
$$

So $u(t, x) = \tilde{u}(t, x)$ for all $t \ge 0$ and $x \in \mathcal{D}$, and

$$ \mathbb{E}^x [f(X_t)] = \frac{1}{2\pi t} \int_{\mathcal{D}} f(y) \sum_{k=0}^{2m-1} e^{-\frac{|x-T_k y|^2}{2t}} dy. $$

□
Proof of Theorem 3 in the general case. As in the above proof we consider

$$u(t, x) = \mathbb{E}^x [f(X_t)], \qquad (B.2)$$

where $t \ge 0$, $x \in \mathcal{D}$, and $f: \mathcal{D} \to [0, +\infty)$ is a non-negative continuous function with compact support. Then $u$ is the solution of the partial differential equation (B.1). Considering the formula obtained for the case $\alpha = \pi/m$ above, we would like to express it so that it does not depend explicitly on $m$. We first assume that $\alpha = \pi/m$ for some $m \in \mathbb{N}$. In order to switch to polar coordinates, let $x = (r_0 \cos(\theta_0), r_0 \sin(\theta_0))$, $y = (r \cos(\theta), r \sin(\theta))$ and recall that $\vartheta_k$ for $0 \le k \le 2m-1$ is the angle of $T_k(y)$ in $[0, 2\pi)$, i.e. $T_k(y) = (r \cos(\vartheta_k), r \sin(\vartheta_k))$ with $\vartheta_{2k} = 2k\alpha + \theta$ and $\vartheta_{2k+1} = 2(k+1)\alpha - \theta$. We then rewrite (4.3):

$$\mathbb{P}^x(X_t \in dy) = \frac{r}{2\pi t} e^{-\frac{r^2+r_0^2}{2t}} \sum_{k=0}^{2m-1} e^{\frac{rr_0}{t}} \cos(\theta_0 - \vartheta_k) drd\theta. \qquad (B.3)$$

We use the following identity (see [GR07, page 933, (8.511.4)]), valid for $\gamma, z \ge 0$:

$$e^{\gamma z} = I_0(z) + 2 \sum_{n=1}^{\infty} T_n(\gamma) I_n(z),$$

where $T_n$ is the $n^{th}$ Tchebychev's polynomial of the first kind and $I_n$ the modified Bessel function of the first kind with order $n$. Then the result (3.4) follows because $T_n(\cos(\theta)) = \cos(n\theta)$ and

$$\sum_{k=0}^{2m-1} \cos(n(\theta_0 - \vartheta_k)) = \begin{cases} 2m \cos(n\theta) \cos(n\theta_0) & \text{if } n \text{ is a multiple of } m, \\ 0 & \text{otherwise.} \end{cases}$$

Using the above properties and $m = \pi/\alpha$, we obtain

$$\mathbb{P}^x(X_t \in dy) = \frac{2r}{t\alpha} e^{-(r^2+r_0^2)/(2t)} \left( \frac{1}{2} I_0\left(\frac{rr_0}{t}\right) + \sum_{n=1}^{\infty} I_{n\pi/\alpha}\left(\frac{rr_0}{t}\right) \cos\left(\frac{n\pi\theta}{\alpha}\right) \cos\left(\frac{n\pi\theta_0}{\alpha}\right) \right) drd\theta, \qquad (B.4)$$

which gives a formula which does not depend on $m$. Now, we have to check that this formula is well defined for any $\alpha \in (0, \pi)$ and that it is the unique solution the partial differential equation (B.1) with $\mathcal{D} = \langle \alpha \rangle$.

This is explained in Corollary 16. In our current situation, $d=2$, the eigenvalues are $\lambda_j := (j\pi/\alpha)^2$ for $j \ge 0$ and the eigenfunctions are $m_j(\theta) = \sqrt{2/\alpha} \cos(j\pi\theta/\alpha)$, $j \ge 1$ and $m_0 = 1/\sqrt{\alpha}$. $\square$

**Corollary 16.** Consider a general $d$-dimensional cone generated by all rays emanating from the origin and passing through a compact subset $D \subset \mathbb{S}^{d-1}$ which has smooth boundary. Consider $X$ to be the normally reflected Brownian motion at the boundary of the cone. Then $X_t$ has a density given by

$$\mathbb{P}^x(X_t \in dy) = \frac{r}{t rr_0^{d/2-1}} e^{-\frac{r^2+r_0^2}{2t}} \left( I_{\alpha_0}\left(\frac{rr_0}{t}\right) m_0(\theta) m_0(\theta_0) + \sum_{n=1}^{\infty} I_{\alpha_n}\left(\frac{rr_0}{t}\right) m_n(\theta) m_n(\theta_0) \right) drd\theta. \qquad (B.5)$$

To explain the elements in the above formula, denote by $L_{\mathbb{S}^{d-1}}$, the Laplace-Beltrami operator on $\mathbb{S}^{d-1}$. With the above assumptions, there exists a complete set of orthonormal eigenfunctions $m_j$ with corresponding eigenvalues $0 \le \lambda_0 < \lambda_1 \le \lambda_2 < \dots$ satisfying

$$\left\{ \begin{array}{ll} L_{\mathbb{S}^{d-1}} m_j(x) = -\lambda_j m_j(x) & \text{for } x \in D \\ \nabla m_j(x) \cdot n(x) = 0 & \text{for } x \in \partial D, \end{array} \right.$$

$$\alpha_j = \sqrt{\lambda_j + \left(\frac{d}{2} - 1\right)^2}.$$
*Proof.* The beginning of the proof is the same as in the statements following (B.1) in what refers to the existence and uniqueness of the associated PDE. Existence and uniqueness for the reflected process in the generalized cone can be deduced from [Bas96] (see also [BBC05, Remark 4.1] and the references therein).

In order to prove that (B.5) satisfies the associated PDE, one follows a similar proof in [BnS97] for the killed case. The proof in the reflected case follows line by line, the proof in [BnS97], except that the required estimates and properties of the eigenvalues and eigenfunctions of the Laplace-Beltrami operator in the case of Neumann boundary conditions have to be referred to the proper literature (for this see, [GN13], [Gri02] and [Kro92]) although the estimates do not change as far as it relates to the proof of [BnS97, Lemma 1] with the corresponding corrections of typos. In order for the required estimates to be satisfied one needs that the generalized $d$-dimensional cone is generated by all rays emanating from the origin and passing through a compact subset $D \subset \mathbb{S}^{d-1}$ which has smooth boundary.

We also remark that the density expressions in [BnS97] are written under polar measures which explains why our expressions have an extra $r$ which appears due to the Jacobian of the change of coordinates. $\square$

We also provide a full elementary proof that (B.4) satisfies the initial conditions in the Supplementary Material.

## C A hint for higher dimensions

The following proposition is provided so as to hint at the possibilities in dimension higher than two. It shows that dimension two is the case which is mathematically difficult to treat.

**Proposition 17.** Let $(R_t)_{t \ge 0}$ be a Bessel process in dimension $d \ge 3$ and let $F(T) = \int_0^t ds/R_s^2$. Then:
$$ \mathbb{E}^{r_0}[F(T)] < +\infty. $$

*Proof.* The density of $R_t$ in $r$ is given in (A.7), so by performing the change of variable $r = ut$, we have

$$ \mathbb{E}^{r_0}[F(T)] = \int_0^T \int_0^\infty \frac{1}{rt} \left(\frac{r}{r_0}\right)^\nu e^{-\frac{r^2+r_0^2}{2t}} I_\nu\left(\frac{rr_0}{t}\right) drdt = \frac{1}{r_0^\nu} \int_0^T t^{\nu-1} \int_0^\infty I_\nu(\nu r_0) e^{-\frac{tu^2}{2}} e^{-\frac{r_0^2}{2t}} u^{\nu-1} du dt. $$

Note that for all $\varepsilon > 0$,

$$ \int_{\varepsilon}^{T} t^{\nu-1} \int_{0}^{\infty} I_{\nu}(ur_{0})e^{-\frac{tu^{2}}{2}-\frac{r_{0}^{2}}{2t}}u^{\nu-1}du dt \leq \left(\int_{\varepsilon}^{T} e^{-\frac{r_{0}^{2}}{2t}}t^{\nu-1}dt\right) \left(\int_{0}^{\infty} I_{\nu}(ur_{0})e^{-\varepsilon u^{2}/2}u^{\nu-1}du\right) < \infty, $$

where for we used the Proposition A.2 for the convergence of the second integral. So to prove the convergence of the integral, we only need to prove the convergence for the integral in $t$ around zero. To do so we use Proposition A.2 again so that for $\nu$, $r_0 > 0$:

$$ \mathbb{E}^{r_0}[F(\varepsilon)] \le \frac{1}{r_0^\nu} \int_0^\varepsilon t^{\nu-1} e^{-\frac{r_0^2}{2t}} \int_0^\infty \left( e^{ur_0} + \frac{1}{\pi(\nu + ur_0)} \right) e^{-\frac{tu^2}{2}} u^{\nu-1} du dt < \infty. $$

In fact, using that $\nu > 0$, we have

$$ \begin{aligned} I_1 &:= \frac{1}{r_0^\nu} \int_0^\varepsilon t^{\nu-1} \int_0^\infty e^{-\frac{t}{2}\left(u-\frac{r_0}{t}\right)^2} u^{\nu-1} du dt \le \frac{1}{r_0^\nu} \int_0^\varepsilon t^{\nu-1} \int_{-\infty}^\infty e^{-\frac{t}{2}u^2} \left(|u|+\frac{r_0}{t}\right)^{\nu-1} du dt \\ &= \frac{2}{r_0^\nu} \int_0^\varepsilon t^{\nu-3/2} \int_0^\infty e^{-u^2/2} \left(\frac{u}{\sqrt{t}}+\frac{r_0}{t}\right)^{\nu-1} du dt \le \frac{2}{r_0^\nu} \int_0^\varepsilon t^{\nu-3/2} \int_0^\infty e^{-u^2/2} \left(\frac{u+r_0}{t}\right)^{\nu-1} du dt \\ &= \frac{2}{r_0^\nu} \int_0^\varepsilon t^{\nu-1/2} \int_0^\infty e^{-u^2/2} (u+r_0)^{\nu-1} du dt < \infty. \end{aligned} $$

As for the second integral, we have

$$ I_2 := \frac{1}{r_0^\nu} \int_0^\varepsilon t^{\nu-1} \int_0^\infty \frac{1}{\pi(\nu + ur_0)} e^{-\frac{tu^2}{2} - \frac{r_0^2}{2t}} u^{\nu-1} du dt \le \frac{1}{r_0^\nu} \int_0^\varepsilon t^{\nu-1} e^{-\frac{r_0^2}{2t}} dt \int_0^\infty \frac{u^{\nu-1}}{\pi(\nu + ur_0)} du < \infty. \square $$
# Acknowledgements

The first author thanks École Normale Supérieure, Département de Mathématiques et Applications for providing
a financial support for a visit to Ritsumeikan University.

The second author is supported in part by KAKENHI 20K03666.

Both authors thank the anonymous referee for a careful review and for correcting mistakes and providing
references.

# References

[Bas96] Richard F. Bass. Uniqueness for the Skorokhod equation with normal reflection in Lipschitz domains. *Electron. J. Probab.*, 1:no. 11, approx. 29 pp. 1996.

[BBC05] Richard F. Bass, Krzysztof Burdzy, and Zhen-Qing Chen. Uniqueness for reflecting Brownian motion in lip domains. *Ann. Inst. H. Poincaré Probab. Statist.*, 41(2):197–235, 2005.

[BC08] Krzysztof Burdzy and Zhen-Qing Chen. Discrete approximations to reflected Brownian motion. *Ann. Probab.*, 36(2):698–727, 2008.

[BC15] Jose Blanchet and Xinyun Chen. Steady-state simulation of reflected Brownian motion and related stochastic networks. *Ann. Appl. Probab.*, 25(6):3209–3250, 2015.

[BGT04] Mireille Bossy, Emmanuel Gobet, and Denis Talay. A symmetrized Euler scheme for an efficient approximation of reflected diffusions. *J. Appl. Probab.*, 41(3):877–889, 2004.

[BM18] Jose Blanchet and Karthyek Murthy. Exact simulation of multidimensional reflected Brownian motion. *J. Appl. Probab.*, 55(1):137–156, 2018.

[BnS97] Rodrigo Bañuelos and Robert G. Smits. Brownian motion in cones. *Probab. Theory Related Fields*, 108(3):299–319, 1997.

[BSCD13] Christophette Blanchet-Scalliet, Areski Cousin, and Diana Dorobantu. Hitting time for correlated three-dimensional Brownian motion. *HAL*, 2013. hal-00846450v2.

[BST10] Christian Bayer, Anders Szepessy, and Raúl Tempone. Adaptive weak approximation of reflected and stopped diffusions. *Monte Carlo Methods Appl.*, 16(1):1–67, 2010.

[CBM15] Marie Chupeau, Olivier Bénichou, and Satya N. Majumdar. Survival probability of a Brownian motion in a planar wedge of arbitrary angle. *Phys. Rev. E (3)*, 91(3):032106, 8, 2015.

[CJ59] H. S. Carslaw and J. C. Jaeger. *Conduction of Heat in Solids*. Oxford University Press, second edition, 1959.

[CPS98] C. Costantini, B. Pacchiarotti, and F. Sartoretto. Numerical approximation for functionals of reflecting diffusion processes. *SIAM J. Appl. Math.*, 58(1):73–102, 1998.

[Die10] A. B. Dieker. Reflected brownian motion. *Wiley Encyclopedia of Operations Research and Management Science*, 2010.

[DL06] Madalina Deaconu and Antoine Lejay. A random walk on rectangles algorithm. *Methodol. Comput. Appl. Probab.*, 8(1):135–151, 2006.

[DM09] A. B. Dieker and J. Moriarty. Reflected Brownian motion in a wedge: sum-of-exponential stationary densities. *Electron. Commun. Probab.*, 14:1–16, 2009.

[Dub04] Julien Dubédat. Reflected planar Brownian motions, intertwining relations and crossing probabilities. *Ann. Inst. H. Poincaré Probab. Statist.*, 40(5):539–552, 2004.
[EFW13] Marcos Escobar, Sebastian Ferrando, and Xianzhang Wen. Three dimensional distribution of Brownian motion extrema. *Stochastics*, 85(5):807–832, 2013.

[EKM97] Paul Embrechts, Claudia Klüppelberg, and Thomas Mikosch. *Modelling extremal events*, volume 33 of *Applications of Mathematics (New York)*. Springer-Verlag, Berlin, 1997. For insurance and finance.

[GN13] D. S. Grebenkov and B.-T. Nguyen. Geometrical structure of Laplacian eigenfunctions. *SIAM Rev.*, 55(4):601–667, 2013.

[GNR86] V. Giorno, A. G. Nobile, and L. M. Ricciardi. On some diffusion approximations to queueing systems. *Adv. in Appl. Probab.*, 18(4):991–1014, 1986.

[Gob01] Emmanuel Gobet. Euler schemes and half-space approximation for the simulation of diffusion in a domain. *ESAIM Probab. Statist.*, 5:261–297, 2001.

[GR07] I. S. Gradshteyn and I. M. Ryzhik. *Table of integrals, series, and products*. Elsevier/Academic Press, Amsterdam, seventh edition, 2007. Translated from the Russian, Translation edited and with a preface by Alan Jeffrey and Daniel Zwillinger, With one CD-ROM (Windows, Macintosh and UNIX).

[Gri02] D. Grieser. Uniform bounds for eigenfunctions of the Laplacian on manifolds with boundary. *Comm. Partial Differential Equations*, 27(7-8):1283–1299, 2002.

[Ha09] Wonho Ha. *Applications of the reflected Ornstein-Uhlenbeck process*. ProQuest LLC, Ann Arbor, MI, 2009. Thesis (Ph.D.)–University of Pittsburgh.

[IKP13] Tomoyuki Ichiba, Ioannis Karatzas, and Vilmos Prokaj. Diffusions with rank-based characteristics and values in the nonnegative quadrant. *Bernoulli*, 19(5B):2455–2493, 2013.

[Iye85] Satish Iyengar. Hitting lines with two-dimensional Brownian motion. *SIAM J. Appl. Math.*, 45(6):983–989, 1985.

[JYC09] Monique Jeanblanc, Marc Yor, and Marc Chesney. *Mathematical methods for financial markets*. Springer Finance. Springer-Verlag London, Ltd., London, 2009.

[Kag07] Wouter Kager. Reflected Brownian motion in generic triangles and wedges. *Stochastic Process. Appl.*, 117(5):539–549, 2007.

[KLR18] Vadim Kaushansky, Alexander Lipton, and Christoph Reisinger. Transition probability of Brownian motion in the octant and its application to default modelling. *Appl. Math. Finance*, 25(5-6):434–465, 2018.

[Kro92] P. Kroger. Upper bounds for the Neumann eigenvalues on a bounded domain in Euclidean space. *J. Funct. Anal.*, 106(2):353–357, 1992.

[KZ16] Steven Kou and Haowen Zhong. First-passage times of two-dimensional Brownian motion. *Adv. in Appl. Probab.*, 48(4):1045–1060, 2016.

[LSU68] O. A. Ladyženskaja, V. A. Solonnikov, and N. N. Ural'ceva. *Linear and quasilinear equations of parabolic type*. Translations of Mathematical Monographs, Vol. 23. American Mathematical Society, Providence, R.I., 1968. Translated from the Russian by S. Smith.

[Met09] Adam Metzler. *Multivariate first-passage models in credit risk*. ProQuest LLC, Ann Arbor, MI, 2009. Thesis (Ph.D.)–University of Waterloo (Canada).

[Met10] Adam Metzler. On the first passage problem for correlated Brownian motion. *Statist. Probab. Lett.*, 80(5-6):277–284, 2010.

[MNP00] Vladimir Maz'ya, Serguei Nazarov, and Boris Plamenevskij. *Asymptotic theory of elliptic boundary value problems in singularly perturbed domains*. Vol. I, volume 111 of *Operator Theory: Advances and Applications*. Birkhäuser Verlag, Basel, 2000. Translated from the German by Georg Heinig and Christian Posthoff.
[Pil14] Andrey Pilipenko. *An introduction to stochastic differential equations with reflection*. Number 1 in Lectures in pure and applied mathematics (1). Universitätsverlag Potsdam, 2014.

[RY99] Daniel Revuz and Marc Yor. *Continuous martingales and Brownian motion*, volume 293 of *Grundlehren der mathematischen Wissenschaften* [Fundamental Principles of Mathematical Sciences]. Springer-Verlag, Berlin, third edition, 1999.

[Sai87] Yasumasa Saisho. Stochastic differential equations for multidimensional domain with reflecting boundary. *Probab. Theory Related Fields*, 74(3):455–477, 1987.

[TT90] Denis Talay and Luciano Tubaro. Expansion of the global error for numerical schemes solving stochastic differential equations. *Stochastic Anal. Appl.*, 8(4):483–509 (1991), 1990.

[VW85] S. R. S. Varadhan and R. J. Williams. Brownian motion in a wedge with oblique reflection. *Comm. Pure Appl. Math.*, 38(4):405–443, 1985.

[Wat44] G. N. Watson. *A Treatise on the Theory of Bessel Functions*. Cambridge University Press, Cambridge, England; Macmillan Company, New York, 1944.
# D Supplementary material

## D.1 Proof of the convergence to the initial condition when $t \to 0$

In this section, we prove the following result:

**Proposition 18.** The formula (4.4) satisfies the initial conditions of the heat equation, i.e. for all functions $f: \mathcal{D} \to \mathbb{R}$ continuous with compact support, $r_0 > 0$ and $\theta_0 \in (0, \alpha)$:

$$ \int_{\mathcal{D}} \hat{f}(r, \theta) \frac{2r}{t\alpha} e^{-(r^2+r_0^2)/2t} \left( \frac{1}{2} I_0 \left( \frac{rr_0}{t} \right) + \sum_{n=1}^{\infty} I_{n\pi/\alpha} \left( \frac{rr_0}{t} \right) \cos \left( \frac{n\pi\theta}{\alpha} \right) \cos \left( \frac{n\pi\theta_0}{\alpha} \right) \right) dr d\theta \underset{t \to 0}{\longrightarrow} \hat{f}(r_0, \theta_0), \quad (\text{D.1}) $$

where $\hat{f}: [0, \infty) \times [0, \alpha] \to \mathbb{R}$ denotes the function $f$ expressed in polar coordinates.

For the proof, we use the following representation:

$$ \forall x > 0, \forall \alpha \ge 0, I_{\alpha}(x) = \frac{1}{\pi} \int_{0}^{\pi} e^{x \cos(u)} \cos(\alpha u) du - \frac{\sin(\alpha \pi)}{\pi} \int_{0}^{\infty} e^{-x \cosh(u) - \alpha u} du. $$

Let us denote by $A_1$ and $A_2$ the two terms appearing in the integral representation in (D.1) after replacing with the above formula, with $d_n = \frac{1}{2}$ if $n = 0$ and $d_n = 1$ for $n \ge 1$:

$$ A_1 := \frac{2}{\alpha\pi} \int_{\mathcal{D}} \frac{r}{t} \hat{f}(r, \theta) e^{-\frac{r^2+r_0^2}{2t}} \sum_{n=0}^{\infty} d_n \cos\left(\frac{n\pi\theta}{\alpha}\right) \cos\left(\frac{n\pi\theta_0}{\alpha}\right) \int_0^{\pi} e^{\frac{rr_0}{t}\cos(u)} \cos\left(\frac{n\pi}{\alpha}u\right) du dr d\theta, $$

$$ A_2 := -\frac{2}{\alpha\pi} \int_{\mathcal{D}} \frac{r}{t} \hat{f}(r, \theta) e^{-\frac{r^2+r_0^2}{2t}} \sum_{n=1}^{\infty} \sin\left(\frac{n\pi^2}{\alpha}\right) \cos\left(\frac{n\pi\theta}{\alpha}\right) \cos\left(\frac{n\pi\theta_0}{\alpha}\right) \int_0^{\infty} e^{-\frac{rr_0}{t}\cosh(u) - \frac{n\pi}{\alpha}u} du dr d\theta. $$

**Proposition 19.** We have $A_2 \underset{t \to 0}{\longrightarrow} 0$.

*Proof. Step 1:* For the moment let us assume that the integration order can be switched. For all $a \in \mathbb{R}$ and $b > 0$, we have:

$$ g(a, b) := \sum_{n=0}^{\infty} \sin(na)e^{-nb} = \Im \left( \sum_{n=0}^{\infty} e^{ina-nb} \right) = \frac{e^{-b} \sin(a)}{(1 - e^{-b} \cos(a))^2 + e^{-2b} \sin^2(a)}. \quad (\text{D.2}) $$

We use the following trigonometric identity written in compact form:

$$ 4 \cos(a) \cos(b) \sin(c) = \sin(a+b+c) + \sin(c-a-b) + \sin(a-b+c) + \sin(c-a+b) \quad (\text{D.3}) $$

$$ A_{jk}(\theta) := \frac{\pi}{\alpha} ((-1)^j \theta + (-1)^k \theta_0 + \pi), $$

$$ \sum_{n=1}^{\infty} \cos\left(\frac{n\pi\theta}{\alpha}\right) \cos\left(\frac{n\pi\theta_0}{\alpha}\right) \sin\left(\frac{n\pi^2}{\alpha}\right) e^{-\frac{n\pi}{\alpha}u} = \frac{1}{4} \sum_{j,k=0}^{1} g\left(A_{jk}(\theta), \frac{\pi}{\alpha}u\right). $$

Now we prove that for all $j, k \in \{0, 1\}$, we have:

$$ J_{j,k} := \frac{1}{t} \int_{\mathcal{D}} r \hat{f}(r, \theta) e^{-\frac{r^2+r_0^2}{2t}} \int_{0}^{\infty} e^{-\frac{rr_0}{t} \cosh(u)} g\left(A_{jk}(\theta), \frac{\pi}{\alpha}u\right) du d\theta dr \underset{t \to 0}{\longrightarrow} 0. $$

We will do the analysis of $J_{j,k}$ dividing the integration region in two: $J_{j,k}^1$ which comprises the integral on the region $\mathcal{D} \times (1, \infty)$ and the remaining which is denoted by $J_{j,k}^2$. With this in mind, note that $g$ is continuous on $[0, \infty) \times (0, \alpha)$ and is locally integrable in $(0, 0)$. In fact, for all $a \in \mathbb{R}$ and $b > 0$:

$$ g(a, b) = \frac{e^{-b} \sin(a)}{(1 - e^{-b} \cos(a))^2 + e^{-2b} \sin^2(a)} \underset{a, b \to 0}{\sim} \frac{a}{\left(\left(1 - e^{-b}\left(1 - \frac{a^2}{2}\right)\right)^2 + a^2\right)} \sim \frac{a}{a^2 + b^2}. \quad (\text{D.4}) $$
Moreover,

$$ \int_{[0,1]^2} \frac{a}{a^2 + b^2} dadb = \int_0^1 \left[ \frac{1}{2} \log(a^2 + b^2) \right]_{a=0}^{a=1} db = \frac{1}{2} \int_0^1 (\log(1+b^2) - \log(b^2)) db < \infty. $$

Next, note that $\partial_b g(a, b)$ is negative, so

$$ \forall u \ge 1, \forall \theta \in \mathbb{R}, \quad \left|g\left(A_{jk}(\theta), \frac{\pi}{\alpha}u\right)\right| \le \frac{e^{-\frac{\pi}{\alpha}}}{(1 - e^{-\frac{\pi}{\alpha}})^2}, $$

and since $u \mapsto u^2/\cosh(u)$ is non-negative and bounded above, there exists $B > 0$ such that

$$ |J_{jk}^1| \le \frac{C_1}{t} \int_0^\infty re^{-\frac{r^2+r_0^2}{2t}} \int_1^\infty e^{-B\frac{rr_0}{t}} u^2 dudr \le \frac{C_2}{\sqrt{tr_0}} e^{-\frac{r_0^2}{2t}} (2t)^{3/4} \int_0^\infty \sqrt{r'} e^{-r'2} dr' \underset{t \to 0}{\longrightarrow} 0. $$

On the other hand, using (D.4) and the continuity of $g$ we obtain

$$ |J_{jk}^2| \le ||f||_\infty \frac{1}{t} \int_0^\infty re^{-\frac{r^2+r_0^2}{2t}} dr \cdot \int_0^\alpha |g\left(A_{jk}(\theta), \frac{\pi}{\alpha}u\right)| d\theta du \underset{t \to 0}{\longrightarrow} 0. $$

So that

$$ |A_2| \le \frac{1}{2\alpha\pi} \sum_{j,k=0}^{1} |J_{jk}| \underset{t \to 0}{\longrightarrow} 0. $$

Step 2: Now, we prove that one can interchange the order of the integrals and sum in $A_2$. Note that Fubini's theorem does not apply here because of the factors $\cos(n\pi\theta/\alpha)$ and $\cos(n\pi\theta_0/\alpha)$. A computation similar to (D.2) leads to, for all $N \in \mathbb{N}$:

$$ g_N(a, b) := \sum_{n=N}^{\infty} \sin(na)e^{-nb} = \Im \left( \sum_{n=N}^{\infty} e^{-nb+ina} \right) = e^{-Nb} \frac{\sin(Na) - e^{-b}\sin((N-1)a)}{(1 - e^{-b}\cos(a))^2 + e^{-2b}\sin^2(a)}. $$

Next, we prove that, for all $j, k \in \{0, 1\}$ and for all $t > 0$, denoting $R_{jk}^N$ the difference between the infinite sum and the partial sum up in $A_2$ to $N$,

$$ R_{jk}^N := \frac{1}{t} \int_D re^{-\frac{r^2+r_0^2}{2t}} \hat{f}(r, \theta) \int_0^\infty e^{-\frac{rr_0}{t}} \cosh(u) g_N \left( A_{jk}(\theta), \frac{\pi}{\alpha}u \right) du d\theta dr \underset{N \to \infty}{\longrightarrow} 0. $$

Let us remark that

$$ e^{Nb}g_N(a, b) = \frac{\sin(Na) - e^{-b}\sin((N-1)a)}{(1 - e^{-b}\cos(a))^2 + e^{-2b}\sin^2(a)} \underset{a,b \to 0}{\sim} \frac{a}{a^2 + b^2}, $$

so that $(a, b) \mapsto e^{Nb}g_N(a, b)$ is integrable in $(0, 0)$. We denote $R_{jk}^{N,1}$ and $R_{jk}^{N,2}$ the two terms obtained after splitting the integral with respect to $u$ on $(0, 1)$ and $(1, \infty)$ respectively. We have then

$$ |R_{jk}^{N,2}| \le \frac{1}{t} ||f||_\infty \int_0^\infty re^{-\frac{r^2+r_0^2}{2t}} dr \int_0^1 e^{-\frac{N\pi u}{\alpha}} \int_0^\alpha e^{-\frac{N\pi u}{\alpha}} |g_N(A_{jk}(\theta), \frac{\pi}{\alpha}u)| d\theta du \underset{N \to \infty}{\longrightarrow} 0, $$

where we use Lemma 20 for the above convergence. Moreover for all $u \ge 1$, $\theta \in \mathbb{R}$,

$$ |g_N(A_{jk}(\theta), \frac{\pi}{\alpha}u)| \le e^{-N\frac{\pi}{\alpha}u} \frac{2}{(1 - e^{-\frac{\pi}{\alpha}})^2}, $$

so we have

$$ |R_{jk}^{N,1}| \le \frac{1}{t} ||f||_\infty \int_0^\infty re^{-\frac{r^2+r_0^2}{2t}} dr \int_1^\infty e^{-N\frac{\pi}{\alpha}u} \frac{2}{(1 - e^{-\frac{\pi}{\alpha}})^2} du \underset{N \to \infty}{\longrightarrow} 0. $$

From the above arguments we obtain the conclusion: $R_{jk}^N \to 0$ for all $j, k \in \{0, 1\}$. $\square$

<figure><img src="image_8.png" /></figure>
**Lemma 20.** Let $b > 0$ and let $f \in L^1((0, b))$ be non-negative and continuous. Then:

$$ \int_0^b e^{-Nx} f(x)dx \underset{N \to \infty}{\longrightarrow} 0. $$

*Proof.* Let $\varepsilon > 0$ and choose $\delta > 0$ such that $\int_0^\delta f(x)dx \le \varepsilon$. Then for $N$ big enough:

$$ \int_0^b e^{-Nx} f(x)dx \le \int_0^\delta f(x)dx + e^{-N\delta}||f||_1 \le 2\varepsilon. $$

□

**Proposition 21.** We have $A_1 \underset{t \to 0}{\longrightarrow} \hat{f}(r_0, \theta_0)$.

*Proof.* First, let us assume that the integration order can be exchanged. This will be further discussed in Step 5.

**Step 1:** We exchange the order of integrals in $A_1$ as:

$$ \begin{align*} A_1 &= \int_0^\pi f_1(u,t)f_2(u,t)du, \\ f_1(u,t) &:= \frac{2}{\sqrt{t\alpha\pi}}e^{-\frac{r_0^2 \sin^2(u)}{2t}}, \\ f_2(u,t) &:= \int_0^\infty \frac{re^{-\frac{(r-r_0 \cos(u))^2}{2t}}}{\sqrt{t}} \sum_{n=0}^\infty d_n \cos\left(\frac{n\pi\theta_0}{\alpha}\right) \cos\left(\frac{n\pi}{\alpha}u\right) \int_0^\alpha \hat{f}(r,\theta) \cos\left(\frac{n\pi\theta}{\alpha}\right) d\theta dr. \end{align*} $$

We will study the limit of $A_1$ in the above integral order. First, we treat the sum inside $f_2(u,t)$ using trigonometric identities for $\cos(\frac{n\pi(\theta_0 \pm u)}{\alpha})$. We see that is enough to find the limit for:

$$ \frac{1}{2} \sum_{n=0}^{\infty} d_n \left( \cos \left( \frac{n\pi(\theta_0 + u)}{\alpha} \right) + \cos \left( \frac{n\pi(\theta_0 - u)}{\alpha} \right) \right) \int_0^\alpha \hat{f}(r, \theta) \cos \left( \frac{n\pi\theta}{\alpha} \right) d\theta. $$

In our case, we will use the following normalization of the classical Fourier inversion formula:

**Theorem 22 (Fourier Formula).** Let $g: \mathbb{R} \to \mathbb{R}$ be a periodic and continuous function of period $2L$. Then the following series converges for all $x \in \mathbb{R}$ and:

$$ \begin{gathered} g(x) = \frac{a_0}{2} + \sum_{n=1}^{\infty} a_n \cos\left(\frac{n\pi x}{L}\right) + b_n \sin\left(\frac{n\pi x}{L}\right), \\ \text{with } a_n = \frac{1}{L} \int_{-L}^{L} g(\theta) \cos\left(\frac{n\pi\theta}{L}\right) d\theta \quad \text{and} \quad b_n = \frac{1}{L} \int_{-L}^{L} g(\theta) \sin\left(\frac{n\pi\theta}{L}\right) d\theta. \end{gathered} $$

Next, we extend the definition of the function $\hat{f}$. For all $r \ge 0$ and $\theta \in [0, \alpha]$, we define $\hat{f}(r, -\theta) := \hat{f}(r, \theta)$, and we make it $2\alpha$-periodic by defining $\hat{f}(r, \theta + 2k\alpha) = \hat{f}(r, \theta)$. This way, $\hat{f}$ is an even and $2\alpha$-periodic function and is still continuous. With this definition, and using the Fourier inversion formula, we get for $\theta_0 \pm u \in (-\alpha, \alpha)$:

$$ \begin{align} & \frac{1}{2} \sum_{n=0}^{\infty} d_n \cos \left( \frac{n\pi(\theta_0 \pm u)}{\alpha} \right) \int_0^{\alpha} \hat{f}(r, \theta) \cos \left( \frac{n\pi\theta}{\alpha} \right) d\theta \nonumber \\ & = \frac{1}{4} \sum_{n=0}^{\infty} d_n \cos \left( \frac{n\pi(\theta_0 \pm u)}{\alpha} \right) \int_{-\alpha}^{\alpha} \hat{f}(r, \theta) \cos \left( \frac{n\pi\theta}{\alpha} \right) d\theta = \frac{\alpha}{4} \hat{f}(r, \theta_0 \pm u). \tag{D.5} \end{align} $$
**Step 2:** Next, let us study $f_2(u, t)$ in two separate cases the integral in $r$:

$$ \int_{0}^{\infty} \hat{f}(r, \theta_{0} \pm u) \frac{r e^{-\frac{(r-r_{0} \cos(u))^{2}}{2t}}}{\sqrt{t}} dr. $$

Case 1: If $r_0 \cos(u) > 0$:

$$ \frac{r e^{-\frac{(r-r_0 \cos(u))^{2}}{2t}}}{\sqrt{t}} = \frac{1}{\sqrt{t}}(r-r_0 \cos(u))e^{-\frac{(r-r_0 \cos(u))^{2}}{2t}} + r_0 \cos(u) \frac{e^{-\frac{(r-r_0 \cos(u))^{2}}{2t}}}{\sqrt{t}}. $$

The total mass of the first term is

$$ \int_{0}^{\infty} \frac{1}{\sqrt{t}}(r-r_0 \cos(u))e^{-\frac{(r-r_0 \cos(u))^{2}}{2t}} dr = \left[ -t \frac{e^{-\frac{(r-r_0 \cos(u))^{2}}{2t}}}{\sqrt{t}} \right]_{r=0}^{\infty} = \sqrt{t} e^{-\frac{r_0^2 \cos^2(u)}{2t}} \underset{t \to 0}{\longrightarrow} 0. $$

And the second term is, up to the multiplicative constant $\sqrt{2\pi}$, an approximation of the unity around $r = r_0 \cos(u) > 0$, so that in this case, the integral in $r$ converges to

$$ \sqrt{2\pi r_0 \cos(u)} \hat{f}(r_0 \cos(u), \theta_0 \pm u). $$

Case 2: If $r_0 \cos(u) \le 0$: The total mass of $\frac{1}{\sqrt{t}}re^{-(r-r_0 \cos(u))^2/(2t)}$ is bounded above by

$$ \frac{1}{\sqrt{t}} \int_{0}^{\infty} r e^{-\frac{r^2}{2t}} dr = \sqrt{t} \underset{t \to 0}{\longrightarrow} 0, $$

so that the integral in $r$ converges to 0. We remark here that convergences in the above two cases are uniform with respect to $u$ within their respective domains.

**Step 3:** From the previous step, we consider now the integral with respect to $u$. Taking into consideration the previous step, we can restrict to the case $\cos(u) > 0$ or equivalently $u \in [0, \pi/2]$. That is, consider

$$ I_{\pm} := \frac{1}{\sqrt{t}} \int_{0}^{\pi/2} \hat{f}(r_0 \cos(u), \theta_0 \pm u) e^{-\frac{r_0^2 \sin^2(u)}{2t}} \cos(u) du. $$

Note that for all $\varepsilon > 0$, by dominated convergence,

$$ \frac{1}{\sqrt{t}} \int_{\varepsilon}^{\pi/2} \hat{f}(r_0 \cos(u), \theta_0 \pm u) e^{-\frac{r_0^2 \sin^2(u)}{2t}} \cos(u) du \underset{t \to 0}{\longrightarrow} 0. $$

Fix $\delta > 0$ and take $\varepsilon$ small enough such that $\cos(u) \simeq 1$ and $\sin(u) \simeq u$ for all $u \in [0, \varepsilon]$. Then

$$ \frac{1}{\sqrt{t}} \int_{0}^{\varepsilon} e^{-\frac{r_{0}^{2} \sin^{2}(u)}{2t}} \cos(u) du \underset{t \to 0}{\approx} \frac{1}{\sqrt{t}} \int_{0}^{\varepsilon} e^{-\frac{r_{0}^{2} u^{2}}{2t}} du \underset{t \to 0}{\longrightarrow} \sqrt{\frac{\pi}{2r_{0}^{2}}}. $$

Thus, up to the multiplicative constant $\sqrt{\pi/(2r_0^2)}$, $u \mapsto e^{-r_0^2 \sin^2(u)/(2t)} \cos(u)$ is an approximation of the unity around $u = 0$, so that

$$ I_+ + I_- \underset{t \to 0}{\longrightarrow} \sqrt{\frac{2\pi}{r_0^2}} \hat{f}(r_0, \theta_0). $$

**Step 4:** Now, we put all previous steps together. In Step 1, we proved that for all $t > 0$:

$$ f_2(u, t) := \int_{0}^{\infty} \frac{r e^{-\frac{(r-r_0 \cos(u))^{2}}{2t}}}{\sqrt{t}} \frac{\alpha}{4} (\hat{f}(r, \theta_0 + u) + \hat{f}(r, \theta_0 - u)) dr. $$
We have proved in Step 2 that for all $u$, $f_2(u,t)$ converges when $t \to 0$ to

$$ \bar{f}_2(u) := \frac{\alpha}{4} \sqrt{2\pi r_0} \cos(u) \left( \hat{f}(r_0 \cos(u), \theta_0 + u) + \hat{f}(r_0 \cos(u), \theta_0 - u) \right), $$

and in Step 3 that $\int_0^\pi f_1(u,t)\bar{f}_2(u)du$ converges when $t \to 0$ to

$$ \begin{aligned} & \lim_{t \to 0} \frac{2}{\alpha \pi \sqrt{t}} \int_0^{\pi/2} \frac{\alpha}{4} \sqrt{2\pi r_0} (\hat{f}(r_0 \cos(u), \theta_0 + u) + \hat{f}(r_0 \cos(u), \theta_0 - u)) e^{-\frac{r_0^2 \sin^2(u)}{2t}} \cos(u) du \\ & = \frac{2}{\alpha \pi} \frac{\alpha}{4} \sqrt{2\pi r_0} \sqrt{\frac{\pi}{2r_0^2}} 2 \hat{f}(r_0, \theta_0) = \hat{f}(r_0, \theta_0). \end{aligned} $$

To end the proof of the convergence, we have to show that

$$ \lim_{t \to 0} \int_0^\pi f_1(u,t)f_2(u,t)du = \lim_{t \to 0} \int_0^\pi f_1(u,t)\bar{f}_2(u)du. $$

We have:

$$ \left| \int_0^\pi f_1(u,t)f_2(u,t)du - f_1(u,t)\bar{f}_2(u)du \right| \le \left( \sup_t \int_0^\pi f_1(u,t)du \right) \|f_2(\cdot,t) - \bar{f}_2(\cdot)\|_\infty, $$

and $\|f_2(\cdot,t) - \bar{f}_2(\cdot)\|_\infty \to 0$ since in Step 2, the convergence is uniform with respect to $u$.

**Step 5:** We now prove that for all $r_0$, $\theta_0$, $t$ fixed, the integration order can be switched. Note that for fixed $n$, by Fubini's theorem, the integration order in $r$, $\theta$ and $u$ can be switched, as

$$ \begin{aligned} & \left| \int_0^\infty \int_0^\alpha \left| \hat{f}(r, \theta) e^{-\frac{r^2+r_0^2}{2t}} \cos\left(\frac{n\pi\theta}{\alpha}\right) \cos\left(\frac{n\pi\theta_0}{\alpha}\right) \right| \int_0^\pi \left| e^{\frac{rr_0}{t}} \cos(u) \cos\left(\frac{n\pi u}{\alpha}\right) \right| du d\theta dr \right| \\ & \le \|f\|_\infty \alpha \pi \int_0^\infty e^{-\frac{r^2+r_0^2}{2t}} e^{\frac{rr_0}{t}} dr < \infty. \end{aligned} $$

So that the order of integration can be exchanged for every partial sum. Now, for $N \in \mathbb{N}$:

$$ \left| \int_0^\pi e^{-\frac{r_0^2 \sin^2(u)}{2t}} \int_0^\infty r e^{-\frac{(r-r_0 \cos(u))^2}{2t}} \sum_{n=0}^N d_n \cos\left(\frac{n\pi\theta_0}{\alpha}\right) \cos\left(\frac{n\pi}{\alpha}u\right) \cdot \int_0^\alpha \hat{f}(r, \theta) \cos\left(\frac{n\pi\theta}{\alpha}\right) d\theta dr du \right. \\ - \int_0^\pi e^{-\frac{r_0^2 \sin^2(u)}{2t}} \int_0^\infty r e^{-\frac{(r-r_0 \cos(u))^2}{2t}} \sum_{n=0}^\infty d_n \cos\left(\frac{n\pi\theta_0}{\alpha}\right) \cos\left(\frac{n\pi}{\alpha}u\right) \cdot \int_0^\alpha \hat{f}(r, \theta) \cos\left(\frac{n\pi\theta}{\alpha}\right) d\theta dr du \\ \left. \le \frac{\alpha}{4} \int_0^\infty r e^{-\frac{r^2+r_0^2}{2t}} \left( \int_{-\pi}^\pi e^{\frac{rr_0}{t}} \cos(u) |\hat{f}_N(r, \theta_0 + u) - \hat{f}(r, \theta_0 + u)| du \right) dr \right| $$

where $\hat{f}_N$ denotes the $N^{th}$ partial Fourier sum of $\hat{f}$. Then, by Cauchy-Schwarz inequality:

$$ \le \frac{\alpha}{4} \int_0^\infty r e^{-\frac{r^2+r_0^2}{2t}} \sqrt{\int_{-\pi}^\pi (\hat{f}_N(r, \theta_0 + u) - \hat{f}(r, \theta_0 + u))^2 du} \sqrt{\int_{-\pi}^\pi e^{\frac{2rr_0}{t}} \cos(u) du} dr $$

$$ \le \frac{\alpha}{4} \sqrt{\left| \frac{\pi}{\alpha} \right|} \sup_r \| \hat{f}_N(r, \cdot) - \hat{f}(r, \cdot) \|_2 \underbrace{\int_0^\infty r e^{-\frac{r^2+r_0^2}{2t}} \sqrt{2\pi e^{\frac{rr_0}{t}}} dr}_{<+\infty} $$

Using Parseval's equality, we obtain

$$ \begin{aligned} \| \hat{f}_N(r, \cdot) - \hat{f}(r, \cdot) \|_2^2 &= \sum_{n=N+1}^{\infty} \int_{-\alpha}^{\alpha} \cos\left(\frac{n\pi\xi}{\alpha}\right)^2 d\xi \frac{1}{\alpha^2} \left( \int_{-\alpha}^{\alpha} \hat{f}(r, \theta) \cos\left(\frac{n\pi\theta}{\alpha}\right) d\theta \right)^2 \\ &= \sum_{n=N+1}^{\infty} \int_{-\alpha}^{\alpha} \cos\left(\frac{n\pi\xi}{\alpha}\right)^2 d\xi \cdot \frac{1}{n^2\pi^2} \left( \int_{-\alpha}^{\alpha} \partial_{\theta} \hat{f}(r, \theta) \sin\left(\frac{n\pi\theta}{\alpha}\right) d\theta \right)^2 \end{aligned} $$
$$ \leq \sum_{n=N+1}^{\infty} \frac{4\alpha^3}{n^2\pi^2} \|\partial_{\theta}\hat{f}\|_{\infty}. $$

Although the extension of $\hat{f}$ on $[-\alpha, \alpha]$ is not an element of $C^1$, we can perform the integration by parts on each interval $[0, \alpha]$ and $[-\alpha, 0]$. Then, $\|\hat{f}_N(r, \cdot) - \hat{f}(r, \cdot)\|_2 \to 0$ as $N \to \infty$ uniformly in $r$. That way we obtain the convergence to 0 of the difference between the partial sum and the series. $\square$