Fast exact simulation of the first passage of a tempered stable subordinator across a non-increasing function

TL;DR

This work provides a fast, exact simulation method for the first-passage time of a tempered stable subordinator across a non-increasing boundary, yielding the triplet $(\tau_b,S_{\tau_b-},S_{\tau_b})$ with running times possessing finite exponential moments. The core innovation is reducing the tempered problem to the stable case and overcoming the challenging undershoot sampling via a domination by a mixture of tractable densities, employing Devroye’s log-concave sampling and careful Newton–Raphson inversions. The authors establish explicit, parameter-dependent bounds on the expected running time and its exponential moments, and show how to achieve near-linear-tempered-time behavior by iteratively capping the boundary; they also provide comprehensive theoretical analyses and practical implementations (in Julia) with applications to barrier options and fractional PDEs. The results advance Monte Carlo methods for non-local operators and tempered fractional derivatives by enabling practically exact, efficiently computable crossing-time samplings. Collectively, the methodology supports robust MC estimation in FPDEs and finance while clarifying the computational trade-offs across the stability parameter, tempering, and boundary behavior.

Abstract

We construct a fast exact algorithm for the simulation of the first-passage time, jointly with the undershoot and overshoot, of a tempered stable subordinator over an arbitrary non-increasing absolutely continuous function. We prove that the running time of our algorithm has finite exponential moments and provide bounds on its expected running time with explicit dependence on the characteristics of the process and the initial value of the function. The expected running time grows at most cubically in the stability parameter (as it approaches either $0$ or $1$) and is linear in the tempering parameter and the initial value of the function. Numerical performance, based on the implementation in the dedicated GitHub repository, exhibits a good agreement with our theoretical bounds. We provide numerical examples to illustrate the performance of our algorithm in Monte Carlo estimation.

Figures (8)

• Figure 2.1: Dependence of algorithms ($B\leftarrow A$ means algorithm $A$ calls algorithm $B$)
• Figure 2.2: Density of the second marginal of $\widetilde{\psi}_s^{(2)}$. The figure shows the graphs of the density of the second marginal of $\widetilde{\psi}_s^{(2)}$, in a custom logarithmic-style scale, for two processes with $\alpha=0.85$ and $\alpha=0.9$, where the parameter $s$ is at its $99^{\text{th}}$ quantile. Thus, on average in one of every $100$ runs of \ref{['alg:triple_stable_conditional_on_time']}, \ref{['alg:psi_2']} is called to sample from a density which, at the mode $z=\sigma_\alpha^{-1}(\alpha s^r)$, takes extremely large value followed by a steep drop towards zero. Moreover, in this case the density is strictly increasing on the interval $(0,z)$ and strictly positive on $(0,1)$. But since its derivatives of all orders are equal to zero at $1$, the graphs appear to be equal to zero on most of the interval $(z,1)$. Furthermore, the behaviour of the density on either side of the mode $z$ is markedly different. Taking a uniform upper bound $\mathrm{e}^{-\alpha}\alpha^\alpha s^{\alpha^2/(1-\alpha)}$ would result in an extremely inefficient algorithm with infinite expected running time (in fact, such a bound was used to sample from this type of function in chi2016exact, see Subsection \ref{['subsec:Chi']} below). For parameter values as in the graph above, \ref{['alg:psi_2']} sets $z_*=1/2$, splits the interval $(0,1)=(0,z_*]\cup(z_*,z]\cup(z,1)$ and uses appropriate proposal densities on each subinterval.
• Figure 2.3: Numerical inversion with numerical integrals. The graph on the left (resp. right) show the function $y\mapsto \sqrt{n}(F_n(y)-G_n(y))$ where $F_n$ and $G_n$ are the empirical distribution functions of the samples produced using \ref{['alg:psi_1']} (resp. \ref{['alg:psi_2']}) and direct numerical inversion, respectively, for $n=10^4$. Such functions indeed resemble Brownian bridges time-changed by the corresponding distribution functions. All samples pass the two-sample Kolmogorov-Smirnov test, confirming that our algorithms and direct numerical inversion produce samples from the same laws.
• Figure 3.1: Implementation of \ref{['alg:triple_stable_conditional_on_time']} with $t_*=\infty$ on Julia. The average time $\mathbb{E}[T]$ is measured in seconds taken for every $10^4$ samples.
• Figure 3.2: Implementation of \ref{['alg:improved_triple_temper_stable']}. The picture depicts an implementation of \ref{['alg:improved_triple_temper_stable']} with parameters $\alpha=0.55$, $\theta=1$, $b(t)\equiv 1$ and $q\in\{\mathrm{e}^0, \mathrm{e}^{1/4},\ldots,\mathrm{e}^{11}\}$. The mean running time $\mathbb{E}[T]$ is measured in seconds taken for every $10^2$ samples. By Corollary \ref{['cor:improved_tempered_expected_time']}, $\mathbb{E}[T]=\mathcal{O}(q)$ as $q\to\infty$. This is supported by the picture: the quotient $\log(\mathbb{E} [T])/\log(\mathrm{e}+q)$ appears to converge to $1$ as $q\to\infty$.

Theorems & Definitions (46)

• Theorem 2.1
• Corollary 2.2
• Theorem 2.3
• Remark 2.1
• Remark 2.2
• Proposition 2.4
• Proposition 2.5
• Remark 2.3
• Proposition 2.6
• Proposition 4.1