2. Stochastic perpetuities

Let Y be a stable process with stability and positivity parameters $\alpha$ and $\rho$ , respectively (see Appendix A below for definition). Since $Y_0=0$ and the scaling property yield $\overline{Y}_{t}=\sup_{s\in [0,t] }Y_{s}\overset{d}{=}t^{1/\alpha}\overline{Y}_{1}$ for all $t\in[0,\infty)$ , we may restrict our attention to $\overline{Y}_1$ . Let $S (\alpha,\rho) $ and $\overline{S} (\alpha,\rho) $ denote the laws of $Y_{1}$ and $\overline{Y}_{1}$ , respectively. Since $\mathbb{P} (Y_{t}>0) =\rho$ for any $t>0$ , the extreme cases $\rho\in \{ 0,1\} $ are excluded from our analysis as they correspond to Y having monotone paths. Let U(0,1) denote the uniform law on (0,1) and define $x^+=\max\{x,0\}$ for any real number $x\in\mathbb{R}$ .

Proposition 2.1. Let $ (\overline{Y}_1,Z,U) \sim \overline{S} (\alpha,\rho) \times S (\alpha,\rho) \times U (0,1) $ . Then the law of $\overline{Y}_{1}$ is the unique solution of the following perpetuity:

(2.1) \begin{equation} \overline{Y}_{1} \overset{d}{=} U^{\frac{1}{\alpha}}\overline{Y}_{1}+ (1-U) ^{\frac{1}{\alpha}}Z^{+}. \label{eqn2} \end{equation}

To prove this result, we need the next definition. For any $a<b$ , the concave majorant of a function $f\colon [a,b]\to\mathbb{R}$ is defined as the smallest concave function $c\colon [a,b]\to\mathbb{R}$ , such that $c (t) \geq x (t) $ for every $t\in [a,b] $ . The proof of Proposition 2.1 exploits the fact that the supremum of a function lies on its concave majorant, at the end of all (if any) faces with positive slope. Following the classical result for the complete description of a concave majorant of random walks, Pitman and Uribe Bravo [Reference Pitman and Uribe Bravo27] described the continuous-time analogue of these results for Lévy processes ([Reference Pitman and Uribe Bravo27] is phrased in terms of the convex minorant, but through a change of sign their results cover the concave majorant). The idea is as follows: fix a sample path of Y and pick a random face of its concave majorant above an independent uniform point in [0,1]. The length of the chosen face is distributed as $V\sim U(0,1)$ and its height is distributed as the increment of a stable process over a time interval of duration V. Moreover, after removing this face (together with the path underneath it) the remainder of the concave majorant behaves like a concave majorant of a stable process over the time interval [0,1-V]; see [Reference Pitman and Uribe Bravo27]. This recursive relation and the scaling property of Y will yield the perpetuity in (2.1).

Proof. A stick-breaking process $ \{ \ell_{n}\} _{n\geq1}$ on [0,1] is defined recursively as follows:

\[ \ell_{n}=V_{n} (1-L_{n-1}) ,\quad n\geq1, \]

where $L_{n-1}=\ell_{1}+\cdots+\ell_{n-1}$ , $L_0=0$ and $ \{ V_{n}\} _{n\geq1}$ is a sequence of i.i.d. random variables with law U(0,1) (independent of Y). Let $C=(C_t)_{t\in[0,1]}$ be the concave majorant of the Lévy process Y. Let $ (d_{n}-g_{n},C_{d_{n}}-C_{g_{n}}) _{n\geq1}$ be the lengths and heights of the faces of C picked at random, uniformly on lengths and without replacement ( $g_n$ and $d_n$ denote the beginning and end times for the nth face). Theorem 1 of [Reference Pitman and Uribe Bravo27] asserts the following equality in law:

\[ (d_{n}-g_{n},C_{d_{n}}-C_{g_{n}}) _{n\geq1}\overset{d}{=} (\ell_{n},Y_{L_{n}}-Y_{L_{n-1}}) _{n\geq1}. \]

The concave majorant $(C_t)_{t\in[0,1]}$ is piecewise linear, with the corresponding slopes forming a non-increasing piecewise constant function in t. Hence $\overline{Y}_{1}$ is always contained in the image of the function C. Moreover, the supremum equals the sum of all the positive heights of C:

\begin{equation*} \overline{Y}_{1} = \sum_{n=1}^{\infty} (C_{d_{n}}-C_{g_{n}}) ^{+}\overset{d}{=}\sum_{n=1}^{\infty} (Y_{L_{n}}-Y_{L_{n-1}}) ^{+}. \end{equation*}

Conditional on $ \{ L_{n}\}_{n\geq1}$ , the random variables $ \{ Y_{L_{n}}-Y_{L_{n-1}}\}_{n\geq1}$ are independent and have the same distribution as the respective $Y_{\ell_{n}}$ . Hence, for an i.i.d. sequence $ \{ Z_{n}\} _{n\geq1}$ with law $S (\alpha,\rho) $ , we have

\[ (\ell_{n},Y_{Z_{n}}-Y_{Z_{n-1}}) _{n\geq1}\overset{d}{=} (\ell_{n},\ell_{n}^{\frac{1}{\alpha}}Z_{n}) _{n\geq1}, \]

implying

(2.2) \begin{equation} \overline{Y}_{1} \overset{d}{=} \sum_{n=1}^{\infty} (Y_{L_{n}}-Y_{L_{n-1}}) ^{+}\overset{d}{=}\sum_{n=1}^{\infty}\ell_{n}^{\frac{1}{\alpha}}Z_{n}^{+}. \label{eqn3} \end{equation}

It is well known that $ \{ \frac{p}{\ell_{n}}{1-\ell_{1}}\} _{n\geq2}$ is a stick-breaking process on [0,1], independent of $\ell_1\sim U(0,1)$ (and $ \{ Z_{n}\} _{n\geq1}$ ). Hence by (2.2) we find the equality in law

\[ \overline{Y}_{1}\overset{d}{=}\sum_{n=2}^{\infty}\bigg(\dfrac{\ell_{n}}{1-\ell_{1}}\bigg)^{\frac{1}{\alpha}}Z_{n}^{+}, \]

which, together with (2.2), implies the perpetuity

\begin{equation*} \overline{Y}_{1} \overset{d}{=} \ell_{1}^{\frac{1}{\alpha}}Z_{1}^{+}+ (1-\ell_{1}) ^{\frac{1}{\alpha}}\overline{Y}_{1}. \end{equation*}

Finally, the uniqueness of solution follows from [Reference Buraczewski, Damek and Mikosch5, Theorem 2.1.3].

Let $S^{+} (\alpha,\rho) $ denote the law of $Y_{1}$ conditioned on being positive. For $n,m\in\mathbb{Z}$ define the sets

(2.3) \begin{equation} \mathcal{Z}^{n} = \{k\in\mathbb{Z}: k<n\},\quad\mathcal{Z}_{m}^{n}=\mathcal{Z}^{n}\backslash \mathcal{Z}^{m}. \label{eqn4} \end{equation}

Proof of Theorem 1.1. Note that the random variable Z ⁺ in Proposition 2.1 behaves like the product of a Bernoulli random variable and a stable random variable conditioned on being positive, that is, if $B\sim \operatorname{\rm{Ber}} (\rho) $ and $S\sim S^{+} (\alpha,\rho) $ are independent, then $Z^{+}\overset{d}{=}BS$ . Since $\mathbb{P} (Z^{+}=0) =1-\rho>0$ , the idea behind the proof of Theorem 1.1 is to iterate perpetuity (2.1) backwards in time until the first time we observe Z ⁺>0.

More precisely, by Proposition 2.1 and Kolmogorov’s consistency theorem, we can construct a stationary Markov chain $ \{ (U_{n},Z_{n},\zeta_{n}) \}_{n\in\mathcal{Z}^{1}}$ with invariant law $U (0,1) \times S (\alpha,\rho) \times\overline{S} (\alpha,\rho) $ , where $ \{ (U_{n},Z_{n}) \} _{n\in\mathcal{Z}^{1}}$ is an i.i.d. sequence with law $U (0,1) \times S (\alpha,\rho) $ and

\begin{equation*} \zeta_{n+1} = U_{n}^{\frac{1}{\alpha}}Z_{n}^{+}+ (1-U_{n}) ^{\frac{1}{\alpha}}\zeta_{n},\quad n\in\mathcal{Z}^{0}. \end{equation*}

Define $V_{0}=1$ and $V_{n}=\prod_{m\in\mathcal{Z}_{n}^{0}} (1-U_{m}) $ for $n\in\mathcal{Z}^{0}$ . Then the following equality holds:

(2.4) \begin{equation} \zeta_{0}=\sum_{m\in\mathcal{Z}_{n}^{0}} (U_{m}V_{m+1}) ^{\frac{1}{\alpha}}Z_{m}^{+}+V_{n}^{\frac{1}{\alpha}}\zeta_{n}\quad \text{for\ all}\ n\in\mathcal{Z}^{0}. \label{eqn5} \end{equation}

Let $\tau=\sup \{ n\in\mathcal{Z}^{0}\colon Z_{n}>0\} $ (with convention $\sup\emptyset=-\infty$ ) be the last time we see a positive value in the sequence $ \{ Z_{n}\}_{n\in\mathcal{Z}^{0}}$ . Substituting $n=\tau$ in (2.4), we get

(2.5) \begin{equation} \zeta_{0}=V_{\tau+1}^{\frac{1}{\alpha}} ( (1-U_{\tau}) ^{\frac{1}{\alpha}}\zeta_{\tau}+U_{\tau}^{\frac{1}{\alpha}}Z_{\tau}) . \label{eqn6} \end{equation}

This equality of course yields the same equality in law. It will hence imply the perpetuity in (1.1), if we prove that the random variables involved have the desired laws and independence structure.

The events $ \{ Z_{n}>0\} $ , $n\in\mathcal{Z}^0$ , are independent with probability $\rho$ , making $\tau$ a geometric random variable on $\mathcal{Z}^0$ with parameter $\rho$ . By construction, the coordinates of the vector $(U_n,Z_n,\zeta_n)$ are independent for any $n\in\mathcal{Z}^0$ . Hence we have $ (U_{\tau},Z_{\tau},\zeta_{\tau}) \sim U (0,1) \times S^{+} (\alpha,\rho) \times\overline{S} (\alpha,\rho) $ . Moreover, $ (U_{\tau},Z_{\tau},\zeta_{\tau}) $ is independent of $ (\tau,V_{\tau+1}) $ . Hence (2.5) will imply the perpetuity in the theorem if we prove that $\Lambda$ has the same law as $V_{\tau+1}$ . Put differently, as $\tau$ and $U_0$ are independent, it is sufficient to prove the following equality in law:

(2.6) \begin{equation} V_{\tau+1}\overset{d}{=}1_{\tau=-1}+1_{\tau\neq-1}U_{0}^{\frac{1}{\rho}}=1+1_{\tau\neq-1} (U_{0}^{\frac{1}{\rho}}-1) . \label{eqn7} \end{equation}

Since $-\log (1-U_{1}) \sim \operatorname{\rm{Exp}} (1) $ is exponential with mean one, $-\log (V_{n}) $ is gamma-distributed with density $x\mapsto x^{-n-1}\,\rm{\mathrm{e}}^{-x}/({-}n-1)!$ for any $n\in\mathcal{Z}^0$ . Hence, on the event $\{\tau\neq -1\}$ , the density of the conditional law $-\log (V_{\tau+1}) |\tau$ is given by $x\mapsto x^{-\tau-2}\,\rm{\mathrm{e}}^{-x}/({-}\tau-2)!$ . Thus, the conditional law $-\log (V_{\tau+1}) |\{\tau\neq -1\}$ is exponential with density

(2.7) \begin{equation} x\mapsto \dfrac{1}{1-\rho}\sum_{k=2}^{\infty}\rho (1-\rho) ^{k-1}\dfrac{x^{k-2}}{ (k-2) !}\,\rm{\mathrm{e}}^{-x}=\rho \,\rm{\mathrm{e}}^{-\rho x},\quad x\gt0. \label{eqn8} \end{equation}

Since $-\log (V_{\tau+1}) $ takes the value 0 when $\tau=-1$ , which happens with probability $\rho$ , and is otherwise exponential with mean $1/\rho$ , the distributional identity in (2.6) follows.

Finally, the uniqueness of the solution for perpetuity (1.1) follows from [Reference Buraczewski, Damek and Mikosch5, Theorem 2.1.3].

3. The Markov chain X and the dominating process D in Algorithm 1.1

Let $\mathcal{A}= (0,\infty) \times (0,1) \times (0,1) \times (0,1]$ and define the function $\phi\colon (0,\infty) \times\mathcal{A}\to (0,\infty) $ by

(3.1) \begin{equation} \phi (x,\theta) = \lambda^{\frac{1}{\alpha}} (u^{\frac{1}{\alpha}}x+ (1-u) ^{\frac{1}{\alpha}}s) ,\quad x\in (0,\infty) ,\quad\theta=(s,u,w,\lambda)\in\mathcal{A}. \label{eqn9} \end{equation}

Note that the map $x\mapsto \phi (x,\theta) $ is increasing and linear in x for all $\theta\in\mathcal{A}$ and does not depend on w. Let $W\sim U (0,1) $ be independent of random variables S, U, and $\Lambda$ defined in Theorem 1.1. Then, by Theorem 1.1, we have $\zeta\overset{d}{=}\phi (\zeta,\Theta) $ , where $\zeta\sim\overline{S} (\alpha,\rho) $ is independent of $\Theta= (S,U,W,\Lambda) $ . Hence a Markov chain with the update function $\phi$ has the correct invariant law but does not allow for coalescence: if for any $x,y\in(0,\infty)$ we have $\phi(x,\Theta)=\phi(y,\Theta)$ , by (3.1) it follows that $x=y$ . But the structure of $\phi$ and the additional randomness in W allow us to modify the update function $x \mapsto \phi(x,\theta)$ so that coalescence can be achieved, while keeping the law of the chain unchanged.

Lemma 3.1. Define the functions $\psi\colon (0,\infty) \times\mathcal{A}\to (0,\infty) $ and $a\colon \mathcal{A}\to (0,\infty) $ by the formulae

(3.2) \begin{align} \psi (x,\theta) & = 1_{ \{ a (\theta) \geq x\} }w^{\frac{1}{\alpha\rho}} (1-u) ^{\frac{1}{\alpha}}s+1_{ \{ a (\theta) \lt{x}\} }\lambda^{\frac{1}{\alpha}} (u^{\frac{1}{\alpha}}x+ (1-u) ^{\frac{1}{\alpha}}s) ,\label{eqn10}\end{align}

(3.3) \begin{align}* a (\theta) & = (\lambda^{-\frac{1}{\alpha}}-1) \bigg(\dfrac{1-u}{u}\bigg)^{\frac{1}{\alpha}}s. \label{eqn11}\end{align}

The map $x\mapsto\psi(x,\theta)$ is non-decreasing in x for all $\theta\in\mathcal{A}$ . Moreover, for $\zeta$ and $\Theta$ as in the paragraph above, we have $\phi (x,\Theta) \overset{d}{=}\psi (x,\Theta) $ for all $x\gt0$ and $\overline{S} (\alpha,\rho) $ is the unique solution of the distributional equation $\zeta\overset{d}{=}\psi (\zeta,\Theta) $ .

Proof. The function $\psi$ takes constant value of $w^{\frac{1}{\alpha\rho}} (1-u) ^{\frac{1}{\alpha}}s$ for $x\in (0,a (\theta) ]$ and increases linearly on the interval $ (a (\theta) ,\infty) $ with the right limit satisfying $\lim_{x\searrow a (\theta) }\psi (x,\theta) = (1-u) ^{\frac{1}{\alpha}}s>\psi (a (\theta) ,\theta) $ . Hence the desired monotonicity follows.

We now prove that $\phi (x,\Theta) \overset{d}{=}\psi (x,\Theta) $ for all $x>0$ , that is, the transition probabilities for the update functions $\phi$ and $\psi$ coincide. Pick $x>0$ and note that

\[ \{ \phi (x,\Theta) =\psi (x,\Theta) \} \supset \{ a (\Theta) <x\} . \]

Thus, for any $y>0$ we have

\[ \mathbb{P} ( \phi (x,\Theta) \leq y, a (\Theta) <x ) =\mathbb{P} ( \psi (x,\Theta) \leq y, a (\Theta) <x ) . \]

Define

\[ v (u,s) =\bigg(\dfrac{ (1-u) ^{\frac{1}{\alpha}}s}{u^{\frac{1}{\alpha}}x+ (1-u) ^{\frac{1}{\alpha}}s}\bigg)^{\alpha\rho}\in (0,1) , \]

and note that $ \{ a (\Theta) \geq x\} = \{ \Lambda^{\rho}\leq v (U,S) \} $ . On this event, the definition of $\Lambda$ in Theorem 1.1 implies the inequality $\Lambda<1$ , in which case $\Lambda^{\rho}$ is uniform on (0,1). Hence the conditional law of $\Lambda$ , given (U, S) and $\{ a (\Theta) \geq x\} $ , is uniform on the interval (0, v (U, S)). Moreover, the conditional law of v (U, S) W, given (U, S) and on $ \{ a (\Theta) \geq x\} $ , is also uniform on (0, v (U, S)). Hence, for any $y>0$ the following equalities hold:

\begin{align*} & \mathbb{P} ( \phi (x,\Theta) \leq y,\ a (\Theta) \geq x \!\mid\! U,S) \\ & \qquad = \mathbb{P} \bigg( \Lambda^{\rho}\leq\bigg(\dfrac{y}{U^{\frac{1}{\alpha}}x+ (1-U) ^{\frac{1}{\alpha}}S}\bigg)^{\alpha\rho},\ a (\Theta) \geq x \!\mid\! U,S\bigg) \\ & \qquad = \mathbb{P} \bigg( v (U,S) W\leq\bigg(\dfrac{y}{U^{\frac{1}{\alpha}}x+ (1-U) ^{\frac{1}{\alpha}}S}\bigg)^{\alpha\rho},\ a (\Theta) \geq x \!\mid\! U,S\bigg) \\ & \qquad = \mathbb{P} ( W^{\frac{p}{1}{\alpha\rho}} (1-U) ^{\frac{1}{\alpha}}S\leq y,\ a (\Theta) \geq x\ \!\mid\! U,S) \\ & \qquad = \mathbb{P} ( \psi (x,\Theta) \leq y,\ a (\Theta) \geq x \!\mid\! U,S) . \end{align*}

Taking expectations in this identity yields the unconditional equality

\[ \mathbb{P} (\phi (x,\Theta) \leq y, a (\Theta) \geq x ) =\mathbb{P} (\psi (x,\Theta) \leq y, a (\Theta) \geq x) . \]

Hence we get

\[ \mathbb{P} (\phi (x,\Theta) \leq y) = \mathbb{P} (\psi (x,\Theta) \leq y) \quad\text{for all $y\gt0$,} \]

implying the equality in law $\phi (x,\Theta) \overset{d}{=}\psi (x,\Theta) $ for arbitrary $x>0$ .

Pick $y\gt0$ . Since $\Theta$ and $\zeta$ are independent, by Theorem 1.1 we have

\begin{align*} \mathbb{P} (\zeta\leq y) & =\mathbb{P} (\phi (\zeta,\Theta) \leq y)\\ & = \int_{[0,\infty)}\mathbb{P} (\phi (x,\Theta) \leq y) \mathbb{P} (\zeta\in \rm{\mathrm{d}} x) \\ & = \int_{[0,\infty)}\mathbb{P} (\psi (x,\Theta) \leq y) \mathbb{P} (\zeta\in \rm{\mathrm{d}} x)\\ & =\mathbb{P} (\psi (\zeta,\Theta) \leq y) , \end{align*}

implying $\zeta\overset{d}{=}\psi (\zeta,\Theta) $ . Moreover, if there exists some $\zeta^{\prime}$ (independent of $\Theta$ ) satisfying $\zeta^{\prime}\overset{d}{=}\psi (\zeta^{\prime},\Theta) $ , this calculation implies the equality $\zeta^{\prime}\overset{d}{=}\phi (\zeta^{\prime},\Theta) $ . By Theorem 1.1 we get $\zeta^{\prime}\overset{d}{=}\zeta$ , as claimed.

By Lemma 3.1 and Kolmogorov’s consistency theorem, there exists a probability space supporting a sequence $ \{ \Theta_{n}\} _{n\in\mathbb{Z}}$ of independent copies of $\Theta$ and a stationary Markov chain $ \{ X_{n}\} _{n\in\mathbb{Z}}$ , satisfying $X_{n+1}=\psi (X_{n},\Theta_{n}) $ for all $n\in\mathbb{Z}$ . In the rest of the paper, $\{(X_n,\Theta_n)\}_{n\in\mathbb{Z}}$ denotes the corresponding Markov chain on $(0,\infty)\times\mathcal{A}$ . In order to detect coalescence in Algorithm 1.1, we now construct a dominating process $\{D_n\}_{n\in\mathbb{Z}}$ .

To this end, fix constants $\delta$ and d satisfying $0<\delta<d<\frac{p}{1}{\alpha\rho}$ . Let $I_{k}^{n}=1_{ \{ S_{k}>\rm{\mathrm{e}}^{\delta (n-1-k) }\} }$ for all $n\in\mathbb{Z}$ , $k\in\mathcal{Z}^{n}$ (see (2.3) above), where $S_k\sim S^{+}(\alpha,\rho)$ is the first component of $\Theta_k$ (see the first paragraph of Section 3). Fix $\gamma>0$ such that $\mathbb{E} S_{1}^{\gamma}<\infty$ (see (A.2)). Markov’s inequality implies

(3.4) \begin{equation} p (m) =\mathbb{P} (S_{1}\leq \rm{\mathrm{e}}^{\delta m}) \geq 1-\rm{\mathrm{e}}^{-\delta\gamma m}\mathbb{E} S_{1}^{\gamma},\quad m\geq0, \label{eqn12} \end{equation}

and hence $\sum_{m=0}^{\infty}(1-p (m) )<\infty$ . Since $\{S_k\}_{k\in\mathbb{Z}}$ are independent, the Borel–Cantelli lemma ensures that, for a fixed $n\in\mathbb{Z}$ , the events $ \{ S_{k}>\rm{\mathrm{e}}^{\delta (n-1-k) }\} = \{ I_{k}^{n}=1\} $ occur for only finitely many $k\in\mathcal{Z}^{n}$ a.s. Let $\chi_{n}$ be the smallest time beyond which the indicators $I_{k}^{n}$ are all zero:

(3.5) \begin{equation} \chi_{n}=(n-1)\wedge\inf \{ k\in\mathcal{Z}^{n}\colon I_{k}^{n}=1\} , \label{eqn13} \end{equation}

with convention $\inf\emptyset =\infty$ . Note that $-\infty\lt\chi_{n}\leq n-1$ holds a.s. for all $n\in\mathbb{Z}$ . Since the integers are countable, we have $n-1\geq \chi_{n}>-\infty$ for all $n\in\mathbb{Z}$ a.s.

Define the i.i.d. sequence $ \{ F_{n}\} _{n\in\mathbb{Z}}$ by $F_{n}=d+\pgfrac{1}{\alpha}\log (\Lambda_{n}U_{n}) $ , where $U_n$ and $\Lambda_n$ are the second and fourth components of $\Theta_n$ , respectively (see the first paragraph of Section 3). Note that $d-F_{n}$ has the same law as a sum of (random) geometrically many independent exponential random variables and is hence exponentially distributed with mean $\mathbb{E}[d-F_{n}]=\frac{p}{1}{\alpha\rho}$ . Let $C= \{ C_{n}\} _{n\in\mathbb{Z}}$ be a random walk defined by $C_{0}=0$ and

(3.6) \begin{equation} C_{n+1}=C_{n}-F_{n},\quad\text{$n\in\mathbb{Z}$.} \label{eqn14} \end{equation}

Recall definition (2.3) and let $R= \{ R_{n}\} _{n\in\mathbb{Z}}$ be the reflected process of the walk $ \{ C_{n}\} _{n\in\mathbb{Z}}$ , that is,

(3.7) \begin{equation} R_{n}=\sup_{k\in\mathcal{Z}^{n+1}}C_{k}-C_{n},\quad n\in\mathbb{Z}. \label{eqn15} \end{equation}

For any $n\in\mathbb{Z}$ , define the following random variables:

(3.8) \begin{align} D_{n} & = \exp (R_{n}) \bigg(\dfrac{\rm{\mathrm{e}}^{ (d-\delta) (\chi_{n}-n) }}{1-\rm{\mathrm{e}}^{\delta-d}}+ \sum_{k\in\mathcal{Z}^n_{\chi_{n}}}\,\rm{\mathrm{e}}^{- (n-1-k) d}S_{k} (1-U_{k}) ^{\frac{1}{\alpha}}\bigg),\label{eqn16}\end{align}

(3.9) \begin{align}* D_{n}^{\prime} & = \exp (R_{n}) \bigg(\dfrac{1}{1-\rm{\mathrm{e}}^{\delta-d}}+D_{n}^{\prime\prime}\bigg), \quad\text{where } D_{n}^{\prime\prime} = \sum_{k\in\mathcal{Z}^n}\,\rm{\mathrm{e}}^{- (n-1-k) d}S_{k}. \label{eqn17}\end{align}

The sum in (3.8) is taken to be zero if $\mathcal{Z}^n_{\chi_n}=\emptyset$ , i.e. if $\chi_n=n$ . Note that the series in $D_n''$ is absolutely convergent by the Borel–Cantelli lemma, but $D_n'$ cannot be simulated directly as it depends on an infinite sum. Finally, define the random element $\Xi_n= (\Theta_{n},R_{n},D_{n}^{\prime}) $ for any $n\in\mathbb{Z}$ .

Proof. (a) Since $\mathbb{E} F_{1}\lt0$ , by the strong law of large numbers we have $C_{-n}\to-\infty$ a.s. as $n\to\infty$ . Hence $R_{n}<\infty$ for all $n\in\mathbb{Z}$ a.s. and a direct termwise comparison yields $D_{n}^{\prime}\geq D_{n}$ for all $n\in\mathbb{Z}$ . It remains to prove that $X_{n}\leq D_{n}$ for all $n\in\mathbb{Z}$ .

Recall that the function $\theta\mapsto a(\theta)$ is defined in (3.3). Let $\tau_{n}=\sup \{ k\in\mathcal{Z}^{n}\colon X_{k}\leq a (\Theta_{k}) \} $ (with convention $\sup\emptyset = -\infty$ ) be the last time the coalescence occurred before $n\in\mathbb{Z}$ . If $\tau_{n}>-\infty$ , the value $X_{1+\tau_{n}}$ does not depend on $X_{\tau_{n}}$ , and neither do the values of the chain taken at subsequent times. In particular,

\[ X_{n}=\psi (X_{n-1},\Theta_{n-1}) =\underbrace{\psi(\cdots\psi}_{n-1-\tau_{n}}(W_{\tau_{n}}^{\frac{1}{\alpha\rho}} (1-U_{\tau_{n}}) ^{\frac{1}{\alpha}}S_{\tau_{n}},\Theta_{\tau_{n}+1}),\ldots,\Theta_{n-1}). \]

In general, by (3.2) and (2.3), $X_n$ can be expressed as

(3.10) \begin{align} X_{n} & = \sum_{k\in\mathcal{Z}^n_{\tau_{n}+1}}\exp\bigg(\dfrac{1}{\alpha} \sum_{j\in\mathcal{Z}^n_{k+1}}\log (\Lambda_{j}U_{j}) \bigg)\Lambda_{k}^{\frac{1}{\alpha}} (1-U_{k}) ^{\frac{1}{\alpha}}S_{k} \notag \\ & \quad\, + 1_{\{\tau_{n}>-\infty\}}\exp\bigg(\dfrac{1}{\alpha}\sum_{j\in\mathcal{Z}^n_{\tau_{n}+1}}\log (\Lambda_{j}U_{j}) \bigg)W_{\tau_{n}}^{\frac{1}{\alpha\rho}} (1-U_{\tau_{n}}) ^{\frac{1}{\alpha}}S_{\tau_{n}}, \label{eqn18}\end{align}

where sums over empty sets in (3.10) are defined to be equal to zero and, if $\tau_{n}=-\infty$ , we define $\mathcal{Z}^n_{\tau_{n}+1} =\mathcal{Z}^n$ . A termwise comparison then yields

(3.11) \begin{align} X_{n} & \leq \sum_{k\in\mathcal{Z}^n}\,\rm{\mathrm{e}}^{C_{k+1}-C_{n}- (n-1-k) d} (1-U_{k}) ^{\frac{1}{\alpha}}S_{k} \notag \\ & \leq \rm{\mathrm{e}}^{R_{n}}\sum_{k\in\mathcal{Z}^n}\,\rm{\mathrm{e}}^{- (n-1-k) d} (1-U_{k}) ^{\frac{1}{\alpha}}S_{k}\quad\text{for all \ensuremath{n\in\mathbb{Z}} a.s.} \label{eqn19}\end{align}

Recall that $S_{k} (1-I_{k}^{n}) \leq \rm{\mathrm{e}}^{\delta (n-1-k) } (1-I_{k}^{n}) $ for all $k\in\mathcal{Z}^n$ . Since $I_{k}^{n}=0$ for $k<\chi_{n}$ , we get

(3.12) \begin{align} & \sum_{k\in\mathcal{Z}^n}\,\rm{\mathrm{e}}^{- (n-1-k) d} (1-U_{k}) ^{\frac{1}{\alpha}}S_{k} \notag \\ &\qquad \leq \sum_{k\in\mathcal{Z}^{\chi_n}}\,\rm{\mathrm{e}}^{- (n-1-k) (d-\delta) } (1-U_{k}) ^{\frac{1}{\alpha}} +\sum_{k\in\mathcal{Z}^n_{\chi_{n}}}\,\rm{\mathrm{e}}^{- (n-1-k) d} (1-U_{k}) ^{\frac{1}{\alpha}}S_{k} \notag \\ & \qquad \leq \dfrac{\rm{\mathrm{e}}^{ (\chi_{n}-n) (d-\delta) }}{1-\rm{\mathrm{e}}^{\delta-d}}+\sum_{k\in\mathcal{Z}^n_{\chi_{n}}}\,\rm{\mathrm{e}}^{- (n-1-k) d} (1-U_{k}) ^{\frac{1}{\alpha}}S_{k}. \label{eqn20}\end{align}

The inequalities in (3.11)–(3.12) and the definition in (3.8) imply $X_{n}\leq D_{n}$ for all $n\in\mathbb{Z}$ a.s.

(b) Note that $C_k-C_n=\sum_{i=k}^{n-1}F_i$ for all $k\in\mathcal{Z}^n$ . Hence $R_n=\sup\{C_k-C_n\colon k\in\mathcal{Z}^{n+1}\}$ and $F_n$ are independent and the Markov property for $ \{ R_{n}\}_{n\in\mathbb{Z}}$ follows from

\[ R_{n}=\max \bigg\{ \sup_{k\in\mathcal{Z}^{n}}C_{k}-C_{n},0\bigg\} =\max \{ R_{n-1}+F_{n-1},0\} . \]

By (3.9) we have $D_{n}^{\prime\prime}=S_{n-1}+\rm{\mathrm{e}}^{-d}D_{n-1}^{\prime\prime}$ . Hence $ (R_{n},D_{n}^{\prime}) $ is a function of

\[ \Xi_{n-1}= (\Theta_{n-1},R_{n-1},D_{n-1}^{\prime}) \]

(recall that $S_{n-1}$ is the first component of the random vector $\Theta_{n-1}$ ). Since the random elements $\Xi_{n-1}$ and $\Theta_{n}$ are independent, the process $ \{\Xi_n\}_{n\in\mathbb{Z}}$ is Markov.

The vector $\Xi_n= (\Theta_{n},R_{n},D_{n}^{\prime}) $ is in a bijective correspondence with $ (\Theta_{n},R_{n},D_{n}^{\prime\prime}) $ .

Since $ \{ \Theta_{n}\}_{n\in\mathbb{Z}}$ are i.i.d., the following equality in law holds:

\begin{equation*} (R_{n+1},D_{n+1}^{\prime\prime}) = \bigg(\sup_{j\in\mathcal{Z}^1}\sum_{k\in\mathcal{Z}^1_j} F_{n+k} ,\sum_{k\in\mathcal{Z}^1}\,\rm{\mathrm{e}}^{kd}S_{n+k}\bigg) \overset{d}{=} \bigg(\sup_{j\in\mathcal{Z}^1}\sum_{k\in\mathcal{Z}^1_j} F_{k} ,\sum_{k\in\mathcal{Z}^1}\,\rm{\mathrm{e}}^{kd}S_{k}\bigg), \end{equation*}

implying the stationarity of $ \{ (\Theta_{n},R_{n},D_{n}^{\prime\prime}) \}_{n\in\mathbb{Z}}$ and hence of R and $\Xi$ .

The process R can jump to 0 in a single step and has positive jumps of size at most $1/(\alpha\rho)-d$ , both with positive probability. Hence it will hit any subinterval of its state space $[0,\infty)$ from any starting point in a finite number of steps with positive probability, making it $\varphi$ -irreducible [Reference Meyn and Tweedie24, p. 82] with respect to its invariant law.

Since $\Theta_n$ is independent of $(R_n,D_n'')$ , the $\varphi$ -irreducibility of $ \{\Xi_{n}\}_{n\in\mathbb{Z}}$ follows if, starting from an arbitrary point, we can prove that the process $\{(R_n,D_n'')\}_{n\in\mathbb{Z}}$ hits any rectangle in the product $[0,\infty)\times (0,\infty)$ with positive probability. Since we already know that R hits intervals and has (arbitrarily) small positive jumps with positive probability, the independence of $\{D''_n\}_{n\in\mathbb{Z}}$ and R, together with the fact that $D_n''$ has a positive density, imply the final statement of the lemma.

Proof of Theorem 1.2. By Lemma 3.2(ii), $\Xi$ is $\pi$ -irreducible, where $\pi$ denotes the invariant law of $\Xi$ . Hence, by [Reference Meyn and Tweedie24, Proposition 10.1.1], $\Xi$ is recurrent, meaning that the expected number of visits of the chain $\Xi$ to any set charged by $\pi$ is infinite for all starting points. By [Reference Meyn and Tweedie24, Theorem 9.0.1], the chain $\Xi$ is Harris-recurrent on a complement of a $\pi$ -null set. Put differently, for any starting point, the number of visits $\Xi$ makes to any set charged by $\pi$ is infinite almost surely.

Consider the Markov chain $\Psi=\{\Psi_n\}_{n\in\mathbb{N}}$ , where $\mathbb{N}=\{0,1,\ldots\}$ and $\Psi_n=\Xi_{-n}$ . In the language of [Reference Revuz28], $\Psi$ is a chain dual to $\Xi$ with respect to $\pi$ . In particular, the invariant law of $\Psi$ equals $\pi$ . Since $\Xi$ is Harris-recurrent on a state space with a countably generated $\sigma$ -algebra, [Reference Revuz28, Theorem 8.1.1] implies that there exists a modification of $\Psi$ (again denoted by $\Psi$ ) that is also Harris-recurrent. Since $\mathbb{P} (a (\Theta_{-n}) \geq D_{-n}^{\prime}) \gt{0}$ for any $n\in\mathbb{N}$ , it follows that the $\Psi$ -stopping time $\sigma^{\prime}=\inf \{ n>0\colon a (\Theta_{-n}) \geq D_{-n}^{\prime}\} $ is finite almost surely. Moreover, by [Reference Meyn and Tweedie24, Theorem 11.1.4] we have $\mathbb{E} [\sigma'|\Psi_0]<\infty$ almost surely.

Recall that $\sigma=\inf \{ n>0\colon a (\Theta_{-n}) \geq D_{-n}\} $ is the number of steps taken backwards in time in Algorithm 1.1. By Lemma 3.2(i) we have $\sigma\leq \sigma'$ . Since, by definition, $\Psi_0=\Xi_0$ , the claim follows.

4. Backward simulation of $ \{ (D_{n},\Theta_{n}) \} _{n\in\mathbb{Z}}$

A key step in Algorithm 1.1 consists of simulating the process $\{(D_n,\Theta_n)\}_{n\in\mathbb{Z}}$ backwards in time until the random time $\sigma=\inf \{ n>0\colon a(\Theta_{-n})\geq D_{-n}\} $ (see (2.3) and (3.3) for the definitions of $\mathcal{Z}^1$ and $a(\theta)$ , respectively). The forthcoming Algorithm 4.1 is responsible for this step. Recall that $\{\Theta_n\}_{n\in\mathbb{Z}}$ is an i.i.d. sequence with $\Theta_n=(S_n,U_n,W_n,\Lambda_n)$ having independent components, where $S_n$ , $U_n$ , and $\Lambda_n$ are distributed as in Theorem 1.1 and $W_n\sim U(0,1)$ .

At time $n\in\mathbb{Z}$ , the dominating process D in (3.8) depends on three components: the sequence $ (\chi_{n}, \{ S_{k}\} _{k\in\mathcal{Z}_{\chi_{n}}^{0}}) $ , the all-time maximum $\sup_{k\in\mathcal{Z}^{n+1}} \{C_k\}$ and $C_n$ (via the reflected process R: see (3.6)–(3.7)) and the uniform random variables $ \{ U_{k}\} _{k\in\mathcal{Z}_{\chi_{n}}^{0}}$ . The time $\chi_n$ in (3.5) is the last time before n the random variables $ \{ S_{k}\}_{k\in\mathcal{Z}^{0}}$ exceed a certain adaptive exponential bound. Algorithm 4.2 for sampling $ (\chi_{n}, \{ S_{k}\} _{k\in\mathcal{Z}_{\chi_{n}}^{0}}) $ is given in Section 4.1 below. A sample for $(R_n,C_n)$ requires the joint forward simulation of the dual random walk −C and its ultimate maximum. This problem was solved in [Reference Blanchet and Sigman4]. The algorithm in [Reference Blanchet and Sigman4], stated for completeness as Algorithm 4.6 of Section 4.2 below for the random walk C in (3.6), requires the simulation of the walk under the exponential change of measure.

Since the increments of C are shifted negative exponential random variables under the original measure, they remain in the same class under the exponential change of measure, making the simulation in Algorithm 4.6 simple. Finally, heaving simulated (R,C) backwards in time, we need to recover the random variables $\Lambda_{k}$ and $U_{k}$ , conditional on the values of increments $F_k=d+(1/\alpha) \log(U_n \Lambda_n)$ we have observed. Algorithm 4.7 in Section 4.3 below describes this step.

The number of steps $N_{-1}$ (resp. $N_{n}$ ) in line 2 (resp. 8) of Algorithm 4.1 is random since Algorithm 4.6, which outputs the all-time maximum of the random walk, may need more values of the random walk than required to recover the previous value of the dominating process $D_{-1}$ (resp. $D_{n}$ ). (In the notation of Section 4.2 below, the integers $N_{n}$ take the form $\Delta (\tau_{m}) $ .) The running time of Algorithm 4.2 is random but has moments of all orders (see Lemma 4.1 in Section 4.1 below). Algorithm 4.7 executes a loop of length equal to the number of steps in the random walk C the algorithm is applied to, with each step sampling one Poisson and one beta random variable (see Section 4.3 below). Hence both Algorithms 4.2 and 4.5 are fast (see Section 5). Algorithm 4.6 of [Reference Blanchet and Sigman4] (see Section 4.2 below) runs Algorithms 4.3, 4.4, and 4.5 sequentially. Each of these algorithms is reliant on rejection sampling and has a finite expected running time, which is easy to quantify in terms of the increments of the walk C.

4.1 Simulation of $ (\chi_{n}, \{ S_{k}\} _{k\in\mathcal{Z}_{\chi_{n}}^{0}}) $

Consider independent Bernoulli random variables $\{J_{n}\}_{n=1}^\infty$ with computable $p_{n}=\mathbb{P}(J_{n}=0)$ , $n\geq 1$ , satisfying $\sum_{n=1}^\infty(1-p_{n})<\infty$ . By the Borel–Cantelli lemma, the random time $\tau=\sup\{n\geq0\colon J_n=1\}^+$ (with convention $\sup\emptyset =-\infty$ ) satisfies $\tau\in\mathbb{N}$ a.s. Clearly, $J_{n}=0$ for all $n>\tau$ , and $\{\tau<n\} =\bigcap_{k=n}^{\infty}\{J_{k}=0\}$ implies $\mathbb{P}(\tau<n)=\prod_{k=n}^{\infty}p_{k}$ . If there exists $n^{\ast}\geq1$ such that for all $n\geq n^{\ast}$ we have a positive computable lower bound $q_{n}\leq\prod_{k=n}^{\infty}p_{k}$ , then we can simulate $(\tau,\{J_k\}_{k\in\{0,\ldots,\tau\}})$ as follows.

Define the auxiliary function $F\colon (0,1) \times (0,1) \to \{ 0,1\} \times (0,1) $ by the formula

\[ F (u,p) =\begin{cases} \bigg(0,\dfrac{u}{p}\bigg) & \text{ if }u\leq p,\\[9pt] \bigg(1,\dfrac{u-p}{1-p}\bigg) & \text{ if }u\gt{p}. \end{cases} \]

The following observation is simple but crucial: for any $p\in(0,1)$ and $U\sim U (0,1) $ , the components of the vector $(J,V) =F (U,p)$ are independent, J is Bernoulli with $\mathbb{P}(J=0)=p$ , and $V\sim U (0,1) $ .

Sample $ \{ J_{n}\} _{n\in\mathcal{Z}_1^{n^\ast}}$ and an independent $U^{ (n^{\ast}) }\sim U (0,1) $ . Let $ (J_{n^\ast},U^{ (n^\ast+1) }) =F (U^{ (n^\ast) },p_{n^\ast}) $ . Hence $J_{n^{\ast}}$ has the correct distribution and is independent of $U^{ (n^{\ast}+1) }\sim U(0,1)$ . Thus, $J_{n^{\ast}}$ is independent of $F (U^{ (n^{\ast}+1) },p_{n^{\ast}+1}) = (J_{n^{\ast}+1},U^{ (n^{\ast}+2) }) $ . Define recursively $ (J_{n},U^{ (n+1) }) =F (U^{ (n) },p_{n}) $ for $n\geq n^\ast+2$ and note that the sequence $ \{ J_{n}\}_{n\in\mathbb{N}}$ of Bernoulli random variables is i.i.d. Moreover, the sequence $\{U^{(n)}\}_{n\geq n^\ast}$ detects the value of $\tau$ since

\[ \{ U^{ (n) }\leq q_{n}\} \subseteq \bigg\{ U^{ (n) }\leq\prod_{k=n}^{\infty}p_{k}\bigg\} = \{ \tau<n\} . \]

Algorithm 4.2 (Simulation of $(\tau,\{J_k\}_{k\in\{1,\ldots,\tau\}})$ )

1. 1. Sample $J_1,\ldots, J_{n^\ast-1}$ and put $n\,:\!=n^\ast-1$

2. 2. Sample $U\sim U(0,1)$

3. 3. loop

4. 4. Put $n \,:\!= n+1$

5. 5. if $U\gt{p_n}$ then

6. 6. Put $J_n\,:\!=1$ and update $U\,:\!={(U-p_n)}/{(1-p_n)}$

7. 7. else if $U \leq q_n$ then

8. 8. Compute $\tau$ from $J_1,\ldots,J_{n-1}$ and exit loop

9. 9. else

10. 10. Put $J_n\,:\!=0$ and update $U\,:\!={U}/{p_n}$

11. 11. end if

12. 12. end loop

13. 13. return $(\tau,\{J_k\}_{k\in\{1,\ldots,\tau\}})$

Algorithm 4.2 samples a single uniform random variable and performs a binary search. Its running time $\varsigma=\inf \{ n\geq n^{\ast}\colon U^{ (n) }\leq q_{n}\} \geq\tau+1$ (with convention $\inf\emptyset =\infty$ ) has the following properties.

Lemma 4.1.

1. (a) If $\lim_{n\to\infty}q_{n}=1$ then $\mathbb{P}(\varsigma\lt\infty)=1$ .

2. (b) If $\sum_{n=n^{\ast}}^{\infty} (1-q_{n}) \lt\infty$ then $\mathbb{E} \varsigma\lt\infty$ .

3. (c) If $\sum_{n=n^{\ast}}^{\infty} (1-q_{n}) \,\rm{\mathrm{e}}^{tn}\lt\infty$ for some $t>0$ , then $\mathbb{E} \,\rm{\mathrm{e}}^{t\varsigma}\lt\infty$ .

   If $q_n p_{n-1}\geq q_{n-1}$ for $n>n^\ast$ , then the converses of (a), (b), and (c) are also true.

Remark 4.1. At the cost of more operations, one may always construct a sequence $\{q_n^\prime\}_{n=n^\ast}^\infty$ that satisfies (d). Indeed, let $q_{n^\ast}^\prime=q_{n^\ast}$ and define recursively $q_n^\prime = \max\{q_n,q^\prime_{n-1}/p_{n-1}\}$ for $n> n^\ast$ ; then these satisfy condition (d), are computable and inductively satisfy $q_{n}^\prime\leq \prod_{k=n}^\infty p_k$ for $n\geq n^\ast$ . This consideration shows that our conditions are sharp.

Proof. (a) For all $n\geq n^{\ast}$ we have $ \{ \varsigma\leq n\} \supseteq \{ U^{ (n) }\leq q_{n}\} $ ; then $\mathbb{P}(\varsigma>n)\leq\mathbb{P}(U^{(n)}>q_n)=1-q_{n}$ . Hence $\mathbb{P}(\varsigma=\infty)=\lim_{n\to\infty}\mathbb{P}(\varsigma>n) \leq\lim_{n\to\infty}(1-q_{n})=0$ and the sufficiency follows.

(b) Similarly,

\[ \mathbb{E} \varsigma=\sum_{n=0}^{\infty}\mathbb{P} (\varsigma>n) \leq n^{\ast}+\sum_{n=n^{\ast}}^{\infty} (1-q_{n}) <\infty \]

and the claim follows.

(c) Note that

\[ (\rm{\mathrm{e}}^{t}-1) \sum_{m=0}^{n-1}\,\rm{\mathrm{e}}^{tm}=\rm{\mathrm{e}}^{tn}-1. \]

Exchanging the order of summation in the third equality of the following estimate implies (c):

\begin{align*} \mathbb{E} \,\rm{\mathrm{e}}^{t\varsigma} & = \sum_{n=0}^{\infty}\mathbb{P} (\varsigma=n) \,\rm{\mathrm{e}}^{tn}\\ &= \sum_{n=0}^{\infty}\mathbb{P} (\varsigma=n) \bigg(1+ (\rm{\mathrm{e}}^{t}-1) \sum_{m=0}^{n-1}\,\rm{\mathrm{e}}^{tm}\bigg)\\ & = 1+ (\rm{\mathrm{e}}^{t}-1) \sum_{m=0}^{\infty}\,\rm{\mathrm{e}}^{tm}\mathbb{P} (\varsigma>m) \\ & \leq \rm{\mathrm{e}}^{tn^{\ast}}+ (\rm{\mathrm{e}}^{t}-1) \sum_{n=n^{\ast}}^{\infty} (1-q_{n}) \,\rm{\mathrm{e}}^{tn}\\ &<\infty. \end{align*}

(d) Condition (d) and the relation $(\tau+1)\vee n^\ast=\inf\{k\geq n^\ast\colon U^{(k)}\leq \prod_{j=k}^\infty p_j\}$ imply for $n\geq k\geq n^\ast$

\[ \{(\tau+1)\vee n^\ast = k, \varsigma\leq n\} = \begin{cases} \displaystyle \bigg\{ U^{(k-1)}\in \bigg[ p_{k-1},p_{k-1}+(1-p_{k-1})q_n \prod_{j\in\mathcal{Z}_{k}^n}p_j\bigg]\bigg\} & k>n^\ast,\\[9pt] \displaystyle\bigg\{ U^{(n^\ast)}\in \bigg[ 0,q_n\prod_{j\in\mathcal{Z}_{n^\ast}^n}p_j\bigg]\bigg\}& k=n^\ast. \end{cases} \]

Thus, a simple calculation yields

\[ \mathbb{P}(\varsigma \leq n) =q_n\prod_{j\in\mathcal{Z}_{n^\ast}^n}p_j+\sum_{k\in\mathcal{Z}_{n^\ast}^{n}}q_n(1-p_k)\prod_{j\in\mathcal{Z}_{k+1}^n}p_j =q_{n}, \]

and the result follows from standard probability theory.

In Algorithm 4.1 we are required to sample $ (\chi_{0}, \{ S_{k}\} _{\mathcal{Z}_{\chi_{0}}^{0}}) $ , and then, iteratively for $n\in\mathcal{Z}^0$ , $\chi_{n}$ and the remaining $ \{ S_{k}\} _{k\in\mathcal{Z}_{\chi_{n}}^{\chi_{n+1}}}$ , given the known values $ (\chi_{n+1}, \{ S_{k}\} _{k\in\mathcal{Z}_{\chi_{n+1}}^{0}}) $ . To apply Algorithm 4.2, we need a computable lower bound on the product of probabilities $p (m) =\mathbb{P}(S_1\leq \rm{\mathrm{e}}^{\delta m})$ , $m\in\mathbb{N}$ . Recall the exponential lower bound on p(m) in (3.4) and define

\[ m^{\ast}=\bigg\lfloor \dfrac{1}{\delta\gamma}\log\mathbb{E} S_1^{\gamma}\bigg\rfloor^+ +1 \]

(here $\lfloor x\rfloor=\sup\{n\in\mathbb{Z}\colon n\leq x\}$ for any $x\in\mathbb{R}$ ). Note that for any $m\geq m^{\ast}$ we have $\rm{\mathrm{e}}^{-\delta\gamma m}\mathbb{E} S_1^{\gamma}<1$ and may hence define

\[ \overline{p} (m) = \exp \bigg({-}\dfrac{1}{1-\rm{\mathrm{e}}^{-\delta\gamma}}\dfrac{\rm{\mathrm{e}}^{-\delta\gamma m} \mathbb{E} S_{1}^{\gamma}}{1-\rm{\mathrm{e}}^{-\delta\gamma m}\mathbb{E} S_{1}^{\gamma}}\bigg) \in (0,1). \]

The inequality in (3.4) implies

\begin{align*} \prod_{j=m}^{\infty}p (j) & \geq \prod_{j=m}^{\infty} (1-\rm{\mathrm{e}}^{-\delta\gamma j}\mathbb{E} S_{1}^{\gamma}) \\ & =\exp\bigg(\sum_{j=m}^{\infty}\log (1-\rm{\mathrm{e}}^{-\delta\gamma j}\mathbb{E} S_{1}^{\gamma}) \bigg)\\ & = \exp\bigg({-}\sum_{j=m}^{\infty}\sum_{k=1}^{\infty}\dfrac{1}{k}\,\rm{\mathrm{e}}^{-\delta\gamma jk} (\mathbb{E} S_{1}^{\gamma}) ^{k}\bigg)\\ & \geq\exp\bigg({-}\sum_{k=1}^{\infty}\dfrac{\rm{\mathrm{e}}^{-\delta\gamma mk} (\mathbb{E} S_{1}^{\gamma}) ^{k}}{1-\rm{\mathrm{e}}^{-\delta\gamma k}}\bigg)\\ & \geq \overline{p} (m) . \end{align*}

Since for any $k\in\mathcal{Z}^{0}_{\chi_0}$ we have

\[ \mathbb{P}(I^0_k=0) = \mathbb{P}(S_k\leq \rm{\mathrm{e}}^{-(k+1)\delta})=p({-}(k+1)), \]

Algorithm 4.2 can be applied (with $n^\ast=m^\ast$ ) to sample the sequence $ \{ I_{k}^{0}\}_{k\in\mathcal{Z}^{0}_{\chi_0}}$ . Moreover, for $m\in\mathbb{N}$ we get

\begin{equation*} \overline{p} (m^{\ast}+m) \geq \exp ({-}r\,\rm{\mathrm{e}}^{-\delta\gamma m}) \geq1-r\,\rm{\mathrm{e}}^{-\delta\gamma m}, \quad \text{where } r = \dfrac{\rm{\mathrm{e}}^{-\delta\gamma m^{\ast}}\mathbb{E} S_{1}^{\gamma}}{ (1-\rm{\mathrm{e}}^{-\delta\gamma}) (1-\rm{\mathrm{e}}^{-\delta\gamma m^{\ast}}\mathbb{E} S_{1}^{\gamma}) }\gt0. \end{equation*}

Hence, for any $t\in(0,\delta\gamma)$ , Lemma 4.1(c) implies that the running time $\varsigma$ satisfies $\mathbb{E} [\rm{\mathrm{e}}^{\varsigma t}]<\infty$ and therefore possesses moments of all orders. Having obtained $ (\chi_{0}, \{ I_{k}^{0}\}_{k\in\mathcal{Z}^{0}_{\chi_0}}) $ , for $k\in\mathcal{Z}^{0}_{\chi_0}$ , we sample $S_{k}$ as $S^{+} (\alpha,\rho) $ conditional on $S_{k}\leq \rm{\mathrm{e}}^{-\delta (k+1) }$ (if $I_{k}^{0}=0$ ) or $S_{k}>\rm{\mathrm{e}}^{-\delta (k+1) }$ (if $I_{k}^{0}=1$ ), yielding a sample of $ (\chi_{n}, \{ S_{k}\} _{k\in\mathcal{Z}_{\chi_{n}}^{0}}) $ .

Assume now that we have already sampled $ (\chi_{n+1}, \{ S_{k}\} _{k\in\mathcal{Z}_{\chi_{n+1}}^{0}}) $ . The adaptive exponential bounds in the indicators $I^{n+1}_k$ and $I^n_k$ are different (see Figure 1) and the relevant probabilities take the form

p^{\prime} (m) =\mathbb{P} (S_{1}\leq \rm{\mathrm{e}}^{\delta m} \!\mid\! S_{1}\leq \rm{\mathrm{e}}^{\delta (m+1) }) ,\quad m\in\mathbb{N}.

Since $\{S_{1}\leq \rm{\mathrm{e}}^{\delta m}\} \subset \{S_{1}\leq \rm{\mathrm{e}}^{\delta (m+1)}\}$ , the inequality $p^{\prime} (m) \geq p (m) $ holds for any $m\in\mathbb{N}$ . Thus

\begin{equation*} \prod_{j=m}^{\infty}p^{\prime} (j) \geq\prod_{j=m}^{\infty}p (j) \geq \overline{p} (m) {,} \end{equation*}

and Algorithm 4.2 can be applied with $n^\ast=\max\{m^\ast,n-\chi_{n+1}\}$ . The same argument as above shows that the running time $\varsigma$ has moments of all orders.

4.2 Simulation of the random walk and its reflected process from [Reference Blanchet and Sigman4]

In this section we present an overview of the algorithm in [Reference Blanchet and Sigman4] for the joint simulation of (C,R) defined in (3.6)–(3.7). We refer to [Reference Blanchet and Sigman4] and [Reference Ensor and Glynn17] for the proofs (the latter paper contains the simulation algorithm for the ultimate maximum of a random walk with negative drift and provides a basis for the simulation algorithm in [Reference Blanchet and Sigman4]).

Let $\eta=\eta(d)$ be the unique positive root of $\psi_{d}(\eta)=0$ , where $\psi_{d}(t)=\log(\mathbb{E} \,\rm{\mathrm{e}}^{tF_0})=dt-\log(1+t/(\alpha\rho))$ . Note that

\[ \psi_{d}^{\prime}(\eta)=d-\dfrac{1}{\alpha\rho+\eta}>0 \]

and $\eta =-\alpha\rho-W_{-1} ({-}\alpha\rho d\,\rm{\mathrm{e}}^{-\alpha\rho d}) /d$ , where $W_{-1}$ is the secondary branch of the Lambert W function. Since $\mathbb{E}[\exp(\eta F_n)]=1$ for all $n\in\mathbb{Z}$ , the process $\{ \exp(\eta C_{n})\}_{n\in\mathcal{Z}^{1}}$ is a positive backward martingale started at one, thus inducing a probability measure $\mathbb{P}^{\eta}$ on $\sigma$ -algebras $\sigma (C_{k};\,k\in\mathcal{Z}_{n}^{1}) $ , $n\in\mathcal{Z}^{1}$ , by the formula $\mathbb{P}^{\eta}(A)=\mathbb{E}[1_{A}\,\rm{\mathrm{e}}^{\eta C_{n}}]$ where $A\in\sigma (C_{k};\,k\in\mathcal{Z}_{n}^{1}) $ . Under $\mathbb{P}^{\eta}$ , the process C remains a random walk with i.i.d. increments satisfying $(\alpha\rho+\eta)(d-F_{n})\sim \operatorname{\rm{Exp}}(1)$ . Hence $\mathbb{E}^{\eta}[C_{-1}]=\psi_{d}^{\prime} (\eta) >0$ , implying $\lim_{n\to-\infty}C_{n}=\infty$ , $\mathbb{P}^\eta$ -a.s. by the strong law of large numbers.

For any $k\in\mathbb{Z}$ define (with convention $\sup \emptyset = -\infty$ )

(4.1) \begin{equation} T_{x}^{k}=\begin{cases} \sup \{ n\in\mathcal{Z}^{k}\colon C_{n}-C_{k}>x\} & \text{ if }x>0,\\ \sup \{ n\in\mathcal{Z}^{k}\colon C_{n}-C_{k}<x\} & \text{ if }x<0. \end{cases} \label{eqn21} \end{equation}

For ease of notation we let $T_{x}=T_{x}^{0}$ . Let E be an independent exponential random variable with mean one. Then, for $x>0$ , we have $\mathbb{P} (R_{0}>x) =\mathbb{P}^{\eta} (L_{E/\eta}>x) $ , where $L_{x}=\inf \{ y\geq0\colon C_{T_{y}}>x\} $ is the right inverse of $x\mapsto C_{T_x}$ ; see e.g. [Reference Ensor and Glynn17]. Hence, for $x\in(0,x')$ , where $x'\leq\infty$ , sampling $1_{ \{ R_{0}>x\} }=1_{ \{ T_{x}>-\infty\} }$ , conditional on $1_{ \{ R_{0}\leq x'\} }=1_{ \{ T_{x'}=-\infty\} }$ , in finite time amounts to sampling E and $C_{-1},\ldots,C_{T_{E/\eta}}$ under $\mathbb{P}^{\eta}$ ; see Algorithm 4.3 below.

Algorithm 4.3. (Simulation of $1_{\{R_{0}\gt{x}\} }$ conditional on $\{ R_{0}\leq x^{\prime}\} $ )

Require: $\infty \geq x^\prime >x >0$

1. 1. loop

2. 2. Sample $E\sim \operatorname{\rm{Exp}}(1)$

3. 3. if $E/\eta\leq x$ then

4. 4. return 0

5. 5. else

6. 6. Sample $C_0=0,C_{-1},\ldots,C_{T_{E/\eta}}$ under $\mathbb{P}^\eta$

7. 7. Compute $L_{E/\eta}$

8. 8. if $L_{E/\eta}\leq x^\prime$ then $\rhd$ Accept sample

9. 9. return $1_{\{L_{E/\eta}\gt{x}\}}$

10. 10. end if

11. 11. end if

12. 12. end loop

Remark 4.2. Since $L_x\leq x$ , then the condition $E/\eta\leq x$ implies $L_{E/\eta}\leq x$ , thus identifying $1_{\{L_{E/\eta}\gt{x}\}}=0$ (see line 3) and saving the computational effort of running all subsequent lines. This algorithm repeats independent experiments with success probability $\mathbb{P}^\eta(L_{E/\eta}\leq x^\prime)>0$ . The expected running time of each iteration in the loop is bounded above by $(\eta^{-1}+d)/\psi'_d(\eta)$ ; see [Reference Ensor and Glynn17, equation (2.3)]. Hence the expected running time of Algorithm 4.3 is finite.

In Algorithm 4.6 below we need to sample the path of the random walk $\{C_k\}_{k\in\mathcal{Z}^1_{T_x}}$ conditioned on the event $ \{R_{0}\in(x,x')\} $ , where $0<x<x'\leq\infty$ . By a rejection sampling method under $\mathbb{P}^\eta$ and Algorithm 4.3 (see [Reference Blanchet and Sigman4, Lemma 3]), this can be achieved as follows.

Algorithm 4.4 (Simulation of $C_{0},\ldots,C_{T_{x}}$ conditional on $\{ T_{x}\gt-\infty=T_{x^{\prime}}\} $ )

Require: $\infty\geq x^\prime >x >0$

1. 1. loop

2. 2. Sample $C_0=0, C_{-1}, \ldots, C_{T_{x}}$ under $\mathbb{P}^\eta$

3. 3. Given $C_{T_{x}}$ , sample independent $1_{\{R_0 ^\prime\leq x^\prime-C_{T_x}\}}$ and $U\sim U(0,1)$ $\rhd$ Algorithm 4.3

4. 4. if $U \leq \exp({-}\eta C_{T_{x}})$ and $1_{\{R_0 ^\prime\leq x^\prime-C_{T_x}\}}=1$ then

   $\rhd$ Accept sample

5. 5. return $\{C_n\}_{n\in\mathcal{Z}_{T_x}^1}$

6. 6. end if

7. 7. end loop

Remark 4.3. Since $L_x\leq x$ , we have $\mathbb{P}(R_0\leq z)\geq \mathbb{P}(E/\eta\leq z)=1-\exp({-}z\eta)$ for all $z\geq0$ . Since the overshoot $C_{T_x}-x$ is in the interval (0,d), the expected running time of Algorithm 4.4 (i.e. one over the acceptance probability) is smaller than $\exp(\eta (x+d))/(1-\exp({-}\eta( x'-x-d))$ if $x'\gt{x+d}$ .

In Algorithm 4.6 we also need to simulate the path of the walk reaching a negative level −x, while staying below a given positive level forever. Algorithm 4.5 achieves this (see [Reference Blanchet and Sigman4, Lemma 3]). Its expected running time is bounded above by

\[ 1/((1-\exp({-}\eta(x'+x)))\mathbb{P}(T_{-x}<T_{x^{\prime}}))<\infty. \]

Algorithm 4.5 (Simulation of $C_{0},\ldots,C_{T_{-x}}$ conditional on $ \{ T_{x^{\prime}}=-\infty\} $ )

Require: $x \in(0,\infty)$ and $x^\prime\in(0, \infty]$

1. 1. loop

2. 2. Sample $C_0=0, C_{-1}, \ldots, C_{T_{-x}}$ under $\mathbb{P}$

3. 3. Given $C_{T_{-x}}$ , sample an independent

   $1_{\{R_0 ^\prime\leq x^\prime-C_{T_{-x}}\}}$ $\rhd$ Algorithm 4.3

4. 4. if $1_{\{R_0 ^\prime\leq x^\prime-C_{T_{-x}}\}}=1$ and $\max_{n\in\mathcal{Z}_{T_{-x}}^1}\{C_n\} \leq x^\prime$ then

   $\rhd$ Accept sample

5. 5. return $\{C_n\}_{n\in\mathcal{Z}_{T_{-x}}^1}$

6. 6. end if

7. 7. end loop

We now give a brief overview of the algorithm in [Reference Blanchet and Sigman4] for the simulation of $ \{ (C_{n},R_{n}) \} _{n\in\mathcal{Z}^{1}}$ . Pick $\kappa\gt\max\{\log(2)/(3\eta),1/(\alpha\rho)\}$ (see assumption in [Reference Blanchet and Sigman4, Proposition 3]). Blanchet and Sigman [Reference Blanchet and Sigman4] constructed sequences $\Delta= \{ \Delta (k) \} _{k\geq0}$ and $\tau= \{ \tau_{k}\}_{k\geq0}$ of decreasing negative and increasing positive times, respectively.

1. (i) At the start of each iteration of the algorithm we are given

   \[ \Big( \{ \tau_{k}\} _{k\in \{ 0,\ldots,m\} },\ \{ \Delta (k) \} _{k\in \{ 0,\ldots,\tau_{m}\} },\ \{ C_{n}\}_{n\in\mathcal{Z}^1_{\Delta (\tau_{m}) }},\ \{ R_{n}\} _{n\in\mathcal{Z}^1_{\Delta (\tau_{m}-1) } }\Big) {.} \]

2. (ii) At each iteration we sample

   \[ \Big(\tau_{m+1}, \{ \Delta (k) \} _{k\in \{ \tau_{m}+1,\ldots,\tau_{m+1}\} },\ \{ C_{n}\} _{n\in\mathcal{Z}^{\Delta (\tau_{m}) }_{\Delta (\tau_{m+1}) }},\ \{ R_{n}\}_{n\in\mathcal{Z}^{\Delta (\tau_{m}-1) }_{\Delta (\tau_{m+1}-1) }} \Big). \]

Note that at the mth iteration we have $\Delta (\tau_{m}) -\Delta (\tau_{m}-1) $ more values of the walk than of the reflected process. More precisely, the algorithm starts by setting $\Delta(0)=0$ and repeats the following steps: given $ \{ \tau_{k}\} _{k\in \{ 0,\ldots,m\} }$ and $ \{ \Delta(k)\} _{k\in \{ 0,\ldots,\tau_{m}\} }$ , then put $\Delta (\tau_{m}+1) =T_{-2\kappa}^{\Delta (\tau_{m}) }$ . Next, if $\Delta (k) $ is the last known value of $\Delta$ and if $R_{\Delta (k) }\gt\kappa$ , then put $\Delta (k+1) =T_{\kappa}^{\Delta (k) }$ and $\Delta (k+2) =T_{-2\kappa}^{\Delta (k+1) }$ . If instead $R_{\Delta (k) }\leq\kappa$ , then put $\tau_{m+1}=k$ . Repeat the previous two steps until we can compute $\tau_{m+1}$ , that is, until $R_{\Delta (k) }\leq\kappa$ . After computing $\tau_{m+1}$ go back and repeat. By construction (see Proposition 3 in [Reference Blanchet and Sigman4]), we have

\[ \sup_{n\in\mathcal{Z}^{\Delta (\tau_{m}) +1}} \{ C_{n}\} \leq C_{\Delta (\tau_{m}-1) }-\kappa,\quad\text{implying } R_{n}=\max_{k\in\mathcal{Z}_{\Delta (\tau_{m}) +1}^{n+1}} \{ C_{k}\} -C_{n}, \ n\in\mathcal{Z}_{\Delta (\tau_{m}-1) }^{\Delta (\tau_{m-1}) }. \]

Hence, we may compute $R_{n}$ , $n\in\mathcal{Z}_{\Delta (\tau_{m}-1) }^{1}$ , from the simulated values

\[ \tau_{m},\quad \Delta_{ (\tau_{m}-1) },\quad \Delta_{ (\tau_{m}) },\quad \{ C_{n}\}_{n\in\mathcal{Z}_{\Delta (\tau_{m}) }^{1}}. \]

Algorithm 4.6. (Simulation of the random walk and its reflected process)

Require $\kappa \gt \max\{{\log(2)}/{(3\eta)}, {1}/{(\alpha\rho)}\}$ , $d\in (0,1)$ , $\infty\geq x\gt{0}$ and $m\geq1$

$\rhd$ x is an upper bound for $R_0$

1. 1. Put $t\,:\!=C_0\,:\!=\Delta(0)\,:\!=\tau_0\,:\!=0$

2. 2. for $k\in\{1,\ldots,m\}$ do

3. 3. Put $t\,:\!=\tau_{k-1}$

4. 4. loop

   Sample $C_{\Delta(t)-1}, \ldots, C_{T^{\Delta(t)}_{-2\kappa}}$ conditioned on $\{R_{\Delta(t)}<x-C_{\Delta(t)}\}$

   $\rhd$ Algorithm 4.5

5. 6. Put $\Delta(t+1)\,:\!=T^{\Delta(t)}_{-2\kappa}$ and $t\,:\!=t+1$

6. 7. Sample $1_{\{R_{\Delta(t)}\gt\kappa\}}$ given $\{R_{\Delta(t)} \lt x-C_{\Delta(t)}\}$

   $\rhd$ Algorithm 4.3

7. 8. if $1_{\{R_{\Delta(t)}>\kappa\}}=1$ then

8. 9. Sample $C_{\Delta(t)-1}, \ldots, C_{T^{\Delta(t)}_\kappa}$ from $\mathbb{P}^\eta$

   $\rhd$ Algorithm 4.4

9. 10. Put $\Delta(t+1)\,:\!=T^{\Delta(t)}_\kappa$ and $t\,:\!= t+1$

10. 11. else

11. 12. Put $x\,:\!=\kappa+C_{\Delta(t)}$ , $\tau_{k}\,:\!= t$ and exit loop

12. 13. end if

13. 14. end loop

14. 15. end for

15. 16. Compute $\{ R_{n}\}_{n\in\mathcal{Z}^1_{\Delta(\tau_{m}-1)} }$

16. 17. return $\Big(\{ \tau_{k}\}_{k\in\{ 0,\ldots,m\}}, \{\Delta(k)\}_{k\in\{ 0,\ldots,\tau_{m}\}}, \{C_{n}\}_{n\in\mathcal{Z}^1_{\Delta(\tau_{m})}}, \{R_{n}\}_{n\in\mathcal{Z}^1_{\Delta(\tau_{m}-1)} }\Big)$

4.3 Sampling $(U_n,\Lambda_n)$ given $F_n$

Algorithm 4.1 requires knowledge of $\{(U_n,\Lambda_n)\}_{n\in\mathcal{Z}^0}$ , given the increments $\{F_n\}_{n\in\mathcal{Z}^0}$ of the random walk C. Since $\log (U_{n}\Lambda_{n}) =\alpha (F_{n}-d) $ for all $n\in\mathbb{Z}$ , by independence, we may restrict attention to $n=1$ . It follows from (2.6) above that $\Lambda_{1}\overset{d}{=}\prod_{i=2}^{T}U_{i}$ for an independent geometric random variable T with parameter $\rho$ on the positive integers (if $T=1$ the right-hand side is defined to equal one). Hence, by independence, we have $U_1 \Lambda_{1}\overset{d}{=}\prod_{i=1}^{T}U_{i}$ . By (2.7), $-\log \Lambda_1$ conditioned on being positive is exponential with mean $1/\rho$ . Hence, for any $n\geq1$ and $y>0$ we obtain

\begin{align*} \mathbb{P}\bigg[ T=n\biggm|-\sum_{i=1}^{T}\log (U_{i} )=y \bigg] &=\gfrac{\bigg(\rho (1-\rho )^{n-1}\dfrac{y^{n-1}\,\rm{\mathrm{e}}^{-y}}{(n-1)!}\bigg)}{(\rho \,\rm{\mathrm{e}}^{-\rho y})}\\ &=\dfrac{[(1-\rho)y]^{n-1}\,\rm{\mathrm{e}}^{-(1-\rho)y}}{(n-1)!}. \end{align*}

Thus the conditional law of $T-1$ given $\sum_{i=1}^{T}\log (U_{i}) =-y$ is Poisson with mean $ (1-\rho) y$ . If $T=1$ , then $-\log(U_1)=y$ and $\Lambda_1=1$ . If $T>1$ , then for $x\in(0,y)$ we get

\begin{align*} \mathbb{P}\bigg[ -\log (U_{1}) \in \rm{\mathrm{d}} x\biggm|T=n,-\sum_{i=1}^{T}\log (U_{i}) =y\bigg] &=\bgfrac{\bigg(\rm{\mathrm{e}}^{-x}\dfrac{ (y-x) ^{n-2}\,\rm{\mathrm{e}}^{- (y-x) }}{ (n-2) !}\bigg)}{\bigg(\dfrac{y^{n-1}\,\rm{\mathrm{e}}^{-y}}{ (n-1) !}\bigg)}\,\rm{\mathrm{d}} x\\ & = (n-1) \dfrac{ (y-x) ^{n-2}}{y^{n-1}}\,\rm{\mathrm{d}} x. \end{align*}

Hence, conditional on $T=n$ and $\log (\prod_{i=1}^{T}U_{i}) =-y$ , the law of $-\pgfrac{1}{y}\log (U_{1}) $ is $\text{{beta}} (1,n-1) $ (understood as the Dirac measure $\delta_{1}$ when $n=1$ ). Finally we set $\Lambda_{1}=\exp (\alpha (F_1-d) ) /U_1$ .

Algorithm 4.7. (Simulation of $\{(U_{k},\Lambda_{k})\}_{k\in\mathcal{Z}_{m}^{n}}$ given $\{F_{k}\}_{k\in\mathcal{Z}_{m}^{n}}$ )

Require: $\{F_{k}\}_{k\in\mathcal{Z}_{m}^{n}}$ for $m,n\in\mathbb{Z}$ and $m< n$ .

1. 1. for $k\in \mathcal{Z}^n_m$ do

2. 2. Sample $T-1\sim \text{Poisson}({-}\alpha(F_{k}-d)(1-\rho))$

3. 3. Sample $L\sim \text{{beta}}(1,T-1)$

4. 4. Let $U_k\,:\!=\exp(L\alpha(F_k-d))$ and $\Lambda_k \,:\!= \exp((1-L)\alpha(F_{k}-d))$

5. 5. end for

6. 6. return $\{(U_{k},\Lambda_{k})\}_{k\in\mathcal{Z}_{m}^{n}}$

5. Implementation

Recall the definitions of the process $\{(C_n,F_n)\}_{n\in\mathbb{Z}}$ in 3.6, of $\{\Theta_n\}_{n\in\mathbb{Z}}$ in the first paragraph of Section 4 and of $\mathbb{P}^\eta$ in the second paragraph of Section 4.2. Before providing a concrete and concise algorithm and testing it, we will introduce a practical improvement based on a simple consideration.

Note that simulating the i.i.d. variables $\{\Theta_n\}_{n\in\mathcal{Z}^0}$ is clearly quicker and easier than employing the full machinery of our algorithms. Recall that the dominating process was introduced only to detect coalescence for the chain $\{X_n\}_{n\in\mathcal{Z}^0}$ . Thus, given $\{\Theta_n\}_{n\in\mathcal{Z}_{\Delta(0)}^0}$ for some burn-in parameter $\Delta(0)\in\mathcal{Z}^0$ and an upper bound $X^\prime_{\Delta(0)}=D_{\Delta(0)}\geq X_{\Delta(0)}$ (recall the definition of $\{D_n\}$ in (3.8)), one could recursively construct $X_{n+1}^\prime=\psi(X_n^\prime,\Theta_n)$ for $n\in\mathcal{Z}_{\Delta(0)}^{0}$ , and if any coalescence were detected, we would be certain that $X_0^\prime=X_0$ . Our objective is hence to take an appropriate $\Delta(0)$ that increases the probability $\mathbb{P}(X_0^\prime=X_0)$ . Algorithm 5.1 is a complete and compact simulation algorithm of $X_0$ , which makes use of this.

It is known that spectrally negative stable processes of infinite variation ( $\alpha>1$ and $\rho=1/\alpha$ ) satisfy $\overline{S}(\alpha,\rho)=S^+(\alpha,\rho)$ [Reference Michna25, Theorem 1]. As a simple application and sanity check, we now present a comparison between the empirical distribution function of $N=10^4$ samples and the actual distribution function in this case. To validate the samples, we compute the Kolmogorov–Smirnov statistic and test the hypothesis. (These graphs can be replicated following the guide available in [Reference González Cázares, Mijatovi and Uribe Bravo19].) In all three cases the null hypothesis of all samples coming from their respective distribution functions is not rejected (see Figure 2).

Figure 2: Empirical distribution functions for spectrally negative infinite variation stable process with parameters $\rho=1/\alpha$ and, from left to right, $\alpha=1.1$ , $\alpha=1.5$ , and $\alpha=1.9$ . The top graphs show the empirical distribution function $F_N$ for $N=10^4$ samples and compare it to the distribution function $F=\overline{S}(\alpha,\rho)=S^+(\alpha,\rho)$ [Reference Michna25, Theorem 1]. The bottom graphs show $s\mapsto\sqrt{N}(F_N\circ F^{-1}(s)-s)$ on [0,1] (or equivalently, the curve $t\mapsto(F(t),\sqrt{N}(F_N(t)-F(t)))$ in $\mathbb{R}^2$ for $t>0$ ), which converges weakly to a Brownian bridge. The dashed lines are the $0.05$ and 0.95 quantiles of the Kolmogorov–Smirnov statistic, derived from the distribution of the signed maximum modulus of the Brownian bridge.

5.1. Parameter choice and numerical performance

As explicitly stated in Algorithm 5.1, and if one allows m* (recall its definition in paragraph 1, p. 14) to vary over $\lfloor\pgfrac{1}{(\delta\gamma)}\log\mathbb{E} S_1^\gamma\rfloor^++\mathbb{N}$ , our simulation procedure has six different parameters. A full theoretical optimization is infeasible as it heavily depends on, among other things, the way the algorithm is coded, the computational cost of simulating each variable, the cost of each calculation, memory accessing cost, the quality and state of the RAM and $(\alpha,\rho)$ . However, for the sake of presenting its practical feasibility, we have implemented the algorithm in the Julia programming language (see [Reference González Cázares, Mijatovi and Uribe Bravo19]) and run it on a macOS Mojave 10.14.3 (18D109) with a 4.2 GHz Intel Core i7 processor and an 8 GB 2400 MHz DDR4 memory. This implementation is far from optimal, but still outputs 104 samples in approximately 1.15 seconds (without multi-threading) for the suggested parameters $(d,\delta,\gamma,\kappa,\Delta(0),m^\ast)=\varpi$ where

\begin{align*} \varpi & =\varpi(\alpha,\rho)\\ &=\bigg(\dfrac{2}{3\alpha\rho},\dfrac{1}{3\alpha\rho},\dfrac{19}{20}\alpha, 4+\max\bigg\{\dfrac{\log(2)}{3\eta(\frac{p}{2}{3\alpha\rho})},\dfrac{1}{\alpha\rho}\bigg\},40, 12+\bigg\lfloor\dfrac{60}{19}\rho\log\mathbb{E} S_1^{\frac{19}{20}\alpha}\bigg\rfloor^+\bigg). \end{align*}

This performance varies slightly for different choices of $(\alpha,\rho)$ . To put things in perspective, Algorithm 4.7 outputs, for the parameter choice $\varpi$ , 10⁶ samples in approximately 0.4322 seconds and drawing 106 samples from $S^+(\alpha,\rho)$ takes 0.1833 seconds. On the other hand, the first iteration of Algorithm 4.2 (which simulates the indicators $\{I_k^0\}_{k\in\mathcal{Z}^0_{\chi_0}}$ and the conditionally positive stable random variables $\{S_k\}_{k\in\mathcal{Z}^0_{\chi_0}}$ ) simulates 10⁴ samples in about 0.8125 seconds and is, although fast, the most computationally costly component of Algorithm 5.1. The main sources of this cost are the calculation of the probabilities {p(m)} (see their definition in (3.4)) and the simulation of $\{S_k\}_{k\in\mathcal{Z}^0_{\chi_0}}$ conditioned on the values of $\{I_k^0\}_{k\in\mathcal{Z}^0_{\chi_0}}$ .

Algorithm 5.1 (Perfect simulation of $X_0\overset{d}{=}\overline{Y}_1$ )

Require: Parameters

\[0\lt\delta\ltd\lt\dfrac{1}{\alpha\rho},\quad \kappa\gt\max\bigg\{\dfrac{\log(2)}{3\eta},\dfrac{1}{\alpha\rho}\bigg\},\quad \gamma\gt0{,}\quad \Delta(0)\in\mathcal{Z}^0 \]

1. 1. Put $x\,:\!=\infty$ , $t\,:\!=1$ ,

   $s\,:\!=\Delta(0)$ , and $m\,:\!=n\,:\!=\Delta(0)+1$

   $\rhd$ x is an upper bound on $\{C_k\}_{k\in\mathcal{Z}^{\Delta(0)}}$

2. 2. Sample $\{\Theta_k\}_{k\in\mathcal{Z}^0_{\Delta(0)}}$

   $\rhd$ Recall its definition in Section 4, paragraph 1

3. 3. loop

4. 4. Sample $(\chi_{m-1}, \{S_k\}_{k\in\mathcal{Z}_{\chi_{m-1}}^{s}})$ $\rhd$ Algorithm 4.2

5. 5. while $n = m$ or $\Delta(t) \gt \chi_{m-1}$ do

6. 6. Sample $C_{\Delta(t)}, \ldots, C_{T^{\Delta(t)}_{-2\kappa}}$ conditional on $\{R_{\Delta(t)}<x-C_{\Delta(t)}\}$

   $\rhd$ Algorithm 4.5

7. 7. Put $\Delta(t+1)\,:\!=T^{\Delta(t)}_{-2\kappa}$ and $t\,:\!=t+1$

   $\rhd$ Recall its definition in (4.1)

8. 8. Sample $1_{\{R_{\Delta(t)}\gt\kappa\}}$ given $\{R_{\Delta(t)} \lt x-C_{\Delta(t)}\}$

   $\rhd$ Algorithm 4.3

9. 9. if $1_{\{R_{\Delta(t)}\gt\kappa\}}=0$ then

10. 10. Compute $\{ R_{k}\}_{k\in\mathcal{Z}^n_{\Delta(t-1)}}$ and put $n\,:\!=\Delta(t-1)$ and $x\,:\!=C_{\Delta(t)}+\kappa$

11. 11. else

12. 12. Sample $C_{\Delta(t)-1}, \ldots, C_{T_\kappa^{\Delta(t)}}$ from $\mathbb{P}^\eta$

    $\rhd$ Algorithm 4.4

13. 13. Put $\Delta(t+1)\,:\!=T_\kappa^{\Delta(t)}$ and $t\,:\!= t+1$

14. 14. end if

15. 15. end while

16. 16. Sample $\{(U_k,\Lambda_k)\}_{k\in\mathcal{Z}_{\chi_{m-1}}^{s}}$ from $\{F_{k}\}_{k\in\mathcal{Z}_{\chi_{m -1}}^{s}}$ and put $s\,:\!=\chi_{m-1}$ $\rhd$ Algorithm 4.7

17. 17. Compute $D_{m-1}$ and put $m\,:\!=m-1$

    $\rhd$ Recall its definition in (3.8)

18. 18. if $D_{m}\leq a(\Theta_{m})$ then

19. 19. return $X_0\,:\!=\psi(\cdots\psi(D_{m},\Theta_{m}),\ldots,\Theta_{-1})$ $\rhd$

    In this case $\sigma=m$

20. 20. else if $m=\Delta(0)$ then

21. 21. Put $X_0\,:\!=\psi(\cdots\psi(D_m,\Theta_m),\ldots,\Theta_{-1})$

22. 22. if coalescence was detected then

23. 23. return $X_0$

24. 24. end if

25. 25. end if

26. 26. end loop

Next we show the local marginal behaviour of the number of samples output with confidence intervals, for a few different choices of parameters $(\alpha,\rho)$ . We will see that, although $\varpi$ may not be optimal, it is a simple and yet efficient choice. Moreover, the variation in performance for parameters close to this one is small, thus showing that this choice is relatively robust.

It should be noted that the data presented in Figure 3 are dependent on the characteristics of the hardware and software used. Hence, these exact numbers are not easily replicated. For instance, these times scale sublinearly as a function of the batch size, and are not replicated despite using the garbage collector and the same random seed. It is readily seen that the exact value of the parameters d, $\delta$ , $\Delta(0)$ , and m* is not too important in so far as they remain at a reasonable distance from their boundaries (where $\infty$ is a right-boundary for $\Delta(0)$ and m*). The value of $\kappa$ is slightly more sensitive, as is $\gamma$ . Other choices of $(\alpha,\rho)$ have slightly different behaviours. The shapes of these curves are similar, but the apparent minima change. Thus, we argue that $\varpi$ is a simple yet sensible choice.

Figure 3: Time taken (in seconds) to simulate $N=10^4$ samples of $\overline{S}(1.3,1/2)$ with parameters moving about $\varpi$ . In each plot, one parameter moves and all others are kept constant at the respective value of $\varpi$ . We took 100 batches of samples, each with $N=10^4$ independent simulations, to construct asymptotic 95% confidence intervals based on the central limit theorem.