Fast exact simulation of the first-passage event of a subordinator

TL;DR

This work develops an exact simulation algorithm for the first-passage event of a general subordinator across a non-increasing boundary, jointly sampling the crossing time and the undershoot/overshoot. The Lévy measure is decomposed as $\nu_Z = \nu_{r,q} + \lambda_r$, enabling a two-component simulation: a truncated tempered-stable part $Y$ and a finite-activity Poisson part $Q$, with the crossing event analyzed via a barrier $b(t) = \min\{r\rho, c(t)\}$. The authors prove that the algorithm’s running time has moments of all orders and provide an explicit bound on the expected running time in terms of model parameters, including $\alpha$, $\vartheta$, $q$, $c_0$, $r$, and $\rho$, together with the costs of the underlying tempered-stable first-passage sampling. The paper also develops efficient sampling of tempered-stable variables conditioned to be small (via Zolotarev representations and Devroye’s method) and discusses practical applications, notably Monte Carlo estimation for FPDEs, with guidance on choosing the truncation parameter and an available Python/Julia implementation. Overall, the approach offers finite-mean running-time guarantees and broad applicability to Monte Carlo schemes requiring exact first-passage sampling for subordinators.

Abstract

This paper provides an exact simulation algorithm for the sampling from the joint law of the first-passage time, the undershoot and the overshoot of a subordinator crossing a non-increasing boundary. We prove that the running time of this algorithm has finite moments of all positive orders and give an explicit bound on the expected running time in terms of the Lévy measure of the subordinator. This bound provides performance guarantees that make our algorithm suitable for Monte Carlo estimation. We provide a GitHub repository with an implementation of the algorithm in Python and Julia.

Figures (6)

• Figure 1.1: Flowchart of Algorithm \ref{['alg:triplet_Z']} for the simulation of the random element $(\tau_c^Z, Z_{\tau_c^Z-}, Z_{\tau_c^Z})$, where the subordinator $Z$ is defined in \ref{['eq:Levy_measure_description']}. The pair $(D,J)$ is the first jump time and size of the compound Poisson process defined via the Lévy measure $\lambda_r$ in \ref{['eq:Levy_measure_description']}, process $Y$ is the driftless truncated tempered stable subordinator with Lévy measure $\nu_{r,q}$ given in \ref{['eq:Levy_measure_description']}, introduced at the beginning of Section \ref{['section:result']}. There are two key steps in Algorithm \ref{['alg:triplet_Z']}: (i) the simulation of the triplet $\left( \tau _{b}^{Y} ,Y_{\tau _{b}^{Y} -} ,Y_{\tau _{b}^{Y}}\right)$ (see line \ref{['step:fpe_truncated_tempered_stable']} of Algorithm \ref{['alg:triplet_Z']}), based on cazares2023fast for the tempered stable first-passage event, and (ii) the simulation from the conditional law $\mathcal{L}( Y_D |\{Y_D < b( D)\})$ via Algorithm \ref{['alg:small_tempered_stable']} below, which is based on Devroye's algorithm devroye2012note for log-concave densities. Algorithm \ref{['alg:small_tempered_stable']} is critical for the bound on the expected running time in Corollary \ref{['cor:time_complexity']} below, because $b(D)$ might be tiny with significant probability, making a naive simulation of $Y_D$, until $Y_D<b(D)$, possibly have infinite expected running time. The parameter $\rho\in(0,1)$ can be optimised in terms of $\alpha$. Its optimal value lies in the interval $[1/\mathrm{e},1/2]$ and we typically set it to $1/2$.
• Figure 1.2: The graph plots the complexity of the implementation in repository of Algorithm \ref{['alg:triplet_Z']} with parameters $\vartheta=2$, $r_0=\infty$, $c(t)\equiv 5$, $\lambda_{r_0}(\mathrm{d} x)=\mathrm{e}^{-x}\mathrm{d} x$, $q=10$ and the truncation $r=\min\{r_0,2\alpha/q\}=2\alpha/q$ for the range $\alpha\in[.05,.95]$. The average running time $\mathbb{E}[\mathcal{T}]$ is measured in seconds taken for every $10^4$ samples.
• Figure 2.1: In (i), $Y$ (and hence $Z$) crosses the boundary $b(t):=\min\{c(t),r\rho\}$ before the first jump of $Q$ at time $D$. Note that in (ii), on the event $\{D<\tau_b^S\}$, we have $Y_D<b(D)$ but it is possible to have either $Y_D+J>b(D)$ or $Y_D+J\leqslant b(D)$ ($J$ is the magnitude of the jump of $Q$ at time $D$).
• Figure 4.1: The pictures show an implementation of Algorithm \ref{['alg:triplet_Z']} with parameters $\vartheta=2$, $r_0=\infty$, $c(t)\equiv 5$, $\lambda_{r_0}(\mathrm{d} x)=\mathrm{e}^{-x}\mathrm{d} x$. The average running time $\mathbb{E}[\mathcal{T}]$ is measured in seconds taken for every $10^4$ samples.
• Figure 4.2: A Monte Carlo estimator $u_n(t,x)$ of the solution $u(t,x)$ of the FPDE in \ref{['eq:fpde_solution']} at time $t-a=5$, based on the representation in \ref{['eq:fpde_solution']}. This computation is based on $n=10^4$ samples.

Theorems & Definitions (18)

• Theorem 2.1
• Corollary 2.2
• proof : Proof of Corollary \ref{['cor:time_complexity']}
• Proposition 3.1
• Remark
• proof : Proof of proposition \ref{['prop:small_temperd_stable']}
• Lemma 5.1
• proof
• Lemma 5.2
• proof