Asymptotic Exponentiality of the First Exit Time of the Shiryaev-Roberts Diffusion with Constant Positive Drift

TL;DR

This paper proves that the first exit time of a Shiryaev-Roberts diffusion with positive drift becomes exponentially distributed as the boundary goes to infinity, with an explicit analytical approach.

Contribution

It provides an explicit analytical proof that the standardized first exit time converges to an exponential distribution as the boundary size increases.

Findings

The moment generating function converges to that of an exponential distribution.

Explicit analytical expression for the first exit time's MGF is derived.

The result extends understanding of quickest change-point detection methods.

Abstract

We consider the first exit time of a Shiryaev-Roberts diffusion with constant positive drift from the interval $[0,A]$ where $A>0$. We show that the moment generating function (Laplace transform) of a suitably standardized version of the first exit time converges to that of the unit-mean exponential distribution as $A\to+\infty$. The proof is explicit in that the moment generating function of the first exit time is first expressed analytically and in a closed form, and then the desired limit as $A\to+\infty$ is evaluated directly. The result is of importance in the area of quickest change-point detection, and its discrete-time counterpart has been previously established - although in a different manner - by Pollak and Tartakovsky (2009).

Equations

$$d{R}_t^r$$

$${\mathcal{S}}_A^r$$

$$M(\alpha ;A,x)$$

$$\mathcal{D}$$

$$\displaystyle\dfrac{\mu^{2}}{2}\dfrac{d^{2}}{dx^{2}}\big{[}x^{2}\,u(x,\lambda)\big{]}-\dfrac{d}{dx}\big{[}u(x,\lambda)\big{]}$$

$$\displaystyle\lim_{x\to 0+}\left\{\dfrac{\mu^{2}}{2}\dfrac{\partial}{\partial x}\big{[}x^{2}\,u(x,\lambda)\big{]}-u(x,\lambda)\right\}$$

$$E\left\{e^{-\alpha {\mathcal{T}}_y^x}\right\}$$

$$\displaystyle\dfrac{1}{2}b(x)\dfrac{\partial^{2}}{\partial x^{2}}\big{[}v(x;\alpha)\big{]}+a(x)\dfrac{\partial}{\partial x}\big{[}v(x;\alpha)\big{]}$$

$$\psi (x,\lambda )$$

$$\xi$$

$$\frac{{\partial }^2}{\partial z^2}w(z)+\left\{-\frac{1}{4}+\frac{a}{z}+\frac{1/4-b^2}{z^2}\right\}w(z)$$

$$M(\alpha ;A,x)$$

$$W_{a,b}(x)$$

$$\frac{{\mu }^{2}A}{2}{e}^{\frac{1}{{\mu }^{2}A}}{W}_{1,\frac{1}{2}\xi (\alpha )}\left(\frac{2}{{\mu }^{2}A}\right)$$

$$\underset{A\to +\infty }{\lim }\left\{\frac{{\mu }^{2}A}{2}{e}^{\frac{1}{{\mu }^{2}A}}{W}_{1,\frac{1}{2}\xi (\alpha )}\left(\frac{2}{{\mu }^{2}A}\right)\right\}$$

$$Q_A(x)$$

$$W_{1,\frac{1}{2}\xi ({\lambda }_A)}\left(\frac{2}{{\mu }^{2}A}\right)$$

$$\underset{A\to +\infty }{\lim }\left\{\frac{\frac{{\mu }^{2}x}{2}{e}^{\frac{1}{{\mu }^{2}x}}{W}_{1,\frac{1}{2}\xi (-\alpha {\lambda }_A)}\left(\frac{2}{{\mu }^{2}x}\right)}{\frac{{\mu }^{2}A}{2}{e}^{\frac{1}{{\mu }^{2}A}}{W}_{1,\frac{1}{2}\xi (-\alpha {\lambda }_A)}\left(\frac{2}{{\mu }^{2}A}\right)}\right\}$$

$$\underset{A\to +\infty }{\lim }\left\{\frac{{\mu }^{2}x}{2}{e}^{\frac{1}{{\mu }^{2}x}}{W}_{1,\frac{1}{2}\xi (-\alpha {\lambda }_A)}\left(\frac{2}{{\mu }^{2}x}\right)\right\}$$

$$\underset{A\to +\infty }{\lim }\left\{\frac{{\mu }^{2}A}{2}{e}^{\frac{1}{{\mu }^{2}A}}{W}_{1,\frac{1}{2}\xi (-\alpha {\lambda }_A)}\left(\frac{2}{{\mu }^{2}A}\right)\right\}$$

$${\mathrm{E}}_{1}1(x)$$

$$G_{n,p}^{\ast ,m}(\begin{array}{c}q\\ a_1,\dots ,a_p\end{array}|b_1,\dots ,b_q)z$$

$$G_{1,2}^{\ast ,3}(\begin{array}{c}3\\ 0,1\end{array}|0,0,0)x$$

$$\displaystyle\begin{split}W_{1,\tfrac{1}{2}\xi(\lambda)}&\left(\dfrac{2}{\mu^{2}x}\right)=\dfrac{2}{\mu^{2}}\,e^{-\tfrac{1}{\mu^{2}x}}\Biggl{\{}\dfrac{1}{x}+\lambda+\dfrac{2}{\mu^{2}}L\left(\dfrac{2}{\mu^{2}x}\right)\lambda^{2}+\\ &\quad+\left(\dfrac{2}{\mu^{2}}\right)^{2}\left[G^{\,*,3}_{1,2}\lparen\begin{smallmatrix}3\\ 0,1\end{smallmatrix}|\,0,0,0\rparen{\dfrac{2}{\mu^{2}x}}-2L\left(\dfrac{2}{\mu^{2}x}\right)\right]\lambda^{3}\Biggr{\}}+\mathcal{O}(\lambda^{4}),\end{split}$$

$$L(x)$$

$$\begin{array}{rl}\frac{{\mu }^{2}A}{2}{e}^{\frac{1}{{\mu }^{2}A}}{W}_{1,\frac{1}{2}\xi (-\alpha {\lambda }_A)}& \left(\frac{2}{{\mu }^{2}A}\right)=1-A(\alpha {\lambda }_A)+A\frac{2}{{\mu }^2}L\left(\frac{2}{{\mu }^{2}A}\right)(\alpha {\lambda }_A{)}^2+\\ & -A{\left(\frac{2}{{\mu }^2}\right)}^2\left[G_{1,2}^{\ast ,3}(\begin{array}{c}3\\ 0,1\end{array}|0,0,0)\frac{2}{{\mu }^{2}A}-2L\left(\frac{2}{{\mu }^{2}A}\right)\right](\alpha {\lambda }_A{)}^3+A{\alpha }^4\mathcal{O}({\lambda }_A^4),\end{array}$$

$$-\frac{1}{A}$$

$$\underset{A\to +\infty }{\lim }[A(-{\lambda }_A)]$$

$$\frac{1}{2}\log \left(1+\frac{2}{x}\right)$$

$$\underset{x\to +\infty }{\lim }\left\{\frac{1}{x}L\left(\frac{1}{x}\right)\right\}$$

Full text

\WithSuffix

[7]G^ #1,#2_#3,#4(#5 #6| #7)

\WarningFilterrefcheckUnused label ‘sec:intro’

Asymptotic Exponentiality of the First Exit Time of the Shiryaev–Roberts Diffusion with Constant Positive Drift

**Aleksey S. Polunchenko

**Department of Mathematical Sciences, State University of New York at Binghamton,

Binghamton, New York, USA

000Address correspondence to A. S. Polunchenko, Department of Mathematical Sciences, State University of New York (SUNY) at Binghamton, 4400 Vestal Parkway East, Binghamton, NY 13902–6000, USA; Tel: +1 (607) 777–6906; Fax: +1 (607) 777–2450; E-mail: aleksey@binghamton.edu.

Abstract: We consider the first exit time of a Shiryaev–Roberts diffusion with constant positive drift from the interval $[0,A]$ where $A>0$. We show that the moment generating function (Laplace transform) of a suitably standardized version of the first exit time converges to that of the unit-mean exponential distribution as $A\to+\infty$. The proof is explicit in that the moment generating function of the first exit time is first expressed analytically and in a closed form, and then the desired limit as $A\to+\infty$ is evaluated directly. The result is of importance in the area of quickest change-point detection, and its discrete-time counterpart has been previously established—although in a different manner—by Pollak & Tartakovsky (2009).

Keywords: First exit times; Generalized Shiryaev–Roberts procedure; Laplace transform; Markov diffusion; Moment generating function; Quickest change-point detection; Whittaker functions.

Subject Classifications: 62L10; 60G40; 60J60.

.  ntroduction

This work centers around the so-called Generalized Shiryaev–Roberts (GSR) stochastic process, a time-homogeneous Markov diffusion well-known in the area of quickest change-point detection. See, e.g., Shiryaev, (1961, 1963, 1978, 2002, 2011), Pollak and Siegmund, (1985), Feinberg and Shiryaev, (2006), Burnaev et al., (2009), Polunchenko and Sokolov, (2016), and Polunchenko, (2016, 2017). More specifically, for a reason to be made clear shortly, the case of interest is that of the GSR process with a constant positive drift. Formally, we shall deal with the solution $(R_{t}^{r})_{t\geqslant 0}$ of the stochastic differential equation (SDE)

[TABLE]

where $\mu\neq 0$ is a given coefficient (whose meaning is explained below), and $(B_{t})_{t\geqslant 0}$ is standard Brownian motion (i.e., $\operatorname{\mathbb{E}}[dB_{t}]=0$, $\operatorname{\mathbb{E}}[(dB_{t})^{2}]=dt$, and $B_{0}=0$); the initial value $r$ is often referred to as the headstart. The process $(R_{t})_{t\geqslant 0}$ governed by (1.1) is a GSR process with unit drift and headstart $r\geqslant 0$; the unit drift can be trivially adjusted to any other constant positive level. The main contribution of this work concerns the distribution of the first exit time of $(R_{t}^{r})_{t\geqslant 0}$ from the interval $[0,A]$ with $A>0$ given, i.e., the stopping time:

[TABLE]

where $A>0$ is a preset level. Correspondingly, the process $(R_{t}^{r})_{t\geqslant 0}$ and its characteristics pose interest only up to the point of “extinction” at time instance $\mathcal{S}_{A}^{r}$, i.e., conditional on $\{\mathcal{S}_{A}^{r}>t\}$ for a given $t\geqslant 0$.

Just as does the GSR diffusion $(R_{t}^{r})_{t\geqslant 0}$—whether with constant positive drift or with a more general affine drift—the stopping time $\mathcal{S}_{A}^{r}$, too, plays a major role in the theory of quickest change-point detection: it is the Run Length of the so-called Generalized Shiryaev–Roberts (GSR) change-point detection procedure, set up to react to a possible shift in the drift of standard Brownian motion monitored “live”. Parameter $\mu$ present in the right-hand side of SDE (1.1) is the anticipated magnitude of the possible change in the drift. More concretely, equation (1.1) describes the dynamics of the GSR statistic $(R_{t}^{r})_{t\geqslant 0}$ in the pre-change regime, i.e., under the assumption that the drift $\mu$ has not yet “kicked in”, so that the observed Brownian motion is still “driftless”. Hence the stopping time $\mathcal{S}_{A}^{r}$ given by (1.2) is the GSR procedure’s Run Length to false alarm: at time instance $\mathcal{S}_{A}^{r}$ the GSR procedure sounds a false alarm, i.e., falsely declares the Brownian motion under surveillance as having gained a drift of size $\mu\neq 0$. The first moment of $\mathcal{S}_{A}^{r}$, i.e., $\operatorname{\mathbb{E}}[\mathcal{S}_{A}^{r}]$, is known in the change-point detection literature as the Average Run Length (ARL) to false alarm, and it is a popular metric of the “cost” of triggering a false alarm. Obviously $\operatorname{\mathbb{E}}[\mathcal{S}_{A}^{r}]$ increases with $A>0$, and, in particular, letting $A$ explode is the same as letting $\operatorname{\mathbb{E}}[\mathcal{S}_{A}^{r}]$ explode, and vice versa.

The GSR procedure was proposed by Moustakides et al., (2011) as a headstarted (hence, more general) version of the classical quasi-Bayesian Shiryaev–Roberts (SR) procedure that emerged from the independent work of Shiryaev (1961; 1963) and that of Roberts (1966). The interest in the GSR procedure (and its variations) is due to its recently discovered strong optimality properties. See, e.g., Burnaev, (2009), Feinberg and Shiryaev, (2006), Burnaev et al., (2009), Pollak and Tartakovsky, (2009), Polunchenko and Tartakovsky, (2010), Tartakovsky and Polunchenko, (2010), Vexler and Gurevich, (2011), and Tartakovsky et al., (2012).

We are now in a position to describe the specific contribution of this work. It is shown in the sequel that a suitably standardized version of the stopping time $\mathcal{S}_{A}^{r}$ given by (1.2) is asymptotically, as $A\to+\infty$, exponentially distributed with unit mean, for any headstart $R_{0}^{r}\triangleq r\geqslant 0$. Put another way, the GSR procedure’s Run Length to false alarm, properly scaled, is asymptotically, as the ARL to false alarm level explodes (i.e., as $\operatorname{\mathbb{E}}[\mathcal{S}_{A}^{r}]\to+\infty$), unit-mean exponential. More specifically, it is shown in the sequel that, as GSR procedure’s ARL to false alarm level gets large, the moment generating function (mgf) or the Laplace transform of a properly scaled version of the GSR stopping time $\mathcal{S}_{A}^{r}$ converges to that of the unit-mean exponential distribution. This implies convergence in distribution. The proof is explicit in that the mgf is first found analytically and in a closed form, and then the desired limit is shown directly to evaluate to the mgf of the unit-mean exponential distribution. It is also of note that the unit-drift assumption imposed on $(R_{t}^{r})_{t\geqslant 0}$ is essential, for it makes $(R_{t}^{r})_{t\geqslant 0}$ a (positive) recurrent process with all the ensuing consequences which ultimately “add up” to the desired asymptotic exponentiality of $\mathcal{S}_{A}^{r}$.

The discrete-time analogue of our result has been previously established—in an entirely different fashion—by Pollak and Tartakovsky (2009); see also Tartakovsky et al., (2008) and Yakir, (1998, 1995). As a matter of fact, Pollak and Tartakovsky (2009) proved the result not only for the GSR procedure, but for an entire class of Markov stopping times, which includes the GSR procedure as well as Page’s (1954) celebrated Cumulative Sum (CUSUM) “inspection” scheme. More importantly, Pollak and Tartakovsky (2009) also illustrated the importance of the result in the context of sequential change-point detection. Specifically, they argued that if the stopping time of a change-point detection procedure is asymptotically exponential under the no-change hypothesis, it is reasonable to expect it to be approximately exponentially distributed (under the no-change hypothesis) whenever the ARL to false alarm is large. Consequently, since the exponential distribution is fully characterized but its mean alone, the ARL to false alarm can be seen as indeed being an exhaustive metric of the false alarm risk. See, e.g., Tartakovsky, (2008) for a more detailed discussion of this issue. Moreover, Pollak and Tartakovsky (2009) also argued that the asymptotic exponentiality (in the pre-change regime) can be used for the evaluation of the change-point detection procedure’s local false alarm probabilities. As pointed out by Tartakovsky, (2005) these probabilities are of importance in a variety of applications. All these considerations obviously apply to the continuous-time setting considered in this work as well.

The rest of the paper is three sections. The first one, Section 2, is the paper’s main section, for this is where we formally state and then prove our main result. The second one, Section 3, is where we offer a short numerical study to complement and confirm our theoretical contribution experimentally. The third one, Section 4, is where we make a few concluding remarks and draw a line under the entire paper.

.  he Main Result

We first formally introduce the main object of study of this work. Let

[TABLE]

denote the mgf (Laplace transform) of the stopping time $\mathcal{S}_{A}^{r}$ given by (1.2). We are interested in the asymptotic behavior of $M(\alpha;A,x)$ as $A\to+\infty$. To that end, an important fact about $M(\alpha;A,x)$ is that, for any $\alpha\geqslant 0$, $A>0$ and $x\in[0,+\infty)$, it can actually be expressed analytically and in closed form through the spectral characteristics of the second-order differential operator

[TABLE]

i.e., the infinitesimal generator of the GSR diffusion $(R_{t}^{r})_{t\geqslant 0}$ governed by the SDE (1.1). More concretely, the operator $\mathscr{D}$ is restricted to the state space of $(R_{t}^{r})_{t\geqslant 0}$, i.e., the interval $[0,A]$, $A>0$, and the relevant spectral characteristics of $\mathscr{D}$ are the solutions $\lambda$ and $u(x,\lambda)$ of the Sturm–Liouville problem \big{[}\mathscr{D}\circ u\big{]}(x,\lambda)=\lambda\,u(x,\lambda) or explicitly

[TABLE]

subject to the boundary conditions

[TABLE]

which, translated into classical Feller’s (1952) boundary classification, cast $x=0$ as an entrance boundary for $(R_{t}^{r})_{t\geqslant 0}$, and $x=A$ as an absorbing boundary for $(R_{t}^{r})_{t\geqslant 0}$, i.e., the process is “killed” once it hits the right end of the interval $[0,A]$; in “differential equations speak”, the former condition is a Neumann-type boundary condition, while the latter condition is a Dirichlet-type boundary condition. It is apparent that the spectrum $\{\lambda\}$ of the operator $\mathscr{D}$ is dependent on $A>0$, and from now on, wherever necessary, we shall emphasize this dependence via the notation $\{\lambda_{A}\}$. Equation (2.3) subject to the boundary conditions (2.4) is a Sturm–Liouville problem, and it has recently received a renewed burst of attention in the literature on mathematical finance and quickest change-point detection. See, e.g., Linetsky, (2004, 2007), Collet et al., (2013), and notably Polunchenko, (2016, 2017). The work of Polunchenko (2017) will be referenced repeatedly throughout the sequel, following, for convenience, Polunchenko’s (2017) original notation.

We now turn to the work of Linetsky, (2007) and recall a general result from the interface between stochastic processes and Sturm–Liouville operator theory (theory of second-order self-adjoint differential operators); see also (Itô and McKean, , 1974, Chapter 4, Section 4.6) and (Borodin and Salminen, , 2002, Chapter II, Section 1.10). Let $(X_{t})_{t\geqslant 0}$ be a one-dimensional, time-homogeneous, regular Markov diffusion whose state space is some interval $(e_{1},e_{2})\subseteq\mathbb{R}$, where $-\infty\leqslant e_{1}<e_{2}\leqslant\infty$, and such that $X_{0}=x\in(e_{1},e_{2})$ is fixed. If $(X_{t})_{t\geqslant 0}$ is generated by the SDE $dX_{t}=a(X_{t})dt+\sqrt{b(X_{t})}\,dB_{t}$ where the diffusion coefficient $b(x)$ is continuous and strictly positive inside $(e_{1},e_{2})$ and the drift coefficient $a(x)$ is continuous on $(e_{1},e_{2})$, then the Laplace transform of the nonnegative random variable $\mathcal{T}_{y}^{x}\triangleq\inf\{t\geqslant 0\colon X_{t}=y\}$ with $\inf\{\varnothing\}=+\infty$ is given by

[TABLE]

where $\alpha>0$, and $\varphi(x;\alpha)$ and $\psi(x;\alpha)$ are two fundamental solutions $v(x;\alpha)$ of the equation

[TABLE]

subject to appropriate boundary conditions. Specifically, these fundamental solutions can be made unique (up to a multiplicative constant factor dependent on $\alpha$ but independent of $x$) by requiring that $\psi(x;\alpha)$ be an increasing function of $x$ subject to a boundary condition at $e_{1}$, while $\varphi(x;\alpha)$ be a decreasing function of $x$, subject a boundary condition at $e_{2}$.

To translate the above to our specific problem (2.3)–(2.4) it suffices to note that equation (2.3) can be easily converted to an equation of the form (2.6) by means of the integrating factor method. As a matter of fact, for our operator $\mathscr{D}$ given by (2.2), it has already been established, e.g., by Polunchenko and Sokolov, (2016) and by Polunchenko, (2016), that

[TABLE]

where

[TABLE]

and where $M_{a,b}(z)$ and $W_{a,b}(z)$ denote the so-called Whittaker $M$ and $W$ functions, respectively. The Whittaker functions are defined as the two fundamental solutions of the classical Whittaker (1904) equation

[TABLE]

where $w(z)$ is the unknown function of $z\in\mathbb{C}$, and $a,b\in\mathbb{C}$ are specified parameters; see, e.g., (Buchholz, , 1969, Chapter I). The Whittaker functions are typically considered in the cut plane $|\arg(z)\,|<\pi$ to ensure they are not multi-valued. For an extensive study of these functions and various properties thereof, see, e.g., Slater, (1960) and Buchholz, (1969).

At this point, in view of (2.5) and (2.7), we can conclude at once that

[TABLE]

where $\xi\equiv\xi(\lambda)$ is as in (2.8). Parenthetically, we remark that, apparently, this result, though relatively simple to obtain, was not previously known to the change-point detection community. It is also of note that the Laplace transform (2.9) can be inverted to yield the (pre-change) distribution of the GSR stopping time $\mathcal{S}_{A}^{r}$, and the inversion has already been performed by Polunchenko, (2017).

To proceed, observe that $\mathcal{S}_{A}^{r}$, by definition (1.2), almost surely explodes as $A\to+\infty$. Hence, it shouldn’t come as a surprise that, for any fixed $\alpha\geqslant 0$ and $x\in[0,+\infty)$, the limit of $M(\alpha;A,x)$ as $A\to+\infty$ is zero. Heuristically, this can be seen directly from the definition (2.1). More formally, one can appeal to the small-argument asymptotic behavior of the Whittaker $W$ function

[TABLE]

where here and onward $\Gamma(z)$ denotes the Gamma function (see, e.g., Abramowitz and Stegun, 1964, Chapter 6), to first get

[TABLE]

so that

[TABLE]

because $\xi(\alpha)\geqslant 1$ for $\alpha\geqslant 0$, and then conclude from the formula (2.9) for the mgf $M(\alpha;A,x)$ that the latter does, in fact, go to zero as $A\to+\infty$.

However, as previously noted by Pollak and Tartakovsky, (2009), there is a way to rescale $\mathcal{S}_{A}^{r}$ so as to get it to converge to a meaningful random variable as $A\to+\infty$; see also Tartakovsky et al., (2008). We now explain the idea.

Let

[TABLE]

denote the GSR statistic’s so-called quasi-stationary cumulative distribution function (cdf) and density, respectively. This time-invariant probability measure is independent of the GSR statistic’s headstart $R_{0}^{r}\triangleq r\in[0,A]$, and its existence can be inferred, e.g., from the work of Cattiaux et al., (2009); see also (Collet et al., , 2013, Section 7.8.2). Exact closed-form formulae for both $Q_{A}(x)$ and $q_{A}(x)$ have been recently obtained by Polunchenko, (2017). If the GSR statistic $(R_{t}^{r})_{t\geqslant 0}$ is started off a random point sampled from its quasi-stationary distribution, i.e., if $R_{0}^{r}\triangleq r\propto Q_{A}(x)$, then the statistical characteristics of the GSR statistic will be time-invariant, until the statistic hits the threshold $A$. Since the probability of hitting $A$ will be time-invariant as well, the distribution of the GSR stopping time will be exponential. More formally, define $(R_{t}^{Q})_{t\geqslant 0}$ as the solution of the SDE $dR_{t}^{Q}=dt+\mu R_{t}^{Q}dB_{t}$ with $R_{0}^{Q}\propto Q_{A}(x)$, and let $\mathcal{S}_{A}^{Q}\triangleq\inf\{t\geqslant 0\colon R_{t}^{Q}=A\}$ with $\inf\{\varnothing\}=+\infty$ and $A>0$. The stopping time $\mathcal{S}_{A}^{Q}$ is known in the quickest change-point detection literature as the randomized Shiryaev–Roberts–Pollak detection procedure, and it was originally proposed (for the discrete-time version of the problem) and first investigated by Pollak, (1985); it was also recently studied by Burnaev et al., (2009). Specifically, since $(0\geqslant)\,\lambda_{A}\triangleq\log\mathbb{P}(R_{t}^{Q}\geqslant A|\mathcal{S}_{A}^{Q}>t)$ is level for all $t\geqslant 0$, one can conclude that $\mathcal{S}_{A}^{Q}$ is exponentially distributed with parameter $-\lambda_{A}\,(>0)$, so that the product $-\lambda_{A}\mathcal{S}_{A}^{Q}$ is unit-mean exponential. As noted by Pollak and Tartakovsky, (2009), intuitively, the large-$A$ behavior of $\mathcal{S}_{A}^{r}$ for each fixed headstart is similar to that of $\mathcal{S}_{A}^{Q}$. Hence, it stands to reason that $-\lambda_{A}\mathcal{S}_{A}^{r}$ is approximately unit-mean exponential, whenever $A$ is large. Put another way, the right scaling factor for $\mathcal{S}_{A}^{r}$ is $-\lambda_{A}\triangleq-\log\mathbb{P}(R_{t}^{Q}\geqslant A|\mathcal{S}_{A}^{Q}>t)$.

The constant $\lambda_{A}\triangleq\log\mathbb{P}(R_{t}^{Q}\geqslant A|\mathcal{S}_{A}^{Q}>t)$ is the largest (nonpositive) eigenvalue of the operator $\mathscr{D}$ given by (2.2). For this kind of an operator it is known from the general Sturm–Liouville theory (see, e.g., Fulton et al., 1999) that its spectrum $\{\lambda\}$ is purely discrete, simple, located to the left of the origin (i.e., nonpositive), and is determined entirely by the Dirichlet condition (2.4), i.e., from the equation $u(A,\lambda)=0$ with $A>0$ fixed. More concretely, from (2.4) and (2.7) it can be readily seen that $\lambda_{A}$ is the largest (nonpositive) solution of the equation

[TABLE]

where $A>0$ is fixed and $\xi(\lambda)$ is as in (2.8). This equation was previously analyzed by Polunchenko, (2016, 2017), who, in particular, obtained order-one, order-two, and order-three asymptotic “large-$A$” approximations to $\lambda_{A}$. As an aside, we note that due to the discrete and simple nature of the spectrum of the operator $\mathscr{D}$ the range of values of $\alpha$ in the above formula (2.9) for the mgf $M(\alpha;A,x)$ can be extended from $\alpha\in[0,+\infty)$ to $\alpha\in(\lambda_{A},+\infty)$ where $\lambda_{A}\leqslant 0$ is largest (nonpositive) eigenvalue of $\mathscr{D}$.

The main contribution of this work can now be succinctly put as follows.

Theorem 2.1**.**

$\lim_{A\to+\infty}M(-\alpha\lambda_{A};A,x)=1/(1+\alpha)$* for any fixed $\alpha\in(-1,+\infty)$ and $x\in[0,+\infty)$; recall that $\lambda_{A}$ here is the largest (nonpositive) solution of equation (2.11).*

The plan for the remainder of this section is to prove this theorem. To that end, in view of (2.3), the problem essentially is to show that

[TABLE]

where $\xi=\xi(\lambda)$ is as in (2.8) and $\lambda_{A}$ is the largest (nonpositive) solution of equation (2.11).

The above limit can be evaluated by treating the numerator and the denominator separately. The key observation for either part is that $\lambda_{A}\nearrow 0$ as $A\to+\infty$, i.e., $\lambda_{A}$ is a monotonically increasing function of $A>0$, converging to 0 from below as $A\to+\infty$; see Polunchenko, (2017) for a proof. An immediate implication of this circumstance is that since $\lim_{A\to+\infty}\lambda_{A}=0$, then from (2.8) we also have $\lim_{A\to+\infty}\xi(\lambda_{A})=1$, and therefore

[TABLE]

because $W_{1,\tfrac{1}{2}}(z)=z\,e^{-\tfrac{z}{2}}$ which is a special case of (Buchholz, , 1969, Identity (28a), p. 23) asserting that $W_{a,a-\tfrac{1}{2}}(z)=z^{a}e^{-\tfrac{z}{2}}$. Hence, the numerator of the fraction under the limit (2.12) goes to unity as $A\to+\infty$.

It remains to take care of the denominator of the fraction under the limit (2.12), i.e., to show that

[TABLE]

which is a more delicate problem. Specifically, the problem is that not only the argument of the Whittaker $W$ function is dependent on $A$, but also its second index $\xi(\lambda_{A})/2$ which goes to $1/2$ as $A\to+\infty$ because $\lim_{A\to+\infty}\lambda_{A}=0$. As a result, the above small-argument asymptotics (2.10) of the Whittaker $W$ function is not “fine” enough and needs to be improved.

To that end, let us again turn to the work of Polunchenko, (2017) where the function $f(\lambda)\triangleq W_{1,\tfrac{1}{2}\xi(\lambda)}(z)$ was expanded into a Taylor series with respect to $\lambda$ around zero up to the third order for any fixed $z\geqslant 0$; it is noteworthy that $W_{a,b}(z)$ is an entire function of $b\in\mathbb{C}$ for any fixed $a\in\mathbb{R}$ and $z>0$. The expansion involves the following two special functions:

• •

The exponential integral

[TABLE]

see, e.g., (Abramowitz and Stegun, , 1964, Chapter 5); and

• •

Meijer’s (1936) celebrated $G$-function defined as the Mellin-Barnes integral

[TABLE]

where $\imath$ denotes the imaginary unit, i.e., $\imath\triangleq\sqrt{-1}$, the integers $m$, $n$, $p$, and $q$ are such that $0\leqslant m\leqslant q$ and $0\leqslant n\leqslant p$, and the contour of integration $\mathcal{C}$ is closed in an appropriate way to ensure the convergence of the integral. It is also required that no difference $a_{j}-b_{k}$ be an integer. The $G$-function is a very general function, and includes, as special cases, not only all elementary functions, but a number of special functions as well. An extensive list of special cases of the Meijer $G$-function can be found, e.g., in the classical special functions handbook of Prudnikov et al., (1990), which also includes a summary of the function’s basic properties. We will need the following particular case of the Meijer $G$-function:

[TABLE]

where $\operatorname{E_{1}}1(z)$ is the exponential integral defined in (2.15); see (Polunchenko, , 2017, Appendix A) for a proof.

We are now in a position present the third-order Taylor expansion obtained by Polunchenko, (2017) for the function $f(\lambda)\triangleq W_{1,\tfrac{1}{2}\xi(\lambda)}(z)$ with $z\geqslant 0$ fixed.

Theorem 2.2** (Polunchenko, 2017).**

For any $x\geqslant 0$ it holds true that

[TABLE]

where $\xi(\lambda)$ is as in (2.8), and

[TABLE]

with $\operatorname{E_{1}}1(x)$ and $G^{\,*,3}_{1,2}\lparen\begin{smallmatrix}3\\ 0,1\end{smallmatrix}|\,0,0,0\rparen{x}$ given by (2.15) and by (2.16), respectively.

This theorem readily gives the expansion

[TABLE]

which, as we shall see shortly, is more “fine” than necessary to pass $A\to+\infty$ and prove (2.14). To do so, we first recall yet another result of Polunchenko, (2017), viz. the double inequality

[TABLE]

where $\mu\neq 0$ is the parameter of the SDE (1.1). Hence $\lambda_{A}=-1/A+\mathcal{O}(A^{-3/2})$ so that

[TABLE]

The only issue is that the functions $L(x)$ and $G^{\,*,3}_{1,2}\lparen\begin{smallmatrix}3\\ 0,1\end{smallmatrix}|\,0,0,0\rparen{x}$ given, respectively, by (2.17) and (2.16), both go to infinity as $x$ goes to zero. However, in view of (Abramowitz and Stegun, , 1964, Inequality 5.1.20, p. 229) which states that

[TABLE]

from (2.15), (2.17) and (2.16) it can be seen that

[TABLE]

i.e., the functions $L(1/x)$ and $G^{\,*,3}_{1,2}\lparen\begin{smallmatrix}3\\ 0,1\end{smallmatrix}|\,0,0,0\rparen{1/x}$ both go to infinity as $x\to+\infty$ slower than $1/x$ goes to zero as $x\to+\infty$.

At this point Theorem 2.1, which is our main result, is straightforward to prove: it is merely a matter of using (2.19) and (2.20) in (2.18) to obtain (2.14), and then combining it with (2.13) to get (2.12), which in view of (2.9) is precisely the desired result.

To draw a line under this section, we remark that the formula (2.9) for the mgf $M(\alpha;A,x)$ of the GSR stopping time $\mathcal{S}_{A}^{r}$ can also be put to a more classical use, viz. to compute $\operatorname{\mathbb{E}}[(\mathcal{S}_{A}^{r})^{n}]$ for $n\geqslant 1$, i.e., to determine the actual moments of the GSR stopping time (under the no-change hypothesis). Specifically, since from definition (2.1) it is evident that

[TABLE]

and because formula (2.9) expresses $M(\alpha;A,x)$ explicitly as a quotient of two Whittaker $W$ functions, getting the $n$-th moment of the GSR stopping time essentially comes down to finding the derivatives, up through the $n$-th order inclusive, of the Whittaker $W$ function $W_{a,b}(z)$ with respect to the second index $b$. More concretely, from (2.9) it is direct to see that the required derivatives are of the following form:

[TABLE]

where $x\in[0,A]$, $A>0$, and $\xi(\lambda)$ is as in (2.8). Since the first three (for $n=1$, $2$ and $3$) of these derivatives are essentially given by Theorem 2.2 due to Polunchenko, (2017), computing $\operatorname{\mathbb{E}}[\mathcal{S}_{A}^{r}]$, $\operatorname{\mathbb{E}}[(\mathcal{S}_{A}^{r})^{2}]$, and $\operatorname{\mathbb{E}}[(\mathcal{S}_{A}^{r})^{3}]$, i.e., the first three moments of the GSR stopping time, is a matter of elementary algebra. The answer is:

[TABLE]

and

[TABLE]

where $r\in[0,A]$, $A>0$, and the functions $L(x)$ and $G^{\,*,3}_{1,2}\lparen\begin{smallmatrix}3\\ 0,1\end{smallmatrix}|\,0,0,0\rparen{x}$ are given, respectively, by (2.17) and (2.16). The first moment formula $\operatorname{\mathbb{E}}[\mathcal{S}_{A}^{r}]=A-r$ is well-known in quickest change-point detection, and was obtained—in an entirely different fashion—by Shiryaev, (1961, 1963) and many others. However, the second and third moment formulae appear to be new results. The fourth and higher moments can be found in a similar fashion, but since the formulae are far more cumbersome, they will be presented elsewhere.

.  umerical Results

To get a better sense as to how fast, as $A\to+\infty$, the random variable $-\lambda_{A}\mathcal{S}_{A}^{r}$ becomes unit-mean exponentially-distributed, we now offer a short numerical study where we assess the proximity of $\mathbb{P}(-\lambda_{A}\mathcal{S}_{A}^{r}\geqslant t)$ to $e^{-t}$ for various values of $A>0$, $r\in[0,A]$, $\mu\neq 0$, and $t\geqslant 0$. Specifically, since from Theorem 2.1 we can deduce that $\lim_{A\to+\infty}\mathbb{P}(-\lambda_{A}\mathcal{S}_{A}^{r}\geqslant t)=e^{-t}$ for any $t\geqslant 0$, or equivalently that $\lim_{A\to+\infty}\log\mathbb{P}(-\lambda_{A}\mathcal{S}_{A}^{r}\geqslant t)=-t$ for any $t\geqslant 0$, it is reasonable to expect the function $f(t)\triangleq\log\mathbb{P}(-\lambda_{A}\mathcal{S}_{A}^{r}\geqslant t)$ to be close to $-t$ for any $t\geqslant 0$, provided, however, that $A>0$ is sufficiently large. It is the proximity of $f(t)\triangleq\log\mathbb{P}(-\lambda_{A}\mathcal{S}_{A}^{r}\geqslant t)$ to the line $-t$ across a range of values of $t\geqslant 0$ that we shall use to judge how close the distribution of $-\lambda_{A}\mathcal{S}_{A}^{r}$ is to unit-mean exponential. The evaluation of $f(t)\triangleq\log\mathbb{P}(-\lambda_{A}\mathcal{S}_{A}^{r}\geqslant t)$ as a function of $t$ for any $A>0$, $r\in[0,A]$, and $\mu\neq 0$ is not a problem at all, because the survival function $\mathbb{P}(\mathcal{S}_{A}^{r}\geqslant t)$ of the GSR stopping time $\mathcal{S}_{A}^{r}$ given by (1.2) was recently found analytically and in a closed-form by Polunchenko, (2016) who also developed a Mathematica script to evaluate $\mathbb{P}(\mathcal{S}_{A}^{r}\geqslant t)$ and $\lambda_{A}$ each to within hundreds of decimal places of accuracy and for any $A>0$, $r\in[0,A]$, and $\mu\neq 0$.

To get started, let us first set $A=100$, which, in practice, would be considered low, so that the asymptotic exponentiality might not be quite in effect yet. Figures 1 show the obtained results for $t\in[0,10]$, $A=100$, $r=0$, and $\mu=\{1/2,1,3/2\}$. Specifically, Figure 1(a) shows $\log\mathbb{P}(-\lambda_{A}\mathcal{S}_{A}^{r}\geqslant t)$ as a function of $t\in[0,10]$, while Figure 1(b) shows the corresponding absolute error $\left|-\log\mathbb{P}(-\lambda_{A}\mathcal{S}_{A}^{r}\geqslant t)-t\right|$. For convenience, Figure 1(a) also includes the line $-t$ which $\log\mathbb{P}(-\lambda_{A}\mathcal{S}_{A}^{r}\geqslant t)$ is to converge to as $A\to+\infty$. An eye examination of these figures suggests that the distribution of $-\lambda_{A}\mathcal{S}_{A}^{r}$ nearly unit-mean exponential, even though $A$ is as low as $100$. The agreement with the asymptotic exponential distribution is even better when $A$ is larger.

Let us now keep $A$ at $100$ but increase the GSR statistic’s headstart $R_{0}^{r}\triangleq r$ to $r=50$. Since $A$ is only $100$, setting $r$ to half that is bringing $(R_{t}^{r})_{t\geqslant 0}$ much closer to $A$, thereby aiding the former to hit the latter sooner (on average). Put another way, increasing the headstart $r$ is, in some sense, akin to lowering the threshold $A>0$. As a result, the asymptotic exponentiality might not “kick in” as fast. This is exactly what we see in Figures 2(a) and 2(b) which show the obtained results for $t\in[0,10]$, $A=100$, $r=50$, and $\mu=\{1/2,1,3/2\}$. Again, Figure 2(a) also includes the line $-t$, but this time around the deviation of $\log\mathbb{P}(-\lambda_{A}\mathcal{S}_{A}^{r}\geqslant t)$ from $-t$ is noticeable with a naked eye. This is evidence that the asymptotic exponentiality isn’t quite there yet, and it is a direct consequence of the higher headstart value $r$.

However, if we keep $r$ at $50$ but increase $A$ to $500$, the distribution will get better aligned with the limiting exponential distribution, as can be seen from Figures 3 which show the results for $t\in[0,10]$, $A=500$, $r=50$, and $\mu=\{1/2,1,3/2\}$. Looking at Figures 3(a) and 3(b) we see that $-\lambda_{A}\mathcal{S}_{A}^{r}$ is fairly close to being a unit-mean exponential random variable.

.  oncluding Remarks

As was mentioned in the introduction, the obtained result, namely Theorem 2.1 which we proved explicitly, is the continuous-time equivalent of a similar result obtained earlier by Pollak and Tartakovsky (2009) in the discrete-time setting; see also Tartakovsky et al., (2008). On a practical level, we were able to confirm experimentally that the GSR stopping time is approximately exponential even if the detection threshold is fairly low. Pollak and Tartakovsky (2009) made the same observation in the discrete-time setting. Although we already elaborated in the introduction on the significance of our findings to theoretical change-point detection, it is also worth adding that some of the new special functions identities utilized in the paper may prove useful in other areas as well, e.g., in stochastic processes, stochastic differential equations, mathematical physics, and mathematical finance, where special functions arise quite often.

cknowledgement

The author’s effort was partially supported by the Simons Foundation via a Collaboration Grant in Mathematics under Award # 304574.