On the Convergence Rate of the Quasi- to Stationary Distribution for the Shiryaev-Roberts Diffusion

TL;DR

This paper analyzes the convergence rate of the quasi-stationary distribution of the Shiryaev-Roberts diffusion process, establishing an explicit bound that improves understanding of its asymptotic behavior as the boundary grows large.

Contribution

The paper provides the first explicit bound on the convergence rate of the quasi-stationary distribution for the Shiryaev-Roberts diffusion, using new bounds based on Bessel function properties.

Findings

Convergence rate is no worse than O(log(A)/A) as A→∞.

Constructed tight bounds for the quasi-stationary cdf using Bessel function properties.

Results are uniform in x≥0.

Abstract

For the classical Shiryaev--Roberts martingale diffusion considered on the interval $[0,A]$, where $A>0$ is a given absorbing boundary, it is shown that the rate of convergence of the diffusion's quasi-stationary cumulative distribution function (cdf), $Q_{A}(x)$, to its stationary cdf, $H(x)$, as $A\to+\infty$, is no worse than $O(\log(A)/A)$, uniformly in $x\ge0$. The result is established explicitly, by constructing new tight lower- and upper-bounds for $Q_{A}(x)$ using certain latest monotonicity properties of the modified Bessel $K$ function involved in the exact closed-form formula for $Q_{A}(x)$ recently obtained by Polunchenko (2017).

Equations

$$d{R}_t^r$$

$$R_t^r$$

$$H(x)$$

$$Q_A(x)$$

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

$$\mathrm{Ei}(x)$$

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

$${\int }_0^{+\infty }{e}^{-ay}\log (b+y)dy$$

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

$$\xi (\lambda )$$

$$q_A(x)$$

$$Q_A(x)$$

$$\frac{1}{A}+\frac{1}{A(1+A)}$$

$$W_{1,0}\left(\frac{2}{\tilde{A}}\right)$$

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

$$Q_{\infty }(x)$$

$$W_{a,a-\frac{1}{2}}(z)$$

$$W_{0,\frac{1}{2}}(z)$$

$$W_{0,b}(2z)$$

$$Q_A(x)$$

$$\frac{K_b(x_1)}{K_b(x_2)}$$

$$e^{\frac{2}{A}}H(x)$$

$$e^{x_2-x_1}{\left(\frac{x_2}{x_1}\right)}^{|b|}$$

$$Q_A(x)$$

$$f(b):=\log \left[\frac{K_b(x_1)}{K_b(x_2)}\right],\text{with}0<x_1<x_2,b\in \mathbb{R}$$

$$\frac{\partial }{\partial A}{Q}_A(x)$$

$$\frac{\partial }{\partial z}\left[e^{-\frac{z}{2}}{W}_{0,b}(z)\right]$$

$$\frac{\partial }{\partial z}\left[e^{-\frac{z}{2}}{z}^a{W}_{a,b}(z)\right]$$

$$\frac{\partial }{\partial A}{Q}_A(x)$$

Full text

\WithSuffix

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

\WarningFilterrefcheckUnused label ‘sec:intro’

On the Convergence Rate of the Quasi- to Stationary Distribution

for the Shiryaev–Roberts Diffusion

**Kexuan Li and 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: For the classical Shiryaev–Roberts martingale diffusion considered on the interval $[0,A]$, where $A>0$ is a given absorbing boundary, it is shown that the rate of convergence of the diffusion’s quasi-stationary cumulative distribution function (cdf), $Q_{A}(x)$, to its stationary cdf, $H(x)$, as $A\to+\infty$, is no worse than $O(\log(A)/A)$, uniformly in $x\geqslant 0$. The result is established explicitly, by constructing new tight lower- and upper-bounds for $Q_{A}(x)$ using certain latest monotonicity properties of the modified Bessel $K$ function involved in the exact closed-form formula for $Q_{A}(x)$ recently obtained by Polunchenko, 2017c .

Keywords: Generalized Shiryaev–Roberts procedure; Markov diffusion; Quickest change-point detection; Whittaker functions.

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

.  ntroduction

This work is an attempt to quantify the relationship between the phenomena of quasi-stationarity and stationarity exhibited by one particular version of the 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, 2017), Pollak and Siegmund, (1985), Feinberg and Shiryaev, (2006), Burnaev et al., (2008), Polunchenko and Sokolov, (2016), and Polunchenko, (2016); Polunchenko, 2017c ; Polunchenko, 2017a ; Polunchenko, 2017b . More specifically, the GSR process’ version of interest is the solution $(R_{t}^{r})_{t\geqslant 0}$ of the stochastic differential equation (SDE)

[TABLE]

where $(B_{t})_{t\geqslant 0}$ is standard Brownian motion in the sense that $\operatorname{\mathbb{E}}[dB_{t}]=0$, $\operatorname{\mathbb{E}}[(dB_{t})^{2}]=dt$, and $B_{0}=0$; the initial value $R_{0}^{r}\coloneqq r$ is often referred to as the process’ headstart. It is straightforward to solve (1.1) and express $(R_{t}^{r})_{t\geqslant 0}$ explicitly as

[TABLE]

so that the set $[0,+\infty)$ is easily seen to be the “natural” state space for $(R_{t}^{r})_{t\geqslant 0}$ because $R_{0}^{r}\coloneqq r\geqslant 0$. Moreover, it is also easily checked that $\operatorname{\mathbb{E}}[R_{t}^{r}-t-r]=0$ for any $t,r\geqslant 0$, i.e., the process $\{R_{t}^{r}-t-r\}_{t\geqslant 0}$ is a zero-mean martingale. Yet, although $(R_{t}^{r})_{t\geqslant 0}$ has a linear upward trend in time, it is actually a recurrent process with a nontrivial probabilistic behavior in the limit, as $t\to+\infty$; cf. (Pollak and Siegmund, , 1985, p. 270). Specifically, if $(R_{t}^{r})_{t\geqslant 0}$ is let run “loose”, i.e., considered on the entire nonnegative half-line, then the limiting (as $t\to+\infty$) behavior of $(R_{t}^{r})_{t\geqslant 0}$ is known as stationarity. The latter is characterized by the invariant probability measure whose cumulative distribution function (cdf) and density (pdf), respectively, are

[TABLE]

provided $r\in[0,+\infty)$. This probability measure has already been found, e.g., by Shiryaev, (1961, 1963), by Pollak and Siegmund, (1985), and more recently also by Feinberg and Shiryaev, (2006); Burnaev et al., (2008); Polunchenko and Sokolov, (2016), to be the momentless (no moments of orders one and higher) distribution

[TABLE]

which is an extreme-value Fréchet-type distribution, and a particular case of the inverse (reciprocal) gamma distribution. See also, e.g., Linetsky, (2004) and Avram et al., (2013). As an aside, note that, in view of (1.4), the stationary distribution of the reciprocal of $(R_{t}^{r})_{t\geqslant 0}$ is exponential with mean $1/2$.

However, if all states from a fixed $A>0$ and up inside the process’ “natural” state space $[0,+\infty)$ are made into absorbing states, then $(R_{t}^{r})_{t\geqslant 0}$ also has a nontrivial probabilistic behavior in the limit as $t\to+\infty$. This behavior is known as quasi-stationarity, and it is characterized by the invariant probability measure whose cdf and pdf, respectively, are

[TABLE]

provided $r\in[0,A)$. The existence of this probability measure was formally established, e.g., by Pollak and Siegmund, (1985), although one can also infer the same result, e.g., from the earlier seminal work of Mandl, (1961). Moreover, analytic closed-form formulae for $Q_{A}(x)$ and $q_{A}(x)$ were recently obtained by Polunchenko, 2017c , apparently for the first time in the literature; see formulae (3.3) and (3.4) in Section 3 below. These formulae were used by Polunchenko and Pepelyshev, (2018) to compute analytically the quasi-stationary distribution’s Laplace transform, and then also by Li et al., (2019) to find the quasi-stationary distribution’s fractional moment of any real order.

Remark 1.1**.**

The phenomenon of quasi-stationarity is also exhibited by $(R_{t}^{r})_{t\geqslant 0}$ in another case, viz. when all states from 0 up through a fixed $A>0$ inclusive inside the process’ “natural” state space $[0,+\infty)$ are made into absorbing states, so that the state space of $(R_{t}^{r})_{t\geqslant 0}$ becomes the set $[A,+\infty)$ with absorbtion at the lower end. This case was first investigated in (Collet et al., , 2013, Section 7.8.2). It was also recently analyzed by Polunchenko et al., (2018) who obtained analytically exact closed-form formulae for the quasi-stationary cdf and pdf.

Not surprisingly, the quasi-stationary distribution (1.5) and the stationary distribution (1.3) are related: as one would expect, the former converges to the latter as $A\to+\infty$. This was formally shown by Pollak and Siegmund, (1986), and not only for the GSR process at hand, but for an entire class of stochastically monotone processes. More specifically, it can be deduced from Pollak and Siegmund, (1986) that $Q_{A}(x)\geqslant H(x)$ for any fixed $A>0$ and any $x\geqslant 0$, and that $\lim_{A\to+\infty}Q_{A}(x)=H(x)$ for any fixed $x\geqslant 0$. The principal question addressed in this work is that of the actual rate of convergence of $Q_{A}(x)$ down to $H(x)$ as $A\to+\infty$ uniformly in $x\geqslant 0$: the answer—obtained and reported in Section 3—is $(0<)\;\sup_{x\geqslant 0}[Q_{A}(x)-H(x)]=O(\log(A)/A)$. This is the main contribution of this work, but not its only contribution: our $Q_{A}(x)$-to-$H(x)$ convergence (as $A\to+\infty$) analysis relies on new lower- and upper-bounds for $Q_{A}(x)$; the upperbounds are three, of varying complexity and accuracy, and all are much tighter than the trivial $Q_{A}(x)\leqslant 1$, while the lowerbound is one and it is much tighter than $H(x)\leqslant Q_{A}(x)$ implied by the earlier work of Pollak and Siegmund, (1986). All of the bounds, which may be considered another contribution of this work, are obtained explicitly with the aid of the formula for $Q_{A}(x)$ latterly offered by Polunchenko, 2017c , and certain recently discovered monotonicity properties of the modified Bessel $K$ function (of the second kind); the latter is involved in the formula for $Q_{A}(x)$ obtained by Polunchenko, 2017c .

The process $(R_{t}^{r})_{t\geqslant 0}$ governed by equation (1.1) arises in quickest change-point detection when the aim is to control the mean of the process $X_{t}\coloneqq t{\operatorname{\mathchoice{\rm 1\mskip-4.2mul}{\rm 1\mskip-4.2mul}{\rm 1\mskip-4.6mul}{\rm 1\mskip-5.2mul}}_{\left\{t>\nu\right\}}}+B_{t}$, where $\nu,t\geqslant 0$, observed “live”. Since $\operatorname{\mathbb{E}}[X_{t}]=t{\operatorname{\mathchoice{\rm 1\mskip-4.2mul}{\rm 1\mskip-4.2mul}{\rm 1\mskip-4.6mul}{\rm 1\mskip-5.2mul}}_{\left\{t>\nu\right\}}}$, it is anticipated that the drift of $(X_{t})_{t\geqslant 0}$ will change from none (zero) to one (per time unit) at time instance $\nu\in[0,+\infty]$ referred to as the change-point. The challenge is that $\nu$ is not known in advance; in particular $\nu=\infty$ is a possibility, i.e., the drift of $(X_{t})_{t\geqslant 0}$ may remain zero indefinitely and never change. The mean of $(X_{t})_{t\geqslant 0}$ is controlled by sounding an alarm should and as soon as the behavior of $(X_{t})_{t\geqslant 0}$ suggest that possibly $\operatorname{\mathbb{E}}[X_{t}]=t\neq 0$, i.e., $t>\nu$; if it is not the case, then the alarm is a false one. More concretely, the so-called GSR quickest change-point detection procedure, set up to control the drift of $(X_{t})_{t\geqslant 0}$, sounds a false alarm at

[TABLE]

where the constant $A>0$ is selected in advance in accordance with the desired false alarm risk level; it is to be understood in the right-hand side of the definition of $\mathcal{S}_{A}^{r}$ that $\inf\{\varnothing\}=+\infty$. Hence $(R_{t}^{r})_{t\geqslant 0}$ is the GSR procedure’s detection statistic in the pre-change regime, i.e., for $t\in[0,\nu]$. It is of note that $\mathbb{P}(\mathcal{S}_{A}^{r}<+\infty)=1$. The definition (1.5) of the quasi-stationary cdf can be rewritten as $Q_{A}(x)=\lim_{t\to+\infty}\mathbb{P}(R_{t}^{r}\leqslant x|\mathcal{S}_{A}^{r}>t)$.

The aforementioned GSR procedure was proposed by Moustakides et al., (2011) as a headstarted (i.e., 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 strong (near-) optimality properties. See, e.g., Burnaev, (2009), Feinberg and Shiryaev, (2006), Burnaev et al., (2008), Polunchenko and Tartakovsky, (2010), Tartakovsky and Polunchenko, (2010), Vexler and Gurevich, (2011), and Tartakovsky et al., (2012). For example, it is known that if the GSR procedure’s headstart is sampled from the quasi-stationary distribution (1.5), then such a randomization of the GSR procedure makes the latter nearly (to within a vanishingly small additive term) minimax in the sense of Pollak, (1985). The idea of such a randomization of the GSR procedure and a proof that the randomized GSR procedure is nearly minimax are due to Pollak, (1985) who was concerned with the discrete-time formulation of the problem. For the problem’s continuous-time formulation, the same result was established by Polunchenko, 2017b who heavily relied on the exact closed-form formulae for $Q_{A}(x)$ and $q_{A}(x)$ obtained by Polunchenko, (2016), as well as on the quasi-stationary distribution’s first two moments, also computed by Polunchenko, (2016).

The rest of the paper is three sections. The first one, Section 2, is to introduce our notation and to provide the necessary preliminary background on the special functions needed for our convergence analysis. The next section, Section 3, is the paper’s main section: this is where we formally state and prove our main result. Lastly, in Section 4 we make a few concluding remarks and wrap up the entire paper.

.  otation and nomenclature

For ease of exposition, we shall follow the standard notation employed uniformly across mathematical literature. This includes not only the usual symbols $\mathbb{R}$, $\mathbb{C}$, $\mathbb{N}$, $\mathbb{Z}$, $\mathrm{i}\coloneqq\sqrt{-1}$, and so on, but, more importantly, also an array of special functions we are to deal with throughout the sequel. These functions, in their most common notation, are:

The Gamma function $\Gamma(z)$, where $z\in\mathbb{C}$, sometimes also regarded as the extension of the factorial to complex numbers, due to the property $\Gamma(n)=(n-1)!$ exhibited for $n\in\mathbb{N}$. See, e.g., (Bateman and Erdélyi, 1953a, , Chapter 1). 2. 2.

The exponential integral function $\operatorname{Ei}(x)$, where $x\in\mathbb{R}\backslash\{0\}$, defined as

[TABLE]

with a singularity at $x=0$. Its basic properties are summarized, e.g., in (Abramowitz and Stegun, , 1964, Chapter 5). More specifically, we will need the function $\operatorname{E_{1}}1(x)\coloneqq-\operatorname{Ei}(-x)$ with $x>0$, so that

[TABLE]

where the second equality is because

[TABLE]

as given, e.g., by (Gradshteyn and Ryzhik, , 2007, Integral 4.337.1, p. 572). 3. 3.

The Whittaker $M$ and $W$ functions, traditionally denoted, respectively, as $M_{a,b}(z)$ and $W_{a,b}(z)$, where $a,b,z\in\mathbb{C}$; the Whittaker $M$ function is undefined when $-2b\in\mathbb{N}$, but can be regularized. These functions were introduced by Whittaker (1904) as the fundamental solutions to the Whittaker differential equation. See, e.g., Slater, (1960) and Buchholz, (1969). 4. 4.

The modified Bessel functions of the first and second kinds, conventionally denoted, respectively, as $I_{a}(z)$ and $K_{a}(z)$, where $a,z\in\mathbb{C}$; the index $a$ is referred to as the function’s order. See (Bateman and Erdélyi, 1953b, , Chapter 7). These functions form a set of fundamental solutions to the modified Bessel differential equation. The modified Bessel $K$ function is also known as the MacDonald function.

.  nalysis of the rate of convergence

To set the ground for our quasi- to stationary distribution convergence analysis we begin by recalling a few earlier results due to Polunchenko, 2017c . Specifically, it can be deduced from (Polunchenko, 2017c, , Theorem 3.1) that if $A>0$ is fixed and $\lambda\equiv\lambda_{A}>0$ is the smallest (positive) solution of the equation

[TABLE]

where

[TABLE]

then the quasi-stationary pdf is given by

[TABLE]

and the respective cdf is given by

[TABLE]

and $q_{A}(x)$ and $Q_{A}(x)$ are each a sufficiently smooth function of $x\geqslant 0$ and $A>0$; observe from (3.1), (3.2), and (3.3) that $q_{A}(A)=0$. The smoothness of $q_{A}(x)$ and $Q_{A}(x)$ is due to certain analytic properties of the Whittaker $W$ function on the right of formulae (3.3) and (3.4). These formulae stem from the solution of a certain Sturm–Liouville problem, and $\lambda$ is the smallest positive eigenvalue of the corresponding Sturm–Liouville operator; if the Sturm–Liouville operator is negated, as was done by Polunchenko, 2017c , then $\lambda$ becomes the operator’s largest negative eigenvalue.

Remark 3.1**.**

The definition (3.2) of $\xi(\lambda)$ can actually be changed to $\xi(\lambda)\coloneqq-\sqrt{1-8\lambda}$ with no effect whatsoever on either equation (3.1), or formulae (3.3) and (3.4), i.e., all three are invariant with respect to the sign of $\xi(\lambda)$. This was previously pointed out by Polunchenko, 2017c , and the reason for this $\xi(\lambda)$-symmetry is because equation (3.1) and formulae (3.3) and (3.4) each have $\xi(\lambda)$ present only as (double) the second index of the corresponding Whittaker $W$ function or functions involved, and the Whittaker $W$ function in general is known (see, e.g., (Buchholz, , 1969, Identity (19), p. 19)) to be an even function of its second index, i.e., $W_{a,b}(z)=W_{a,-b}(z)$.

It is evident that equation (3.1) is a key component of formulae (3.3) and (3.4), and consequently, of all of the characteristics of the quasi-stationary distribution as well. As a transcendental equation, it can only be solved numerically, although to within any desired accuracy, as was previously done, e.g., by Linetsky, (2004); Polunchenko, (2016); Polunchenko, 2017c ; Polunchenko, 2017a , with the aid of Mathematica developed by Wolfram Research: Mathematica’s special functions capabilities have long proven to be superb. Yet, it is known (see, e.g., Linetsky, 2004 and Polunchenko, 2016) that for any fixed $A>0$, the equation has countably many simple solutions $0<\lambda_{1}<\lambda_{2}<\lambda_{3}<\cdots$, such that $\lim_{k\to+\infty}\lambda_{k}=+\infty$. All of them depend on $A$, but since we are interested only in the smallest one, we shall use either the “short” notation $\lambda$, or the more explicit $\lambda_{A}$ to emphasize the dependence on $A$. It was shown by Polunchenko, (2016) that $\lambda_{A}$ is a monotonically decreasing function of $A$, such that

[TABLE]

whence $\lim_{A\to+\infty}\lambda_{A}=0$, and more specifically $\lambda_{A}=A^{-1}+O(A^{-3/2})$; cf. (Polunchenko, , 2016, p. 136 and Lemma 3.3). See also Polunchenko and Pepelyshev, (2018) for a discussion of potential ways to improve the foregoing double inequality.

Remark 3.2**.**

*Since $\lambda\equiv\lambda_{A}$ is monotonically decreasing in $A$, and such that $\lim_{A\to+\infty}\lambda_{A}=0$, one can conclude from (3.2) that $\xi(\lambda_{A})$, for any finite $A>0$, is either

(a) purely imaginary (i.e., $\xi(\lambda)=\mathrm{i}\alpha$ where $\mathrm{i}\coloneqq\sqrt{-1}$ and $\alpha\in\mathbb{R}$) if $A$ is sufficiently small, or (b) purely real and between 0 inclusive and 1 exclusive (i.e., $0\leqslant\xi(\lambda)<1$) otherwise . The borderline case is when $\xi(\lambda)=0$, i.e., when $\lambda_{A}=1/8$, and the corresponding critical value of $A$ is the solution $\tilde{A}>0$ of the equation*

[TABLE]

as can be established by a basic numerical calculation. Hence, if $A<\tilde{A}\approx 10.240465$, then $\lambda_{A}>1/8$ so that $\xi(\lambda)$ is purely imaginary; otherwise, if $A\geqslant\tilde{A}\approx 10.240465$, then $\lambda_{A}\in(0,1/8]$ so that $\xi(\lambda)$ is purely real and such that $\xi(\lambda)\in[0,1)$ with $\lim_{A\to+\infty}\xi(\lambda_{A})=1$.

Since the SR process is a stochastically monotone Markov process, it can be concluded at once from Pollak and Siegmund, (1986) that $Q_{A}(x)\downarrow H(x)$ as $A\to+\infty$ for every fixed $x\geqslant 0$. However, since $Q_{A}(x)$ is given explicitly by formula (3.4), it is of interest to see if the same conclusion can be reached in a more explicit fashion. To this end, recall first the limit

[TABLE]

which was previously proved by Polunchenko, 2017a . On account of this limit it can be seen from (3.4) that

[TABLE]

because $\xi(\lambda_{A})\to 1$, as $A\to+\infty$, by Remark 3.2. Now, since

[TABLE]

as given, e.g., by (Buchholz, , 1969, Identity (28a), p. 23), so that

[TABLE]

the (expected) conclusion that $Q_{\infty}(x)=H(x)$ for every $x\geqslant 0$ becomes apparent.

To show that the convergence of $Q_{A}(x)$ to $H(x)$ as $A\to+\infty$ is from above, it is convenient to first rid the formula (3.4) for $Q_{A}(x)$ of the two Whittaker $W$ functions involved in it, and instead express it through a ratio of two modified Bessel $K$ functions (of the second kind): the Bessel functions and ratios thereof appear to have been studied more extensively than the Whittaker $W$ function. Specifically, observe that since

[TABLE]

as given, e.g., by (Abramowitz and Stegun, , 1964, Identity 9.6.48, p. 377), formula (3.4) can be rewritten equivalently as

[TABLE]

and this formula, though completely equivalent to formula (3.4), will prove to be more convenient for our purposes than (3.4).

We are now ready to show that $Q_{A}(x)\geqslant H(x)$ for every $x\geqslant 0$. Specifically, let us focus only on the case of $x\in(0,A)$, because the result is trivial for $x\in\{0\}\cup[A,+\infty)$. The idea is to appeal to the inequality

[TABLE]

which was first stated as a conjecture by (Baricz, , 2010, pp. 589–590), but its first formal proof is apparently due to Yang and Zheng, (2017); cf. (Yang and Zheng, , 2017, Inequality (3.2), p. 2951). From inequality (3.10), formula (3.9), and Remark 3.2 whereby $\xi(\lambda_{A})\in[0,1]$ for any $A\geqslant\tilde{A}\approx 10.240465$, one can conclude that

[TABLE]

which is a much tighter lowerbound for $Q_{A}(x)$ than just $H(x)$ itself, and may thus be considered a new result in its own right.

Similarly, by (Baricz, , 2010, Inequality (3.3), p. 580), which states that

[TABLE]

and formula (3.9), we find that

[TABLE]

which, as will be demonstrated shortly, is a fairly stringent upperbound, especially if $A$ is large; this bound is a new result in itself, too.

It is evident that together inequalities (3.11) and (3.13) imply that $Q_{A}(x)\downarrow H(x)$ as $A\to+\infty$ for every fixed $x\geqslant 0$. However, it is actually possible to show more, thanks to (Yang and Zheng, , 2017, Theorem 2.6(iii), p. 2948), whereby the function

[TABLE]

is convex. As we shall now prove explicitly, this recently established property of the Bessel $K$ function enables one to conclude that

[TABLE]

for every $x\in(0,A)$, so that $Q_{A_{1}}(x)\geqslant Q_{A_{2}}(x)$ for every $x\geqslant 0$ whenever $A_{2}\geqslant A_{1}\geqslant\tilde{A}\approx 10.240465$; recall Remark 3.2 and equation (3.6) which defines $\tilde{A}\approx 10.240465$.

One way to find the partial derivative of $Q_{A}(x)$ with respect to $A>0$ with $x\in(0,A)$ assumed constant is via direct differentiation of formula (3.4). The necessary Whittaker $W$ function derivative identity is

[TABLE]

which is a special case of the more general identity

[TABLE]

given by, e.g., (Slater, , 1960, Formula (2.4.24), p. 25); incidentally, identity (3.15) is also a way to get formula (3.3) for the quasi-stationary distribution’s pdf $q_{A}(x)$ from formula (3.4) for $Q_{A}(x)$, and vise versa. By virtue of identities (3.15) and (3.8), and equation (3.1), the derivative of $Q_{A}(x)$ with respect to $A$ simplifies to

[TABLE]

and we note the similarity of the $\log$ of the ratio of two Bessel $K$ function to the setting of Yang and Zheng, (2017). Now observe that $\xi(\lambda_{A})$ defined by (3.2) is a monotonically increasing function $A>0$, for $\lambda_{A}>0$ is a monotonically decreasing function $A>0$. Hence, if we could show that

[TABLE]

then the desired conclusion (3.14) would follow. This is where (Yang and Zheng, , 2017, Theorem 2.6(iii), p. 2948) we mentioned earlier comes in. Specifically, since $\xi(\lambda_{A})$ is between 0 and 1 for any $A\geqslant\tilde{A}$, we have

[TABLE]

but

[TABLE]

and the next question is to evaluate each of the two derivatives: one at $b=0$ and one at $b=1/2$. The derivative at $b=0$ is zero, and the reasons it is zero are two. The first is the fact that $K_{b}(z)>0$ for all $z>0$ and $b\in\mathbb{R}$; cf., e.g., (Yang and Zheng, , 2017, p. 2944). The second reason is the identity

[TABLE]

as given, e.g., by (Brychkov, , 2008, Identity 1.14.2.1, p. 40). For the derivative at $b=1/2$ we find

[TABLE]

where $\operatorname{E_{1}}1(z)$ denotes the exponential integral function (2.2). This result is an immediate consequence of (Watson, , 1922, Identity (13), p. 80) which states that

[TABLE]

and (Brychkov, , 2008, Identity 1.14.2.3, p. 40) which states that

[TABLE]

where $\operatorname{Ei}(z)$ denotes the exponential integral function (2.1).

At this point we can conclude that

[TABLE]

for every $x\in(0,A)$. This shows (3.14), i.e., gives the desired conclusion that $Q_{A}(x)$ is a nonincreasing function of $A$ for every fixed $x\in(0,A)$, at least for $A\geqslant\tilde{A}\approx 10.240465$; recall, again, that $\tilde{A}\approx 10.240465$ was introduced in Remark 3.2.

Remark 3.3**.**

While the assumption made earlier that $A\geqslant\tilde{A}\approx 10.240465$ is required to ensure the index shared by the two Bessel $K$ functions involved in formula (3.9) is purely real, so that inequalities (3.10) and (3.12) do apply, and yield the bounds (3.11) and (3.13) for $Q_{A}(x)$, these bounds actually seem to remain valid for $A<\tilde{A}\approx 10.240465$ as well. Specifically, we carried out an extensive numerical experiment to test the double inequality

[TABLE]

and could not find a single value of $A>0$ for which the double inequality would not hold. For $A\geqslant\tilde{A}\approx 10.240465$, the right half of this inequality coincides with (3.13), and for $A<\tilde{A}\approx 10.240465$, it reduces down to

[TABLE]

which we determined through the numerical experiment to be somewhat conservative.

We implemented the obtained lower- and upper-bounds (3.11) and (3.13) in a Mathematica script, and used the script to produce Figures 1, 2, and 3. These figures show the performance of the bounds as functions of $x\in[0,A]$ for $A=10$, $50$, and $100$; for $A=10$, which is less than $\tilde{A}\approx 10.240465$, the upperbound was computed using (3.18). Specifically, Figures 1(a), 2(a), and 3(a) show the quasi-stationary distribution’s cdf $Q_{A}(x)$, the stationary cdf $H(x)$, and the lower- and upper-bounds (3.11) and (3.13)—all as functions of $x\in[0,A]$, for $A=10$, $50$, and $100$, respectively. To solve equation (3.1) and compute $\lambda_{A}$, and, subsequently, also evaluate $Q_{A}(x)$, we relied on the Mathematica script written earlier by Polunchenko, 2017c . The figures leave the impression that the bounds perform reasonably well, even though in quickest change-point detection any value of $A$ that is less than $100$ is generally considered low because the corresponding false alarm risk level in high. The only “exception” is the case of $A=10$: the bounds in this case are somewhat loose, especially the upperbound, whose corresponding relative deviation from $Q_{A}(x)$ peaks almost $100\,\%$. However, expectedly, as $A$ increases, the quality of the bounds improves, and the gap between $Q_{A}(x)$ and $H(x)$ diminishes as well, uniformly for all $x\in[0,A]$. For $A=50$ and higher, either bound’s largest margin of error is adequate. Juxtaposed side-by-side to Figures 1(a), 2(a), and 3(a) are Figures 1(b), 2(b), and 3(b) which show the corresponding lower- and upper-bound errors. Specifically, the lowerbound error is defined as the excess of $Q_{A}(x)$ over the bound, while the upperbound error is defined as the excess of the bound over $Q_{A}(x)$. All three error figures suggest that the upperbound is generally tighter than the lowerbound, unless $A$ is very low, in which case both bounds are off by a large margin.

We are now in position to explain the main contribution of this work. It concerns the actual speed at which $Q_{A}(x)$ descends down to $H(x)$, as $A\to+\infty$, for every $x\geqslant 0$, and can be formally and succinctly stated thus:

[TABLE]

and it is shown next.

Our proof of the foregoing asymptotics relies on the inequality

[TABLE]

valid for any $x\in[0,A]$ and at least for $A\geqslant\tilde{A}\approx 10.240465$; recall Remark 3.2 and that $\tilde{A}$ is the solution of equation (3.6). From this inequality and the lowerbound (3.11) it can be seen that $Q_{A}(x)$ descends down to $H(x)$ no faster than $O(1/A)$ but no slower than $O(\log(A)/A)$. The upper estimate for the convergence speed is because

[TABLE]

which is trivial. The lower estimate is because

[TABLE]

and to see this observe that $1-\xi(\lambda_{A})=O(1/A)$, as can be concluded from (3.5) plugged into (3.2), and recall that \lim_{x\to+\infty}\big{(}\sqrt[x]{x}\,\big{)}=1 and that

[TABLE]

as given by (Abramowitz and Stegun, , 1964, Inequality 5.1.20, p. 229);

To prove the above inequality for $Q_{A}(x)-H(x)$, fix $A\geqslant\tilde{A}\approx 10.240465$ and $x\in(0,A)$, and Taylor-expand $Q_{A}(x)$ given by (3.9) regarded as a function of $\xi(\lambda)/2\in[0,1/2)$ near $1/2$ up to the linear term; the assumption made here that $\xi(\lambda_{A})$ is nonnegative is not restrictive, because $Q_{A}(x)$ is an even function of $\xi(\lambda_{A})$, as can be inferred from formula (3.4) and Remark 3.1. On account of (3.12) and its consequence (3.16) this gives

[TABLE]

for some $b^{*}=b^{*}(x,A)\in(\xi(\lambda)/2,1/2)$. Now, since the second term in the obtained Taylor expansion is nonneative, the expansion can be easily turned into a new upperbound for $Q_{A}(x)$: it suffices to apply (3.13) along with (Yang and Zheng, , 2017, Theorem 2.6(iii), p. 2948) to upperbound the expression dependent on $b^{*}$ on the right of the Taylor expansion. Specifically, the new upperbound for $Q_{A}(x)$ takes the form

[TABLE]

and it is generally tighter than the upperbound (3.13) we obtained earlier.

To conclude, we remark that one can improve the upperbound for $Q_{A}(x)-H(x)$ by Taylor-expanding $Q_{A}(x)$ given by (3.9) regarded as a function of $\xi(\lambda)/2\in[0,1/2)$ near $1/2$ further, viz. up to the quadratic term. Specifically, a simple calculation using (3.13) and (Yang and Zheng, , 2017, Theorem 2.6(iii), p. 2948) gives

[TABLE]

where again $A\geqslant\tilde{A}\approx 10.240465$ with $\tilde{A}\approx 10.240465$ defined by equation (3.6) and Remark 3.2.

.  onclusion and outlook

The obtained new bounds for the cdf of the quasi-stationary distribution of the Generalized Shiryaev–Roberts process do have applications in quickest change-point detection. For one, the upperbounds can be used to tightly cap Pollak’s (1985) maximal Average Detection Delay (ADD) delivered by Pollak’s (1985) randomized Shiryaev–Roberts–Pollak (SRP) procedure, and subsequently show that the SRP procedure is nearly minimax in the sense of Pollak (1985); the rate of convergence to the unknown optimal minimax performance can be estimated as well. See, e.g., Polunchenko, 2017b who upperbounded the quasi-stationary cdf by unity, which is, of course, a very conservative upperbound. More importantly, the upperbounds and the lowerbound offered in this work may also be instrumental in showing that the Generalized Shiryaev–Roberts procedure with a carefully designed headstart is, too, nearly minimax in the sense of Pollak (1985). We are of the view that any proof of such a strong claim would certainly merit a separate paper. For example, Tartakovsky et al., (2012) confirm the claim to be valid in the discrete-time setting. In the continuous-time setting, the main challenge is the fact that while all of the relevant performance characteristics can be computed analytically and in a closed form, the expressions are quite involved, and depend on special functions; see, e.g., Polunchenko, (2016). However, the bounds obtained in this work may enable one to find an explicit, simpler yet sufficiently tight upperbound for Pollak’s (1985) maximal ADD delivered by the GSR procedure, and then use the upperbound to show that its excess over the unknown optimal minimax ADD is negligible in the limit, as the Average Run Length to false alarm goes to infinity. Some preparatory work in this direction is already underway.

cknowledgement

The effort of A.S. Polunchenko was partially supported by the Simons Foundation via a Collaboration Grant in Mathematics under Award # 304574.