3.1. Skovgaard double-saddlepoint approximations

The first method is simple and easy to describe as it uses the Skovgaard double-saddlepoint approximation [Reference Skovgaard18]. This method amounts to taking the double GF $\mathcal{R}(\mathrm{y},\mathrm{z})$ and treating it as a bivariate PGF in variables R(T) and T, where R(T) given $T=n$ has the distribution of R(n), and T has an improper uniform distribution over $\mathbb{N}$ specified by counting measure. This is exactly how $\mathcal{R}(\mathrm{y},\mathrm{z})$ is computed in (1) as the GF of $\{\mathcal{R}_{n}(\mathrm{y})\}$ in $n\leftrightsquigarrow\mathrm{z}$ , with variable n associated with counting measure. This is not a proper mass function since $\mathcal{R}(1, 1)=\infty$ , but the point made in [Reference Butler6, Section 2] is that this doesn’t matter: the Skovgaard approximation remains valid despite this. Thus we proceed by treating $\mathcal{M}(\mathrm{r},\mathrm{s})=\mathcal{R}({\mathrm{e}}^{\mathrm{r}},{\mathrm{e}}^{\mathrm{s}})$ as the joint MGF of the random pair $\{R(T),T\}$ , and approximate the survival function of R(n) as the Skovgaard double-saddlepoint approximation for the conditional survival function $\mathbb{P}\{R(T)\geq m\mid T=n\}$ . We treat $\mathcal{K}(\mathrm{r},\mathrm{s})=\ln\mathcal{R}({\mathrm{e}}^{\mathrm{r}},{\mathrm{e}}^{\mathrm{s}})$ as the joint cumulant-generating function (CGF) and from it determine two continuity-corrected double-saddlepoint approximations $\mathbb{P}_{D1}$ or $\mathbb{P}_{D2}$ as described in the Supplementary Material and further discussed in [Reference Borovkov1, Section 4.2]. Approximation $\mathbb{P}_{D1}$ uses the solution $(\hat{r},\hat{s})$ to the saddlepoint equations

(5) \begin{equation} m = \mathcal{K}^{\prime}_{\mathrm{r}}(\hat{r},\hat{s}), \qquad n=\mathcal{K}^{\prime}_{\mathrm{s}}(\hat{r},\hat{s}), \end{equation}

where $\mathcal{K}^{\prime}_{\mathrm{r}}=\partial\mathcal{K}(\mathrm{r},\mathrm{s})/\partial\mathrm{r}$ . Approximation $\mathbb{P}_{D2}$ uses the saddlepoint solution $(\tilde{r},\tilde{s})$ to $m-\frac12 = \mathcal{K}^{\prime}_{\mathrm{r}}(\tilde{r},\tilde{s})$ and $n=\mathcal{K}^{\prime}_{\mathrm{s}}(\tilde{r},\tilde{s})$ , with m offset to the value $m-\frac12$ .

Table 2. Exact values $(\mathbb{P})$ for $\mathbb{P}\{\mathcal{R}(40)\geq m\}$ are compared with two double-saddlepoint approximations ( $\mathbb{P}_{D1}$ and $\mathbb{P}_{D2}$ ), two residue-saddlepoint approximations ( $\mathbb{P}_{R1}$ and $\mathbb{P}_{R2})$ from Section 3.2, and a continuity-corrected Normal $(3.79,3.69)$ approximation (Norm) from Section 4. The bold digit indicates the last ‘accurate’ digit or the last digit in agreement with the exact result when both are rounded to the same number of digits. In the notation used here, $0.0^{2}43=0.0043$ , etc. A Binomial $(40,0.1)$ upper bound (Binomial) is shown.

Example 5. (Viral RNA tests.) This example provides the means for judging the accuracy of the double-saddlepoint approximations since it is possible to compute the exact survival function of R(40) up to $m=10$ . The exact value of $\mathcal{R}_{n}(\mathrm{y})$ is

(6) \begin{align} \mathcal{R}_{n}(\mathrm{y}) & = [\mathrm{z}^{n}]\frac{\mathcal{F}_{S}(\mathrm{z})}{1-\mathcal{H}(\mathrm{y,z})} \nonumber \\ & =v[\mathrm{z}^{n}]\frac{(1-p_{0\cdot}\mathrm{z})^{-1}} {1-\mathrm{z}B_{1}(\mathrm{y})/\{1-\mathrm{z}B_{0}(\mathrm{y})\}} = \frac{B_{1}(\mathrm{y})}{B(\mathrm{y})-p_{0\cdot}}B(\mathrm{y})^{n} + \frac{B_{0}(\mathrm{y})-p_{0\cdot}}{B(\mathrm{y})-p_{0\cdot}}p_{0\cdot}^{n}, \end{align}

where $B(\mathrm{y})=B_{0}(\mathrm{y})+B_{1}(\mathrm{y})=p_{\cdot1}\mathrm{y}+p_{\cdot0}$ . Expression (6) follows from a partial fraction expansion with exact inversion in z $\leftrightsquigarrow n$ . The true survival function can now be computed by expanding $\{1-\mathrm{y}\mathcal{R}_{40}(\mathrm{y})\}/(1-\mathrm{y})$ in a power series in y using Maple, since

\begin{align*} \mathbb{P}\{\mathcal{R}(40)\geq m\} = [\mathrm{y}^{m}]\frac{1-\mathrm{y}\mathcal{R}_{40}(\mathrm{y})}{1-\mathrm{y}}. \end{align*}

Such expansions in Maple failed for $n>40$ and this motivates our choice of $n=40$ . The exact values are shown in Table 2 as $\mathbb{P}$ , where they are compared with the two double-saddlepoint approximations $\mathbb{P}_{D1}$ and $\mathbb{P}_{D2}$ .

For this example, we might have mistakenly thought that R(n) has a Binomial $(n,p_{\cdot1})$ distribution so that $\mathcal{R}_{n}(\mathrm{y})$ should be $B(\mathrm{y})^{n}$ rather than the more complicated expression $\mathcal{R}_{n}(\mathrm{y})$ in (6). As such, then at time (test) n we should expect a Binomial $(n,p_{\cdot1})$ count of true positives. This reasoning is incorrect, and this is related to how R(n) packages the true positives. After n tests, there are N(n) packages each with a Geometric $(p_{1\cdot})$ number of tests so we are not counting rewards for exactly n tests. The counts from an unfulfilled renewal span at time n do not get counted so that if $C\sim\mathrm{Binomial}(n,p_{\cdot1})$ , then $R(n)\leq C$ with probability 1. Thus, the last table entry in Table 1 provides an upper bound to the true probabilities which becomes a sharper bound as n increases.

Figure 1 plots the percentage relative error for the approximations, given as

\begin{align*} \%\text{ rel. err } = 100\frac{\mathbb{P}_{D1}\mathcal{\{}R(40)\geq m\}-\mathbb{P}(m)} {\min\{\mathbb{P}(m),1-\mathbb{P}(m)\}}, \end{align*}

where $\mathbb{P}(m)=\mathbb{P}\{\mathcal{R}(40)\geq m\}$ . This percentage respects both the sign of the difference as well as the tails by dividing by the smaller value of $\mathbb{P}(m)$ and $1-\mathbb{P}(m)$ . Approximations $\mathbb{P}_{D1}$ (crosses) and $\mathbb{P}_{D2}$ (diamonds) are shown along with the residue-saddlepoint approximations $\mathbb{P}_{R1}$ (solid red line) and $\mathbb{P}_{R2}$ (dashed blue line) discussed in Section 3.2 and a continuity-corrected normal approximation (boxes) developed in Section 4.

Figure 2 shows the boundary for the convergence domain of the double MGF $\mathcal{R}({\mathrm{e}}^{\mathrm{r}},{\mathrm{e}}^{\mathrm{s}})$ as a solid line. Plotted within its convergence domain below as circles are the saddlepoints $(\hat{r},\hat{s})$ solving (5) for the $\mathbb{P}_{D1}$ approximation.

Figure 1. Percentage relative errors respecting the tails for $\mathbb{P}_{D1}$ (crosses), $\mathbb{P}_{D2}$ (diamonds), $\mathbb{P}_{R1}$ (solid red line), $\mathbb{P}_{R2}$ (dashed blue line), and a continuity-corrected Normal $(3.79,3.69)$ approximation (boxes). Missing boxes for the normal approximation are out of the range $-1\%$ –11%.

Figure 2. Plot of the convergence boundary (solid line) for the double MGF $\mathcal{F}_{S}({\mathrm{e}}^{\mathrm{s}})/\{1-\mathcal{H}({\mathrm{e}}^{\mathrm{r}},{\mathrm{e}}^{\mathrm{s}})\}$ in the viral RNA example. Circles within the convergence domain show saddlepoints for the $\mathbb{P}_{D1}$ approximation for $m=1$ (upper left) up to $m=10$ (lower right).

In general, the convergence domain of the double MGF $\mathcal{R}({\mathrm{e}}^{\mathrm{r}},{\mathrm{e}}^{\mathrm{s}})$ is $\mathcal{B}\cap\mathcal{C}$ where $\mathcal{B}=\{(\mathrm{r},\mathrm{s})\in\mathbb{R}^{2}\,:\,\mathcal{H}({\mathrm{e}}^{\mathrm{r}},{\mathrm{e}}^{\mathrm{s}})<1\}$ and $\mathcal{C}=\{(\mathrm{r},\mathrm{s})\in\mathbb{R}^{2}\,:\,\mathrm{s}<\ln\mathfrak{r}\}$ . This leads to a piecewise convergence boundary which combines the boundaries of the two regions $\mathcal{B}$ and $\mathcal{C}$ . Only the boundary $\partial\mathcal{B}$ is shown in Figure 2 since $\partial\mathcal{C}=\{\mathrm{s}=\ln\mathfrak{r}\}$ is a horizontal line which is out of the range of the graph. The two boundaries meet at $(r_{0},s_{0})=({-}3.91,0.103)=\partial\mathcal{B}\cap\partial\mathcal{C}$ , so the boundary is flat for $\mathrm{r}<-3.91$ , and for $\mathrm{r}\geq-3.91$ it is concave and determined by $\partial\mathcal{B}$ . The proof of the next result is given in the Supplementary Material.

3.2. Residue-saddlepoint approximations

The second approach to the inversion of $\mathcal{R}(\mathrm{y},\mathrm{z})$ is to perform the inversion in two steps. First, we approximately invert $\mathcal{R}(\mathrm{y},\mathrm{z})$ in $\mathrm{z} \leftrightsquigarrow n$ using a residue expansion as in Lemma 1 to get an approximation $\mathcal{\hat{R}}_{n}(\mathrm{y})$ for $\mathcal{R}_{n}(\mathrm{y})$ . Then, using $\mathcal{\hat{R}}_{n}(\mathrm{y})$ as a surrogate for the true PGF of R(n), we apply a single-saddlepoint approximation to the approximate CGF $\ln\mathcal{\hat{R}}_{n}({\mathrm{e}}^{\mathrm{r}})$ to determine what we call a residue-saddlepoint approximation. There are two such single-saddlepoint approximations, which we denote as $\mathbb{P}_{R1}$ and $\mathbb{P}_{R2}$ , that relate to the continuity correction applied, as presented in the Supplementary Material.

For the first stage, the residue approximation $\mathcal{\hat{R}}_{n}(\mathrm{y})$ for $\mathcal{R}_{n}(\mathrm{y})$ extends the development in [Reference Butler6, Section 3] by using residue expansions developed in [Reference Butler5, Corollary 1], which we summarise as follows.

Lemma 1. Consider the generating function $\mathcal{A}(\mathrm{z})=\sum_{n=0}^{\infty}a_{n}\mathrm{z}^{n}$ with $a_{n}\geq0$ for all n. Assume $\mathcal{A}(\mathrm{z})$ is analytic on $\{\mathrm{z}\in\mathbb{C}\,:\,|\mathrm{z}|\leq\mathfrak{r}\}$ apart from a simple pole at $\mathfrak{r}\in(0,\infty)$ . Then, for some $\varepsilon>0$ ,

(7) \begin{equation} a_{n} = -\frac{\xi_{-1}}{\mathfrak{r}^{n+1}} + o\{(\mathfrak{r}+\varepsilon)^{-n}\}, \quad n\rightarrow\infty, \end{equation}

where $\xi_{-1}=\lim_{\mathrm{z}\rightarrow\mathfrak{r}}\{(\mathrm{z}-\mathfrak{r})\mathcal{A}(\mathrm{z})\}$ is the residue of $\mathcal{A}$ at $\mathfrak{r}$ .

The value of $\varepsilon$ is limited by any singularity $\mathfrak{r}_{2}\in\mathbb{C\cup\{\infty\}}$ of $\mathcal{A}(\mathrm{z})$ which has the next largest modulus, say $|\mathfrak{r}_{2}|$ . Then $\varepsilon$ may be any value in $(0,|\mathfrak{r}_{2}|-\mathfrak{r)}$ . If no other singularities exist (including at $\infty$ ), then the approximation in (7) is exact. For a PGF, this means that $\mathcal{A}(\mathrm{z})=(1-1/\mathfrak{r})/(1-\mathrm{z}/\mathfrak{r}).$

Suppose $\mathcal{R}({\mathrm{e}}^{\mathrm{r}},{\mathrm{e}}^{\mathrm{s}})$ has a convergence region as described for Figure 2 and fix ${\mathrm{e}}^{\mathrm{r}}=\mathrm{y}>y_{0}={\mathrm{e}}^{r_{0}}$ . If $\hat{z}(\mathrm{y})$ denotes the smallest-modulus pole of $\mathcal{R}(\mathrm{y},\mathrm{z})$ in z for fixed y, then $(\mathrm{y},\hat{z}(\mathrm{y}))$ lies on the convergence boundary of $\mathcal{R}(\mathrm{y},\mathrm{z})$ and $\hat{z}(\mathrm{y})$ is the smallest positive root of $1-\mathcal{H}(\mathrm{y},\mathrm{z})=0$ and a dominant pole for $\mathcal{R}(\mathrm{y},\mathrm{z})$ . Thus, the range of allowable values of y for using Lemma 1 is $\mathrm{y}\in(y_{0},y_{1})$ , where $y_{0}={\mathrm{e}}^{r_{0}}$ and $y_{1}$ is the largest value of y for which $(\mathrm{y},\hat{z}(\mathrm{y}))$ remains on the convergence boundary of $\mathcal{R}(\mathrm{y},\mathrm{z})$ .

Theorem 2. ( $\mathcal{\hat{R}}_{n}(\mathrm{y})$ ). Assume $\mathcal{H}(\mathrm{y},\mathrm{z})$ has a maximal convergence region which is an open neighbourhood of (1, 1), and $\mathcal{F}(\mathrm{z})=\mathcal{H}(1,\mathrm{z})$ is aperiodic with convergence region $\{\mathrm{z}\in\mathbb{C}\,:\,|\mathrm{z}|<\mathfrak{r}\}$ for $\mathfrak{r}>1$ . A residue approximation for $\mathcal{R}_{n}(\mathrm{y})$ is

\begin{align*} \mathcal{\hat{R}}_{n}(\mathrm{y}) = \frac{1}{\hat{z}(\text{y})^{n+1}} \frac{\mathcal{F}_{S}\{\hat{z}(\text{y})\}}{\mathcal{H}^{\prime}_{\mathrm{z}}\{\mathrm{y},\hat{z}(\text{y})\}}, \quad \mathrm{y}\in(y_{0},y_{1}). \end{align*}

The order of the expansion as $n\rightarrow\infty$ is $o\{(\hat{z}(\text{y})+\varepsilon)^{-n}\}$ for any $0<\varepsilon<\min\{\mathfrak{r}-\hat{z}(\mathrm{y}),|\hat{z}_{2}(\mathrm{y})|-\hat{z}(\mathrm{y})\}$ , where $|\hat{z}_{2}(\mathrm{y})|$ is the modulus of the second largest zero of $1-\mathcal{H}(\mathrm{y},\mathrm{z})$ , if such a zero exists, or $\infty$ otherwise.

Proof. While holding $\mathrm{y}\in(y_{0},y_{1})$ fixed, $\mathcal{H}(\mathrm{y},\mathrm{z})$ is monotone increasing from 0 to $\infty$ as z increases from 0 until $(\mathrm{y},\mathrm{z})$ approaches the convergence boundary of $\mathcal{R}(\mathrm{y},\mathrm{z})$ . Thus, a smallest positive real root $\hat{z}(\mathrm{y})$ exists. To show that $\hat{z}(\mathrm{y})$ is the dominant pole of $\mathcal{R}(\mathrm{y},\mathrm{z})$ as in Lemma 1, we show that it is the unique zero of $1-\mathcal{H}(\mathrm{y},\mathrm{z})$ with smallest modulus. It suffices to show that, for $\mathrm{z}\in(0,\hat{z}(\mathrm{y})]$ ,

(8) \begin{equation} |\mathcal{H}(\mathrm{y},\mathrm{z}{\mathrm{e}}^{\mathrm{i}\theta})| < \mathcal{H}(\mathrm{y},\mathrm{z}), \quad \theta\neq0,\ (\mathrm{y},\mathrm{z})\in\mathcal{N}. \end{equation}

Then $ |\mathcal{H}(\mathrm{y},\mathrm{z}{\mathrm{e}}^{\mathrm{i}\theta})| < \mathcal{H}(\mathrm{y},\mathrm{z}) \leq \mathcal{H}\{\mathrm{y},\hat{z}(y)\} = 1, $ with equality only occurring when $\mathrm{z}=\hat{z}(y)$ and $\theta=0$ . Inequality (8) is proved in the Supplementary Material. It is the strict inequality in (8) that matters, as the non-strict inequality $(\leq)$ for (8) easily follows from the triangle inequality (see the Supplementary Material).

From (8), we see that as z increases from 0 for fixed y, the value of $|\mathcal{H}(\mathrm{y},\mathrm{z}{\mathrm{e}}^{\mathrm{i}\theta})|$ increases the fastest in the direction $\theta=0$ . Hence, for y fixed, $\mathcal{H}(\mathrm{y},\mathrm{z})$ is analytic on $\{\mathrm{z}\in\mathbb{C}\,:\,|\mathrm{z}|\leq\hat{z}(\mathrm{y})\}$ apart from the pole at $\hat{z}(\mathrm{y})$ , as required by Lemma 1. The pole is simple since $\mathcal{H}^{\prime}_{\mathrm{z}}\{\mathrm{y},\hat{z}(\mathrm{y})\}>0$ . The order of the expansion follows from the discussion of order for Lemma 1.

The assumption that $\mathcal{H}(\mathrm{y},\mathrm{z})$ has an open convergence domain can be relaxed as long as $\hat{z}(\mathrm{y})$ is uniquely defined for each y and the value of y is appropriately restricted.

Example 6. (Viral RNA tests.) The value $\hat{z}(\mathrm{y})$ which solves $\mathcal{H}(\mathrm{y},\mathrm{z})=1$ in z using (4) is $\hat{z}(\mathrm{y})=1/B(\mathrm{y})$ , and

\begin{equation*} \mathcal{\hat{R}}_{n}(\mathrm{y}) = B(\mathrm{y}\mathbf{)}^{n} \frac{B_{1}(\mathrm{y})}{B(\mathrm{y})-p_{0\cdot}}, \quad \mathrm{y}>y_{0}=0.02, \end{equation*}

with $y_{1}=\infty$ . This replicates the dominant term in the partial fraction expansion (note that $p_{0\cdot}^{n}$ is small) in (6). We now use $\mathcal{\hat{K}}_{n}(\mathrm{r})=\ln\mathcal{\hat{R}}_{n}({\mathrm{e}}^{\mathrm{r}})$ as a surrogate for the true CGF $\mathcal{K}_{n}(\mathrm{r})=\ln\mathcal{R}_{n}({\mathrm{e}}^{\mathrm{r}})$ for R(n). Table 2 displays the results of this strategy where $\mathbb{P}_{R1}$ and $\mathbb{P}_{R1}$ are continuity-corrected single-saddlepoint approximations as presented in [Reference Butler3, Section 1.2.3]. To approximate $\mathbb{P}\{R(n)\geq m\}$ , $\mathbb{P}_{R1}$ uses saddlepoint $\hat{r}$ solving $\mathcal{\hat{K}}^{\prime}_{n}(\mathrm{r})=m$ and $\mathbb{P}_{R2}$ uses saddlepoint $\tilde{r}$ solving $\mathcal{\hat{K}}^{\prime}_{n}(\mathrm{r})=m-\frac12$ . Table 2, as well as the percentage relative error plots in Figure 1, show that $\mathbb{P}_{R1}$ (solid red line) and $\mathbb{P}_{R2}$ (dashed blue line) achieve remarkable accuracy.

At $\mathrm{y}=y_{0}$ , $B(y_{0})-p_{0\cdot}=0$ and $\mathcal{\hat{R}}_{n}(y_{0})$ is undefined. For $\mathrm{y}\in(0,y_{0})$ , $\mathcal{\hat{R}}_{n}(\mathrm{y})<0$ but in the true PGF expansion in (6) with both terms, $\mathcal{R}_{n}(\mathrm{y})>0$ . This limitation in the domain for $\mathcal{\hat{R}}_{n}(\mathrm{y})$ had no practical relevance in this or the other examples.

3.3. Continuous rewards

For continuously valued rewards $R\geq0$ , the distribution of R(n) is a mixture distribution which puts a point mass at $x=0$ with probability $\mathbb{P}\{F>n\}$ and has a continuous density over $(0,\infty)$ with probability $\mathbb{P}\{F\leq n\}$ . In order to best use saddlepoint and residue approximations, the point mass should be removed so we are dealing only with the continuous portion of the distribution.

We assume that R and F have a mixed MGF-PGF $\mathcal{H}_{\mathrm{c}}(\mathrm{r},\mathrm{z})=\mathbb{E}\{{\mathrm{e}}^{\mathrm{r}R}\mathrm{z}^{F}\}$ , and the maximal convergence region is an open neighbourhood of $(\mathrm{r},\mathrm{z})=(0,1)$ . The portion of the double transform of reward R(n) associated with the point mass is

\begin{align*} \sum_{n=0}^{\infty}\mathbb{E}\big\{{\mathrm{e}}^{\mathrm{r}R(n)}\mathbf{1}_{\{R(n)=0\}}\big\}\mathrm{z}^{n} = \sum_{n=0}^{\infty}\mathbb{P}\{F>n\}\mathrm{z}^{n} = \mathcal{F}_{S}(\mathrm{z}).\end{align*}

Removing this from $\mathcal{R}_{\mathrm{c}}(\mathrm{r},\mathrm{z})$ gives

(9) \begin{equation} \mathcal{R}_{\mathrm{c}0}(\mathrm{r},\mathrm{z}) = \sum_{n=0}^{\infty}\mathbb{E}\big\{{\mathrm{e}}^{\mathrm{r}R(n)}\mathbf{1}_{\{R(n)>0\}}\big\}\mathrm{z}^{n} = \mathcal{R}_{\mathrm{c}}(\mathrm{r},\mathrm{z})-\mathcal{F}_{S}(\mathrm{z}) = \frac{\mathcal{F}_{S}(\mathrm{z})\mathcal{H}_{\mathrm{c}}(\mathrm{r},\mathrm{z})} {1-\mathcal{H}_{\mathrm{c}}(\mathrm{r},\mathrm{z})}, \end{equation}

and the conditional MGF of R(n) given $R(n)>0$ is $[\mathrm{z}^{n}]\mathcal{R}_{\mathrm{c}0}(\mathrm{r},\mathrm{z})/c_{n}$ , with $c_{n}=\mathbb{P}\{F\leq n\}$ .

To approximate the marginal probability $\mathbb{P}\{R(n)\geq x\}$ for $x>0$ , we may apply the continuous version of the Skovgaard double-saddlepoint approximation $\mathbb{P}_{\mathrm{C}}$ in (54) as outlined in the Supplementary Material. As noted in [Reference Butler6, Section 6.1], we may take the joint CGF as $\mathcal{K}(\mathrm{r},\mathrm{s})=\ln\{\mathcal{R}_{\mathrm{c}0}(\mathrm{r},{\mathrm{e}}^{\mathrm{s}})\}$ and the approximation $\mathbb{P}_{\mathrm{C}}$ is for $\mathbb{P}\{R(n)\geq x\mid R(n)>0\}$ . This needs to be scaled by $c_{n}=\mathbb{P}\{F\leq n\}$ to give the marginal survival probability as

\begin{align*} \mathbb{P}\{R(n)\geq x\}\approx c_{n}\mathbb{P}_{\mathrm{C}}\{R(T)\geq x\mid T=n\}, \quad x>0.\end{align*}

The joint CGF $\mathcal{K}(\mathrm{r},\mathrm{s})=\ln\{\mathcal{R}_{\mathrm{c}0}(\mathrm{r},{\mathrm{e}}^{\mathrm{s}})\}$ is convergent on a region, as seen in Figure 2.

For the residue-saddlepoint approximation, the first stage gives an approximate marginal MGF for R(n). We apply Lemma 1 to (9) and take $\hat{z}_{\mathrm{c}}(\textit{r})$ as the smallest positive zero of $1-\mathcal{H}_{\mathrm{c}}(\mathrm{r},\mathrm{z})$ in z. Then, using $\mathcal{H}_{\mathrm{c}}\{\mathrm{r},\hat{z}_{\mathrm{c}}(\mathrm{r})\}=1$ ,

(10) \begin{equation} \mathcal{\hat{R}}_{\mathrm{c}n}(\mathrm{r}) = \frac{1}{\hat{z}_{\mathrm{c}}(\mathrm{r})^{n+1}} \frac{\mathcal{F}_{S}\{\hat{z}_{\mathrm{c}}(\mathrm{r})\}} {\mathcal{H}_{\mathrm{c}\mathrm{z}}^{\prime}\{\mathrm{r},\hat{z}_{\mathrm{c}}(\mathrm{r})\}}, \quad \mathrm{r}\in(r_{0},r_{1}). \end{equation}

The range for r has $r_{0}$ such that $(r_{0},\mathfrak{r})$ is the junction point in the boundary of convergence for $\mathcal{R}_{\mathrm{c}}(\mathrm{r},\mathrm{z})$ , i.e. $1-\mathcal{H}_{\mathrm{c}}(r_{0},\mathfrak{r})=0$ , where $\mathfrak{r}$ is the convergence radius of $\mathcal{F}_{S}$ . The upper range $r_{1}$ is the supremum of values r such that $\mathcal{H}_{\mathrm{c}}(\mathrm{r},\hat{z}_{\mathrm{c}}(\mathrm{r}))=1$ so that $(\mathrm{r},\hat{z}_{\mathrm{c}}(\mathrm{r}))$ remains on the convergence boundary of $\mathcal{R}_{\mathrm{c}}(\mathrm{r},\mathrm{z})$ .

At the second stage, $\mathcal{\hat{R}}_{\mathrm{c}n}(\mathrm{r})$ is used as a surrogate for the true MGF $\mathcal{R}_{n}(\mathrm{r})=\mathbb{E}\{{\mathrm{e}}^{\mathrm{r}R(n)}\}$ in a single-saddlepoint approximation of the Lugannani–Rice [Reference Lugannani and Rice15] type as given [Reference Borovkov1, Section 1.2.1]. No rescaling by $c_{n}$ is required here as with the double saddlepoint method.

Computations for continuous rewards can be found in Example 13, and in Example 15 of the Supplementary Material.

4.1. Uniform expansion and central limit theorem

The CGF for R(n) is $\mathcal{K}_{n}(\mathrm{r})=\ln\mathcal{R}_{n}({\mathrm{e}}^{\mathrm{r}})$ , which is approximated by using $\mathcal{\hat{K}}_{n}(\mathrm{r})=\ln\mathcal{\hat{R}}_{n}({\mathrm{e}}^{\mathrm{r}})$ . The residue approximation is uniformly $o\{(1+\varepsilon)^{-n}\}$ as $n\rightarrow\infty$ over compact subsets of r, as stated in Corollary 3. Building upon this, we show that the standardised value of R(n) using the moments of Corollary 2 has a CLT. Since the distribution of R(n) tends to be quite skewed even for moderately large n, these normal limits only achieve reasonable accuracy when n is quite large.

Corollary 3. (Uniform asymptotic expansion.) For any compact interval $\mathcal{D}$ on which $\mathcal{\hat{K}}_{n}(\mathrm{r})$ is well-defined, there exists an $\varepsilon>0$ such that

(13) \begin{equation} \max_{\mathrm{r}\in\mathcal{D}}|\mathcal{K}_{n}(\mathrm{r})-\mathcal{\hat{K}}_{n}(\mathrm{r})| = o\{(1+\varepsilon)^{-n}\}, \quad n\rightarrow\infty. \end{equation}

A central limit is now derived for $Z_{n}=[R(n)-\mathbb{E}\{R(n)\}]/\sqrt{\mathbb{V}\{R(n)\}}$ , the standardised value of R(n), with $\mathbb{E}\{R(n)\}$ and $\mathbb{V}\{R(n)\}$ as given in Corollary 2.

6.1. The double GF of $\{P_{m}\,:\, m\geq0\}$

The PGF of $P_{m}$ is $\mathcal{P}_{m}(\mathrm{z})=\sum_{n=1}^{\infty}\mathbb{P}\{P_{m}=n\}\mathrm{z}^{n}$ for $m\geq1$ . For $m=0$ , $P_{0}=0$ with probability 1. The double GF in variables $(\mathrm{y},\mathrm{z})\leftrightsquigarrow(m,n)$ is

\begin{align*} \mathcal{P}(\mathrm{y},\mathrm{z})=\sum_{m=0}^{\infty}\mathcal{P}_{m}(\mathrm{z})\mathrm{y}^{m}, \quad |\mathrm{y}|<1,\ |\mathrm{z}|<1.\end{align*}

Let $\mathcal{E}(\mathrm{y})=\mathcal{H}(\mathrm{y},1)$ be the reward PGF and write $\mathcal{E}_{\mathcal{S}}(\mathrm{y})=\{1-\mathcal{E}(\mathrm{y})\}/(1-\mathrm{y})$ .

Theorem 3. (GF for the first-passage reward mass function.) Suppose (R,F) has a joint PGF $\mathcal{H}(\mathrm{y},\mathrm{z})$ which is convergent in a neighbourhood of (1, 1). Then

(16) \begin{equation} \mathcal{P}(\mathrm{y},\mathrm{z}) = \frac{1}{1-\mathrm{y}}\bigg[1-\frac{\mathrm{y}\{1-\mathcal{F}(\mathrm{z})\}} {1-\mathcal{H}(\mathrm{y},\mathrm{z})}\bigg] , \end{equation}

so that $\mathcal{P}_{m}(\mathrm{z})=[\mathrm{y}^{m}]\mathcal{P}(\mathrm{y},\mathrm{z})$ and $\mathbb{P}\{P_{m}=n\}=[\mathrm{y}^{m}\mathrm{z}^{n}]\mathcal{P}(\mathrm{y},\mathrm{z})$ . Furthermore,

\begin{align*} \mathbb{P}\{P_{m}\geq n\} & = [\mathrm{y}^{m-1}\mathrm{z}^{n-1}]\bigg\{\frac{\mathcal{F}_{\mathcal{S}}(\mathrm{z})} {(1-\mathrm{y})\{1-\mathcal{H}(\mathrm{y},\mathrm{z})\}}\bigg\}, \\ \mathbb{E}\{P_{m}\} & = [\mathrm{y}^{m-1}]\bigg\{\frac{\mu}{(1-\mathrm{y})^{2}\mathcal{E}_{\mathcal{S}}(\mathrm{y})\}}\bigg\}. \end{align*}

Note that $\mathrm{y}=1$ is a removable singularity of $\mathcal{P}(\mathrm{y},\mathrm{z})$ in (16) for values of z in the convergence domain of $\mathcal{F}(\mathrm{z})$ .

Proof. Starting with (15), compute the GF in $n\leftrightsquigarrow\mathrm{z}$ of the left side and the GF in $m\leftrightsquigarrow\mathrm{y}$ of the right side expressed as $\mathbb{P}\{R(n)\leq m-1\}$ . The equality in (15) says that the respective inversions of these transforms are equal and that

\begin{align*} [\mathrm{z}^{n}]\bigg\{\frac{1-\mathcal{P}_{m}(\mathrm{z})}{1-\mathrm{z}}\bigg\} = [\mathrm{y}^{m-1}]\bigg\{\frac{\mathcal{R}_{n}(\mathrm{y})}{1-\mathrm{y}}\bigg\} = [\mathrm{y}^{m}]\bigg\{\frac{\mathrm{y}\mathcal{R}_{n}(\mathrm{y})}{1-\mathrm{y}}\bigg\} , \end{align*}

where $\mathcal{R}_{n}(\mathrm{y})$ is the PGF of R(n). Now take the GF in $m\leftrightsquigarrow\mathrm{y}$ of the left side and the GF in $n\leftrightsquigarrow\mathrm{z}$ on the right side to get

\begin{align*} [\mathrm{y}^{m}\mathrm{z}^{n}]\bigg\{\frac{1}{(1-\mathrm{y})(1-\mathrm{z})} - \frac{\mathcal{P}(\mathrm{y},\mathrm{z})}{1-\mathrm{z}}\bigg\} = [\mathrm{z}^{n}\mathrm{y}^{m}]\bigg\{\frac{\mathrm{y}}{1-\mathrm{y}} \frac{\mathcal{F}_{\mathcal{S}}(\mathrm{z})}{1-\mathcal{H}(\mathrm{y},\mathrm{z})}\bigg\}, \end{align*}

where

\begin{align*} \mathcal{R}(\mathrm{y},\mathrm{z})=\sum_{n=0}^{\infty}\mathcal{R} _{n}(\mathrm{y})\mathrm{z}^{n}=\frac{\mathcal{F}_{\mathcal{S}}(\mathrm{z} )}{1-\mathcal{H}(\mathrm{y},\mathrm{z})} \end{align*}

has been used on the right-hand side. Equating the two expressions in curly braces and solving for $\mathcal{P}(\mathrm{y},\mathrm{z})$ leads to (16).

To derive the survival function double GF in (16), start with

(17) \begin{align} \mathbb{P}\{P_{m}\geq n\} = [\mathrm{z}^{n}]\bigg\{\frac{1-\mathrm{z}\mathcal{P}_{m}(\mathrm{z})}{1-\mathrm{z}}\bigg\} & = [\mathrm{z}^{n}]\bigg\{\frac{1-\mathrm{z}[\mathrm{y}^{m}]\mathcal{P}(\mathrm{y},\mathrm{z})} {1-\mathrm{z}}\bigg\} \nonumber \\ & = [\mathrm{y}^{m}\mathrm{z}^{n}]\bigg\{\frac{1}{(1-\mathrm{y})(1-\mathrm{z})} - \frac{\mathrm{z}\mathcal{P}(\mathrm{y},\mathrm{z})}{1-\mathrm{z}}\bigg\} \nonumber \\ & = [\mathrm{y}^{m}\mathrm{z}^{n}]\bigg\{\frac{1}{1-\mathrm{y}} + \frac{\mathrm{yz}\{1-\mathcal{F}(\mathrm{z})\}} {(1-\mathrm{y})(1-\mathrm{z})\{1-\mathcal{H}(\mathrm{y},\mathrm{z})\}}\bigg\} \nonumber \\ & = [\mathrm{y}^{m-1}\mathrm{z}^{n-1}]\bigg\{\frac{\mathcal{F}_{S}(\mathrm{z})} {(1-\mathrm{y})\{1-\mathcal{H}(\mathrm{y},\mathrm{z})\}}\bigg\}. \end{align}

To get the third line we have substituted the value of $\mathcal{P}(\mathrm{y},\mathrm{z})$ . For the fourth line, note that $[\mathrm{z}^{n}]\{1/(1-\mathrm{y})\}=0$ .

Since $\mathbb{E}\{P_{m}\}=\sum_{n=1}^{\infty}\mathbb{P}\{P_{m}\geq n\}$ , the inverse in $\mathrm{y}\leftrightsquigarrow m-1$ of the double GF in (17) when evaluated at $\mathrm{z}=1$ gives the GF for the sequence of means as

\begin{equation*} \sum^{\infty}_{n=0} \mathbb{P}\{P_x \geq n\} -1 = [\mathrm{y}^{m}]\left\{\frac{1}{1-y}+\frac{\mathrm{y}\mathcal{F}_{S}(1)}{(1-\mathrm{y})\{1-\mathcal{H}(\mathrm{y},1)\}}\right\} -1 = [\mathrm{y}^{m-1}]\left\{\frac{\mu}{(1-\mathrm{y})^ {2}\mathcal{E}_{\mathcal{S}}(\mathrm{y})}\right\}. \end{equation*}

6.2. Asymptotic expansions

Expansions for the moments and distribution of $P_{m}$ as $m\rightarrow\infty$ are given below.

Corollary 7. (Mean and variance expansions for first-passage reward moments.) If (R,F) has a joint PGF $\mathcal{H}(\mathrm{y},\mathrm{z})$ which is convergent in a neighbourhood of (1, 1), then

(18) \begin{align} \mathbb{E}\{P_{m}\} & = \bigg(m-\frac12\bigg)\frac{\mu}{\rho} + \frac{\mu}{2\rho^{2}}\mathbb{E}\{R^{2}\} + o(1), \nonumber \\ \mathbb{V}\{P_{m}\} & = m\frac{\mu^{2}}{\rho^{3}}\mathbb{E}\bigg\{R-\frac{\rho}{\mu}F\bigg\}^{2} + O(1) = m\bigg(\frac{\mu}{\rho}\bigg)^{3}\sigma_{R.F}^{2} + O(1), \end{align}

as $m\rightarrow\infty$ .

The passage time $P_{m}$ has a limiting normal distribution as $m\rightarrow\infty$ after standardisation using the mean and variance of Corollary 7 to order O(1). The proof is a standard result and given in the Supplementary Material.

Corollary 8. (Central limit for $P_{m}$ .) Subject to the conditions of Corollary 7,

\begin{align*} \frac{P_{m}-\mathbb{E}\{P_{m}\}}{\sqrt{\mathbb{V}\{P_{m}\}}} \quad \overset{w}{\longrightarrow} \quad \mathrm{Normal}(0,1)\quad \mathrm{as}\ m\rightarrow \infty. \end{align*}

The convergence in Corollary 8 tends to be very slow as $P_{m}$ has a skewed distribution even for moderately large m.

Example 8. (Bivariate geometric rewards.) Expansion of $\mathcal{P}(\mathrm{y},\mathrm{z})$ in powers of y provides the first few PGFs for $P_{0},\ldots,P_{3}$ . With $q_{i}=1-p_{i}$ , then

(19) \begin{equation} \mathcal{P}(\mathrm{y},\mathrm{z}) = 1 + \frac{p_{2}\mathrm{z}}{1-q_{2}\mathrm{z}}(\mathrm{y}+\mathrm{y}^{2}) + \frac{q_{0}\mathrm{z}}{1-p_{0}\mathrm{z}}\frac{p_{2}\mathrm{z}}{1-q_{2}\mathrm{z}}\mathrm{y}^{3} + O(\mathrm{y}^{4}). \end{equation}

From this, we see that $P_{0}=1$ with probability 1, and $P_{1}$ and $P_{2}$ are Geometric $(p_{2})$ . The stopping time for a reward of at least one or two only happens when a renewal occurs in which the two fail simultaneously, hence this is a Geometric $(p_{2})$ waiting time as seen in (19). As a check, take $p_{0}=\frac{7}{10}$ , $p_{1}=\frac{2}{10}$ , and $p_{2}=\frac{1}{10}$ . The geometric probability $\mathbb{P}\{P_{2}=2\}=p_{2}q_{2}=\frac{1}{10}\times\frac{9}{10}=\frac{9}{100}$ and occurs in two steps with outcomes $\{02\}\cup\{12\}$ , where $\{12\}$ means one failure followed by two failures in the first two ticks of the clock. These probabilities add to $\frac{7}{100}+\frac{2}{100}=\frac{9}{100}$ , which agrees.

The stopping time for three or more rewards is $P_{3}\sim \mathrm{Geometric}(q_{0}) + \mathrm{Geometric}(p_{2})$ , with the two geometric distributions independent. This amounts to first waiting for one or two failures in Geometric $(q_{0})$ time followed by Geometric $(p_{2})$ time for a renewal event to realise the reward. To check this numerically, we compute the convolution as

\begin{align*} \mathbb{P}\{P_{3}=3\} = \sum_{k=1}^{2}\bigg(\frac{7}{10}\bigg)^{k-1}\frac{3}{10} \bigg(\frac{9}{10}\bigg)^{3-k-1}\frac{1}{10} = \frac{6}{125}. \end{align*}

Since $\{P_{3}=3\}=\{012\}\cup\{102\}\cup\{202\}\cup\{022\}\cup\{112\}\cup\{212\}$ , summing these mutually exclusive outcomes gives

\begin{align*} \mathbb{P}\{P_{3}=3\} = 2\bigg(\frac{14}{1000}\bigg) + 2\bigg(\frac{7}{1000}\bigg) + \frac{4}{1000} + \frac{2}{1000} = \frac{6}{125}. \end{align*}

The expansion terms of Corollary 7 and the true moments for $P_{5}$ are

\begin{align*} 18.148 = \mathbb{E}\{P_{5}\} = 18.125 + o(1), \qquad 111.76 = \mathbb{V}\{P_{5}\} = 34.375 + O(1). \end{align*}

While the mean expansion is quite accurate, the variance expansion of order O(m) is very inaccurate. This is partly due to a small value of $m=5$ , but it also suggests that the quite complicated term of order O(1) which has been left out is large and most likely even larger than $34.375$ . Without a more accurate variance approximation, the normal approximation would be very inaccurate and misleading. The double-saddlepoint and residue-saddlepoint methods of the next section provide much greater accuracy.

7.1. Residue-saddlepoint inversion using $\mathcal{R}(\mathrm{y},\mathrm{z})$

Inversion of the right-hand side of (20) starts with the double GF for mass function $\mathbb{P}\{R(n)=m\}$ in (2) of Theorem 1. If $\hat{z}(\mathrm{y})$ is the smallest positive root of $1-\mathcal{H}(\mathrm{y},\mathrm{z})=0$ , then the PGF of R(n) has residue approximation

(21) \begin{equation} \mathcal{\hat{R}}_{n}(\mathrm{y}) = \frac{1}{\hat{z}(\mathrm{y})^{n+1}} \frac{\mathcal{F}_{S}\{\hat{z}(\mathrm{y})\}}{\mathcal{H}^{\prime}_{\mathrm{z}}\{\mathrm{y},\hat{z}(\mathrm{y})\}}, \quad \mathrm{y}\in(y_{0},y_{1}), \end{equation}

where $0\leq y_{0}<1<y_{1}$ . The value $y_{0}$ solves $\hat{z}(y_{0})=\mathfrak{r}$ and the value $y_{1}$ is the supremum of the values of y which solve $1-\mathcal{H}\{\mathrm{y},\hat{z}(\mathrm{y})\}=0$ . Substituting $\mathrm{y}={\mathrm{e}}^{\mathrm{r}}$ , the approximate CGF for R(n) is $\mathcal{\hat{K}}_{n}(\mathrm{r})=-(n+1)\ln\{\hat{z}({\mathrm{e}}^{\mathrm{r}})\}+l(\mathrm{r})$ for

\begin{align*} l(\mathrm{r}) = \ln\frac{\mathcal{F}_{S}\{\hat{z}({\mathrm{e}}^{\mathrm{r}})\}} {\mathcal{H}^{\prime}_{\mathrm{z}}\{{\mathrm{e}}^{\mathrm{r}},\hat{z}({\mathrm{e}}^{\mathrm{r}})\}},\end{align*}

where $\mathrm{r}\in(r_{0},r_{1})$ with $r_{0}=\ln y_{0}^{+}<0<\ln y_{1}=r_{1}$ . To compute $\mathbb{P}_{R1}$ , the saddlepoint $\hat{r}=\hat{r}(m)$ solves the saddlepoint equation

(22) \begin{equation} m = \mathcal{\hat{K}}^{\prime}_{n}(\hat{r}) = -(n+1)\frac{\hat{z}^{\prime}({\mathrm{e}}^{\hat{r}}){\mathrm{e}}^{\hat{r}}}{\hat{z}({\mathrm{e}}^{\hat{r}})} + l^{\prime}(\hat{r}). \end{equation}

This determines the single-saddlepoint approximation $\mathbb{\hat{P}}_{R1}\{R(n)\geq m\}$ for the right-hand side of (20), which is the survival approximation of $P_{m}$ at $n+1$ .

A parametric plot in Maple and Mathematica avoids having to solve (22). For fixed m, the value for n can be expressed explicitly in terms of $\hat{r}$ by solving (22) so that $\hat{r}$ becomes the parametric index for the plot. If $n(\hat{r})$ solves (22) in n for given $\hat{r}$ , then the parametric plot is $(n(\hat{r})+1, 1-\mathbb{P}_{R1}\{R\{n(\hat{r})\}\geq m\})$ , $\hat{r}\in(r_{0},r_{1})$ . This provides a smooth approximation for the survival function of $P_{m}$ valid at integer arguments. The form of this parametric curve follows from the identity (20) in which the second component is the survival probability at $n(\hat{r})+1$ for $P_{m}$ .

A small value of $\hat{r}$ just above $r_{0}$ computes the right tail of $P_{m}$ , and values approaching $r_{1}$ compute the left tail. To see this, note that $\hat{z}(\mathrm{y})$ is decreasing in $\mathrm{y}={\mathrm{e}}^{\mathrm{r}}$ , and hence r. Thus, a small value of $\hat{r}$ reflects a large value of $n(\hat{r})$ or the right tail for $P_{m}$ and the left tail of $R\{n(\hat{r})\}$ .

The second continuity-corrected saddlepoint approximation is much the same except that one solves saddlepoint equation (22) with m offset to $m^{-}=m-\frac12$ . With saddlepoint $\tilde{r}$ , and $n^{-}(\tilde{r})$ the solution using $m^{-}$ , the parametric plot is $(n^{-}(\tilde{r})+1, 1-\mathbb{P}_{R2}\{R\{n^{-}(\tilde{r})\}\geq m\})$ , $\tilde{r}\in(r_{0},r_{1})$ .

Example 9. (Bivariate geometric rewards.) Solving for $\hat{z}(\mathrm{y})$ leads to

\begin{align*} \hat{z}(\mathrm{y}) = \frac{1}{p_{0}+p_{1}\mathrm{y}+p_{2}\mathrm{y}^{2}} =: \frac{1}{D(\mathrm{y})}, \end{align*}

where $D(\mathrm{y})$ is the GF for the number of failures with each tick of the clock. Since $\mathcal{F}_{S}(\mathrm{z})=1/(1-9\mathrm{z}/10)$ ,

\begin{align*} \big\{\mathrm{y}>0\,:\,\hat{z}(\mathrm{y})<\tfrac{10}{9}\big\} = \big\{\mathrm{y}>0\,:\, D(\mathrm{y}) > \tfrac{9}{10}\big\} = \{\mathrm{y}>0.732\}, \end{align*}

so $y_{0}=0.732$ . As $\mathrm{y}\rightarrow\infty$ , $\hat{z}(\mathrm{y})\rightarrow1/p_{0}$ so $y_{1}=\infty$ . This leads to the unbounded saddlepoint range $\hat{r}>\ln(0.732)=-0.312=r_{0}$ as the index range of the parametric plot. The $\mathbb{P}_{R1}$ survival approximation of $P_{5}$ is given in the left panel of Figure 3 (solid red line) with exact values shown as circles. The $\mathbb{P}_{R2}$ approximation is overlaid (dashed blue line) and is almost indistinguishable. Percentage relative errors are shown with matching line types in the right panel, showing under $-2.0\%$ and under $-3.5\%$ respectively above the median. As n increases, the residue approximation $\mathcal{\hat{R}}_{n}(\mathrm{y})$ approaches $\mathcal{R}_{n}(\mathrm{y})$ for each y, which helps to explain its greater accuracy for larger n.

Exact computations were based on inverting $\mathcal{P}(\mathrm{y},\mathrm{z})$ using symbolic Taylor expansion. The double GF $\mathcal{P}(\mathrm{y},\mathrm{z})$ was expanded in z to extract $\mathcal{P}_{5}(\mathrm{y})$ as the coefficient of $\mathrm{z}^{5}$ . This was followed by Taylor expansion of the GF $\{1-\mathrm{y}\mathcal{P}_{5}(\mathrm{y})\}/(1-\mathrm{y})$ in y. This is only feasible due to the small value of $m=5$ .

Figure 3. Left: Approximations for the survival function of $P_{5}$ by inverting $\mathcal{R}(\mathrm{y},\mathrm{z})$ . The $\mathbb{P}_{R1}$ (solid red line) and $\mathbb{P}_{R2}$ (dashed blue line) approximations are shown with exact values (circles) for odd integers of n. Right: Percentage relative errors respecting the tails are shown for the $\mathbb{P}_{R1}$ and $\mathbb{P}_{R2}$ approximations using matching line types. Also shown are the errors for the double-saddlepoint $\mathbb{P}_{D1}$ (red crosses) and $\mathbb{P}_{D2}$ (blue diamonds) approximations when inverting $\mathcal{R}(\mathrm{y},\mathrm{z})$ .

Example 10. (Reward as an interarrival random walk.) In the context of Example 4 in which the joint PGF takes the form $\mathcal{H}(\mathrm{y},\mathrm{z})=\mathcal{F}\{\mathrm{z}\mathcal{G}_{0}(\mathrm{y})\}$ , the residue expansion for the PGF of R(n) in (21) takes a particularly simple form. We assume that the convergence domains of both $\mathcal{F}(\mathrm{z})$ and $\mathcal{G}_{0}(\mathrm{y})$ are open neighbourhoods of 0 with radii $\mathfrak{r}>1$ and $\mathfrak{r}_{0}>1$ respectively. The value $\hat{z}(\mathrm{y})=1/\mathcal{G}_{0}(\mathrm{y})$ , and (21) reduces to

\begin{align*} \mathcal{\hat{R}}_{n}(\mathrm{y}) = \mathcal{G}_{0}(\mathrm{y})^{n}\frac{\mathcal{F}_{S}\{1/\mathcal{G}_{0}(\mathrm{y})\}}{\mu}, \quad \mathrm{y}\in(y_{0},\mathfrak{r}_{0}). \end{align*}

Here, $y_{0}<1$ solves $1/\mathcal{G}_{0}(y_{0})=\mathfrak{r}$ . For Example 1, $\mathcal{F}_{S}(\mathrm{z})=1=\mu$ so that $\mathcal{\hat{R}}_{n}(\mathrm{y})=(p\mathrm{y}+q)^{n}$ is exact.

7.2. Double-saddlepoint inversion of $\mathcal{R}(\mathrm{y},\mathrm{z})$

We use the double-saddlepoint method in Section 3.1 and treat

\begin{equation*} \mathcal{K}(\mathrm{r},\mathrm{s}) = \ln\mathcal{R}({\mathrm{e}}^{\mathrm{r}},{\mathrm{e}}^{\mathrm{s}}) = \ln\frac{\mathcal{F}_{S}({\mathrm{e}}^{\mathrm{s}})}{1-\mathcal{H}({\mathrm{e}}^{\mathrm{r}},{\mathrm{e}}^{\mathrm{s}})}, \quad (\mathrm{r},\mathrm{s})\in\mathcal{B}\cap\mathcal{C}, \end{equation*}

as the CGF for $\{R(T),T\}$ whose joint measure is $\mu\{R(T)=m,T=n\}=\mathbb{P}\{R(n)=m\}$ with $(m,n)\leftrightsquigarrow(\mathrm{y},\mathrm{z})$ . This CGF is well-defined over the open convex region $\mathcal{B}\cap\mathcal{C}$ , where $\mathcal{B}=\{(\mathrm{r},\mathrm{s})\,:\,\mathcal{H}({\mathrm{e}}^{\mathrm{r}},{\mathrm{e}}^{\mathrm{s}})<1\}$ and $\mathcal{C}=\{(\mathrm{r},\mathrm{s})\,:\,\mathrm{s}<\ln\mathfrak{r}\}$ . Region $\mathcal{B}$ excludes the first quadrant and has the form in Figure 2. Then $\mathbb{P}\{P_{m}\geq n+1\}=1-\mathbb{P}\{R(n)\geq m\}=1-\mathbb{P}\{R(T)\geq m\mid T=n\}$ , where the last conditional probability is approximated to determine the survival function of $P_{m}$ at $n+1$ .

For Example 9, the $\mathbb{P}_{D1}$ approximation (red crosses) and the $\mathbb{P}_{D2}$ approximation (blue diamonds) are shown in Figure 4 compared with exact values (solid line). The percentage relative errors using matching symbol types are shown in the right panel of Figure 3. The residue-saddlepoint approximations are somewhat more accurate, but these double-saddlepoint approximations are still quite accurate.

Figure 4. Double-saddlepoint approximations for the survival function of $P_{5}$ by inverting $\mathcal{R}(\mathrm{y},\mathrm{z})$ . The $\mathbb{P}_{D1}$ (red crosses) and $\mathbb{P}_{D2}$ (blue diamonds) approximations are shown at every third value. Exact values are on the solid line.

7.3. Residue-saddlepoint inversion of $\mathcal{P}(\mathrm{y},\mathrm{z})$

From Theorem 3,

(23) \begin{equation} \mathcal{P}(\mathrm{y},\mathrm{z}) = \sum_{m=0}^{\infty}\mathcal{P}_{m}(\mathrm{z})\mathrm{y}^{m} = \frac{1}{1-\mathrm{y}}\bigg[1 - \frac{\mathrm{y}\{1-\mathcal{F}(\mathrm{z})\}}{1-\mathcal{H}(\mathrm{y},\mathrm{z})}\bigg]. \end{equation}

A residue expansion as in Lemma 1, which inverts in $\mathrm{y}\leftrightsquigarrow m$ to approximate $\mathcal{P}_{m}(\mathrm{z})$ from the right-hand expression, is

(24) \begin{equation} \mathcal{\hat{P}}_{m}(\mathrm{z}) = \frac{1}{\hat{y}(\mathrm{z})^{m}}\frac{1-\mathcal{F}(\mathrm{z})} {\{\hat{y}(\mathrm{z})-1\}\mathcal{H}^{\prime}_{\mathrm{y}}\{\hat{y}(\mathrm{z}),\mathrm{z}\}}, \quad \mathrm{z}<\mathfrak{r}, \end{equation}

where $\hat{y}(\mathrm{z})$ is the smallest positive root of $1-\mathcal{H}(\mathrm{y},\mathrm{z})=0$ in y. This compares with the residue expansion in Section 7.1 for inverting the double GF $\mathcal{R}(\mathrm{y},\mathrm{z})$ in z to get $\mathcal{R}_{n}(\mathrm{y})$ . In that expansion, $\hat{z}(\mathrm{y})$ is the smallest positive root of $1-\mathcal{H}(\mathrm{y},\mathrm{z})=0$ in z rather than y. Thus, the functions $\hat{z}(\mathrm{y})$ and $\hat{y}(\mathrm{z})$ are inverse functions.

At $\mathrm{z}=1$ , $\mathcal{\hat{P}}_{m}(1)=1$ and $\mathcal{\hat{P}}_{m}(\mathrm{z})$ has a removable singularity, as may be shown using l’Hôpital’s rule. To show this, we compute $\hat{y}(1)=1$ and note that

\begin{align*} \frac{{\mathrm{d}}\hat{y}(\mathrm{z})}{{\mathrm{d}}\mathrm{z}} = -\frac{\mathcal{H}^{\prime}_{\mathrm{z}}\{\hat{y}(\mathrm{z}),\mathrm{z}\}} {\mathcal{H}^{\prime}_{\mathrm{y}}\{\hat{y}(\mathrm{z}),\mathrm{z}\}}<0.\end{align*}

Thus, $\hat{y}(\mathrm{z})-1\gtrless0$ if and only if $\mathrm{z}\lessgtr1$ , which ensures that $\mathcal{\hat{P}}_{m}(\mathrm{z})>0$ . The factor $1-\mathcal{F}(\mathrm{z})$ in $\mathcal{\hat{P}}_{m}(\mathrm{z})$ constrains $\mathrm{z}\in(0,\mathfrak{r})$ , where $\mathfrak{r}>1$ is the convergence radius of $\mathcal{F}(\mathrm{z})$ . The right tail for the distribution of $P_{m}$ uses a saddlepoint approaching $\ln\mathfrak{r}$ . As $\mathrm{z}\uparrow\mathfrak{r}$ , $\hat{y}(\mathrm{z})\downarrow\hat{y}(\mathfrak{r})<\hat{y}(1)=1$ .

For Example 9, the left panel of Figure 5 shows the $\mathbb{P}_{R1}$ and $\mathbb{P}_{R2}$ residue-saddlepoint approximations from inverting $\mathcal{P}(\mathrm{y},\mathrm{z})$ , which are remarkably accurate. The right panel shows the percentage relative errors with a maximum error of $-2\%$ above the median for both approximations. The two double-saddlepoint approximations $\mathbb{P}_{D1}$ and $\mathbb{P}_{D2}$ from inverting $\mathcal{P}(\mathrm{y},\mathrm{z})$ are also shown, and both have a maximum percentage relative error of $-4\%$ above the median.

Figure 5. Left: Approximations for the survival function of $P_{5}$ by inverting $\mathcal{P}(\mathrm{y},\mathrm{z})$ . The $\mathbb{P}_{R1}$ (solid red) and $\mathbb{P}_{R2}$ (dashed blue) approximations are shown with exact values (circles) for the odd integers n. Right: Percentage relative errors respecting the tails are given for the $\mathbb{P}_{R1}$ and $\mathbb{P}_{R2}$ approximations with matching line types. Also shown are the errors for the double-saddlepoint $\mathbb{P}_{D1}$ (red crosses) and $\mathbb{P}_{D2}$ (blue diamonds) approximations from $\mathcal{P}(\mathrm{y},\mathrm{z})$ .

Example 11. (Reward as an interarrival random walk.) Here, $\mathcal{H}(\mathrm{y},\mathrm{z})=\mathcal{F}\{\mathrm{z}\mathcal{G}_{0}(\mathrm{y})\}$ and we assume both $\mathcal{F}$ and $\mathcal{G}_{0}$ convergence in open circles of radius $\mathfrak{r}$ and $\mathfrak{r}_{0}$ . A residue approximation for the first-passage PGF $\mathcal{P}_{m}(\mathrm{z})$ of $P_{m}$ , as given in (24), uses the value $\hat{y}(\mathrm{z})$ which solves $\mathcal{G}_{0}(\mathrm{y})=1/\mathrm{z}$ . The expression in (24) reduces to

\begin{align*} \mathcal{\hat{P}}_{m}(\mathrm{z}) = \frac{1}{\hat{y}(\mathrm{z})^{m}}\frac{1-\mathcal{F}(\mathrm{z})} {\{\hat{y}(\mathrm{z})-1\}\mu\mathrm{z}\mathcal{G}^{\prime}_{0}\{\hat{y}(\mathrm{z})\}}, \quad 0<\mathrm{z}<\mathfrak{r}. \end{align*}

As $\mathrm{z}\downarrow0$ , $\mathcal{\hat{P}}_{m}(\mathrm{z})\rightarrow0=\mathcal{P}_{m}(0)$ for $m\geq1$ , as shown in the Supplementary Material.

For Example 1, $P_{m}$ has a negative Binomial(m, p) distribution. The pole is $\hat{y}(\mathrm{z})=(1-q\mathrm{z})/(p\mathrm{z})$ and the residue approximation is exact:

\begin{align*} \mathcal{\hat{P}}_{m}(\mathrm{z}) = \bigg(\frac{p\mathrm{z}}{1-q\mathrm{z}}\bigg)^{m} = \mathcal{P}_{m}(\mathrm{z}).\end{align*}

7.4. Double saddlepoint inversion of $\mathcal{P}(\mathrm{y},\mathrm{z})$

This approach for inverting $\mathcal{P}(\mathrm{y},\mathrm{z})$ in (23) to determine $\mathbb{P}\{P_{m}\geq n\}$ fixes the value of $m\leftrightsquigarrow\mathrm{y}$ in the double-saddlepoint expressions given in the Supplementary Material. The resulting $\mathbb{P}_{D1}$ (red crosses) and $\mathbb{P}_{D2}$ (blue diamonds) approximations are shown in the left panel of Figure 6 with exact values (solid line). Percentage relative errors are shown in the right panel of Figure 6 as red crosses and blue diamonds respectively.

Figure 6. Left: Approximations for the survival function of $P_{5}$ by inverting $\mathcal{P}(\mathrm{y},\mathrm{z})$ . Double-saddlepoint approximations $\mathbb{P}_{D1}$ (red crosses) and $\mathbb{P}_{D2}$ (blue diamonds) are shown at every third value. Exact values are on the solid line. Right: Percentage relative errors respecting the tails for the residue-saddlepoint approximations $\mathbb{P}_{R1}$ (solid red) and $\mathbb{P}_{R2}$ (dashed blue) which invert $\mathcal{P}(\mathrm{y},\mathrm{z})$ , the residue-saddlepoint approximations $\mathbb{P}_{R1}$ (dotted green) and $\mathbb{P}_{R2}$ (dot-dashed magenta) which invert $\mathcal{R}(\mathrm{y},\mathrm{z})$ , and the double-saddlepoint approximations $\mathbb{P}_{D1}$ (red crosses) and $\mathbb{P}_{D2}$ (blue diamonds) from inverting $\mathcal{P}(\mathrm{y},\mathrm{z})$ .

7.5. Summary of methods

Eight plots estimating the first-passage distribution of $P_{5}$ have been provided and two conclusions are apparent for this example and which we suggest will apply more generally. First, inverting $\mathcal{P}(\mathrm{y},\mathrm{z})$ rather than $\mathcal{R}(\mathrm{y},\mathrm{z})$ tends to result in slightly greater accuracy regardless of the particular method used. Secondly, the residue-saddlepoint approximations tend to provide greater accuracy than the double-saddlepoint approximations regardless of the type of continuity correction and regardless of the double GF which is inverted. The right panel of Figure 6 overlays six relative error plots for the four residue-saddlepoint approximation plots inverting both $\mathcal{P}(\mathrm{y},\mathrm{z})$ and $\mathcal{R}(\mathrm{y},\mathrm{z})$ and the two double-saddlepoint approximations from inverting $\mathcal{P}(\mathrm{y},\mathrm{z})$ . This shows that the residue-saddlepoint $\mathbb{P}_{R1}$ and $\mathbb{P}_{R2}$ plots (solid red, dashed blue) for inverting $\mathcal{P}(\mathrm{y},\mathrm{z})$ are the most accurate. The double-saddlepoint methods $\mathbb{P}_{D1}$ and $\mathbb{P}_{D2}$ still achieve remarkable accuracy even though they are not the most accurate.

8.1. Double-saddlepoint inversions

To approximate $\mathbb{P}\{R(n)\geq x\}$ for $x>0$ by inverting $\mathcal{R}_{\mathrm{c}0}(\mathrm{r},\mathrm{z})$ in (9), the double-saddlepoint Skovgaard approximation $\mathbb{P}_{\mathrm{C}}$ in (54) of the Supplementary Material is used as described in Section 3.3. The expression $\ln\mathcal{R}_{\mathrm{c}0}(\mathrm{r},{\mathrm{e}}^{\mathrm{s}})$ is used as the CGF for the conditional joint distribution of $\{R(T),T\}$ given $R(T)>0$ . Inversion in $(\mathrm{r},\mathrm{s})\leftrightsquigarrow(x,n)$ gives

\begin{align*} \mathbb{P}\{R(n)\geq x\} \approx c_{n}\mathbb{P}_{\mathrm{C}}\{R(T)\geq x\mid T=n,\ R(T)>0\}, \quad x>0,\end{align*}

where $c_{n}=\mathbb{P}\{R(T)>0\mid T>0\}=\mathbb{P}\{F\leq n\}$ .

Double-saddlepoint inversion of $\mathcal{P}_{\mathrm{c}}(\mathrm{r},\mathrm{z})$ in (25) inverts the other way round. Associate $(\mathrm{r},\mathrm{z}={\mathrm{e}}^{\mathrm{s}})\leftrightsquigarrow(x,n)$ and treat $\mathcal{K}(\mathrm{r},\mathrm{s})=\ln\mathcal{P}_{\mathrm{c}}({-}\mathrm{r},{\mathrm{e}}^{\mathrm{s}})$ as the joint CGF of $\{T,P_{T}\}$ . Here, $P_{T}\mid T=x$ has the discrete distribution of $P_{x}$ and marginally T has improper Lebesgue measure on $(0,\infty)$ . This joint distributional structure is reflected in the joint MGF

(29) \begin{equation} \mathcal{P}_{\mathrm{c}}({-}\mathrm{r},{\mathrm{e}}^{\mathrm{s}}) = \int_{0}^{\infty}\mathbb{E}\{{\mathrm{e}}^{\mathrm{s}P_{x}}\}{\mathrm{e}}^{\mathrm{r}x}\,{\mathrm{d}} x. \end{equation}

Then, $\mathbb{P}\{P_{x}\geq n\}=\mathbb{P}\{P_{T}\geq n\mid T=x\}$ , with $(x,n)\leftrightsquigarrow(\mathrm{r},\mathrm{s})$ . This uses either of the two continuity-corrected Skovgaard approximations $\mathbb{P}_{D1}$ or $\mathbb{P}_{D2}$ in the Supplementary Material.

8.2. Residue-saddlepoint inversions

To compute the distribution of R(n), a residue approximation first inverts $\mathcal{R}_{\mathrm{c}0}(\mathrm{r},\mathrm{z})$ in (9) in the variable $\mathrm{z}\leftrightsquigarrow n$ as described in Section 3.3. This gives $\mathcal{\hat{R}}_{\mathrm{c}n}(\mathrm{r})$ in (10) as an approximation to the MGF $\mathcal{R}_{n}(\mathrm{r})=\mathbb{E}\{{\mathrm{e}}^{\mathrm{r}R(n)})$ . Now, $\mathcal{\hat{R}}_{\mathrm{c}n}(\mathrm{r})$ is used as a surrogate for $\mathcal{R}_{n}(\mathrm{r})$ in a single-saddlepoint approximation of the Lugannani–Rice [Reference Lugannani and Rice15] type for a continuous distribution to compute an approximation to $\mathbb{P}\{R(n)>x\}$ . See Example 13.

When inverting $\mathcal{P}_{\mathrm{c}}(\mathrm{r},\mathrm{z})$ in (25) instead of $\mathcal{R}_{\mathrm{c}0}(\mathrm{r},\mathrm{z})$ , the residue-saddlepoint approximation first inverts $\mathcal{P}_{\mathrm{c}}(\mathrm{r},\mathrm{z})$ in the other continuous variable $\mathrm{r}\leftrightsquigarrow x$ and not in the discrete variable $\mathrm{z}\leftrightsquigarrow n$ . From (29), the PGF for $P_{x}$ is

\begin{align*} \mathcal{P}_{x}(\mathrm{z}) = [{\mathrm{e}}^{-x\mathrm{r}}] \int_{0}^{\infty}\mathbb{E}\{\mathrm{z}^{P_{x}}\}{\mathrm{e}}^{-\mathrm{r}x}\,{\mathrm{d}} x = [{\mathrm{e}}^{-x\mathrm{r}}]\mathcal{P}_{\mathrm{c}}(\mathrm{r},\mathrm{z}).\end{align*}

Residue expansions which invert LTs and MGFs in this continuous context were developed in [Reference Butler6] and involve more conditions than those for inverting GFs as in Lemma 1. The derivation of the next theorem is given in the Supplementary Material. Define $\hat{r}(\mathrm{z})$ to be the largest real root of $0=1-\mathcal{H}_{\mathrm{c}}({-}\mathrm{r},\mathrm{z})$ in r. Note that $-\hat{r}(\mathrm{z})$ and $\hat{z}(\mathrm{r})$ are both roots of the same equation and are, in fact, inverse functions. Implicit differentiation shows that $\hat{r}(\mathrm{z})$ is a strictly increasing function of z. Since $\hat{r}(1)=0$ , then $\hat{r}(\mathrm{z})\gtrless0$ for $\mathrm{z}\gtrless1$ .

Theorem 5. (Residue approximation for $\mathcal{P}_{x}(\mathrm{z})$ .) Suppose $\mathcal{H}_{\mathrm{c}}(\mathrm{r},\mathrm{z})=\mathbb{E}\{{\mathrm{e}}^{\mathrm{r}R}\mathrm{z}^{F}\}$ has a maximal convergence region which is an open neighbourhood of $(\mathrm{r},\mathrm{z})=(0,1)$ . Let $\mathcal{F}(\mathrm{z})$ have convergence radius $\mathfrak{r}\leq\infty$ , and suppose condition $\mathcal{Z}$ below holds. Then, a residue approximation for $\mathcal{P}_{x}(\mathrm{z})$ is

(30) \begin{equation} \mathcal{\hat{P}}_{x}(\mathrm{z}) = {\mathrm{e}}^{\hat{r}(\mathrm{z})x}\frac{1-\mathcal{F}(\mathrm{z})} {-\hat{r}(\mathrm{z})\mathcal{H}^{\prime}_{\mathrm{c}\mathrm{r}}\{-\hat{r}(\mathrm{z}),\mathrm{z}\}}, \quad \mathrm{z}<\mathfrak{r}. \end{equation}

$(\mathcal{Z})$ For each $\mathrm{z}<\mathfrak{r}$ , there exist $\eta>0$ and $\varepsilon>0$ such that $|\mathcal{H}^{\prime}_{\mathrm{c}\mathrm{r}}\{-\hat{r}(\mathrm{z})-\eta+\mathrm{i}w, \mathrm{z}\}| = o(|w|^{-\varepsilon})$ as $|w|\rightarrow\infty$ .

Note that $-\{1-\mathcal{F}(\mathrm{z})\}/\hat{r}(\mathrm{z})>0$ for all z and has a removable singularity at $\mathrm{z}=1$ ; thus $\mathcal{\hat{P}}_{x}(\mathrm{z})>0$ for all z. As $\mathrm{z}\rightarrow1$ , then $\hat{r}^{\prime}(\mathrm{z})\rightarrow\mu/\rho$ and, by l’Hôpital’s rule, $\mathcal{\hat{P}}_{x}(\mathrm{z})\rightarrow1$ . The condition $\mathcal{Z}$ placed on $\mathcal{H}_{\mathrm{c}}(\mathrm{r},\mathrm{z})$ is extremely weak and holds for all the practical joint distributions found in [Reference Johnson, Kotz and Balakrishnan9].

Example 12. (Reward as an interarrival random walk.) Suppose each tick of the clock returns a continuous reward that has MGF $\mathcal{G}_{0}(\mathrm{r})$ and which is awarded at the end of the interarrival. The mixed MGF-PGF for (R, F) is $\mathcal{H}_{\mathrm{c}}(\mathrm{r},\mathrm{z})=\mathcal{F}\{\mathrm{z}\mathcal{G}_{0}(\mathrm{r})\}$ , and is defined over $\{(\mathrm{r},\mathrm{z})\in\mathbb{R}\times[0,\infty)\,:\, \mathrm{z}\mathcal{G}_{0}(\mathrm{r})<\mathfrak{r},\ \mathrm{r}<b_{0}\}$ , where $b_{0}>0$ is the convergence bound of $\mathcal{G}_{0}$ .

The residue approximation for $\mathcal{R}_{n}(\mathrm{r})$ , the MGF of reward at time n, is

(31) \begin{equation} \mathcal{\hat{R}}_{\mathrm{c}n}(\mathrm{r}) = \mathcal{G}_{0}(\mathrm{r})^{n}\frac{\mathcal{F}_{S}\{1/\mathcal{G}_{0}(\mathrm{r})\}}{\mu}, \quad \mathcal{G}_{0}^{-1}(1/\mathfrak{r})<\mathrm{r}<b_{0}. \end{equation}

This follows from the fact that $\hat{z}_{\mathrm{c}}(\mathrm{r})=1/\mathcal{G}_{0}(\mathrm{r})$ and $\mathcal{H}_{\mathrm{c}\mathrm{z}}^{\prime}\{\mathrm{r},\hat{z}_{\mathrm{c}}(\mathrm{r})\} = \mathcal{F}^{\prime}(1)\mathcal{G}_{0}(\mathrm{r})$ .

The residue approximation for $\mathcal{P}_{x}(\mathrm{z})$ , the PGF of the first-passage time to reward x, is

\begin{equation*} \mathcal{\hat{P}}_{x}(\mathrm{z}) = {\mathrm{e}}^{\hat{r}(\mathrm{z})x}\frac{1-\mathcal{F}(\mathrm{z})} {-\hat{r}(\mathrm{z})\mu\mathrm{z}\mathcal{G}^{\prime}_{0}\{-\hat{r}(\mathrm{z})\}}, \quad \mathrm{z}<\mathfrak{r}, \end{equation*}

where $-\hat{r}(\mathrm{z})=\mathcal{G}_{0}^{-1}(1/\mathrm{z})$ and $\mathcal{H}^{\prime}_{\mathrm{c}\mathrm{r}}\{-\hat{r}(\mathrm{z}),\mathrm{z}\} = \mathcal{F}^{\prime}(1)\mathrm{z}\mathcal{G}^{\prime}_{0}\{-\hat{r}(\mathrm{z})\}$ has been used.

Example 13. (Reward as an interarrival random walk.) Suppose each tick of the clock returns an inverse Gaussian(2, 4) reward with mean 2, variance 2, and MGF $\mathcal{G}_{0}(\mathrm{r})=$ $\exp(2-2\sqrt{1-2\mathrm{r}})$ . Then, R given $F=n$ is inverse Gaussian $(2n,4n^{2})$ with MGF $\mathcal{G}_{0}(\mathrm{r})^{n}$ . Taking

\begin{align*}\mathcal{F}(\mathrm{z})=\frac{\mathrm{z}/3}{(1-2\mathrm{z}/3)}\end{align*}

with $\mu=3$ and $\mathbb{V}\{F\}=6$ , then

(32) \begin{equation} \mathcal{H}_{\mathrm{c}}(\mathrm{r},\mathrm{z}) = \frac{(\mathrm{z}/3){\mathrm{e}}^{2-2\sqrt{1-2\mathrm{r}}}}{1-(2\mathrm{z}/3){\mathrm{e}}^{2-2\sqrt{1-2\mathrm{r}}}}, \end{equation}

and $\rho=\mathbb{E}\{R\}=6$ , $\mathbb{V}\{R\}=30$ . The residue approximation for the MGF of R(n), the reward at time n, has $\hat{z}(\mathrm{r})=1/\mathcal{G}_{0}(\mathrm{r})$ so that (31) is

\begin{align*} \mathcal{\hat{R}}_{\mathrm{c}n}(\mathrm{r}) = \mathcal{G}_{0}(\mathrm{r})^{n+1}\frac{1}{3\mathcal{G}_{0}(\mathrm{r})-2}, \quad-0.223<\mathrm{r}<\frac12. \end{align*}

For $n=10$ we are able to compute the exact MGF of R(10) as $\mathcal{R}_{10}(\mathrm{r})=[\mathrm{z}^{10}]\mathcal{R}_{\mathrm{c}}(\mathrm{r},\mathrm{z})$ , which is obtained by symbolically computing $\big(\partial^{10}\mathcal{R}_{\mathrm{c}}(\mathrm{r},\mathrm{z})/\partial\mathrm{z}^{10}\big|_{\mathrm{z}=0}\big)/10!$ . From this, we compute the exact mean and variance as $16.07$ and $37.15$ . We remove the point mass of $\big(\frac23\big)^{-10}$ at 0 and work with the conditional MGF of R(10) given $R(10)>0$ of the form $\mathcal{R}_{10+}(\mathrm{r}) = \big\{\mathcal{R}_{10}(\mathrm{r})-\big(\frac23\big)^{-10}\big\}/\big\{1-\big(\frac23\big)^{-10}\big\}$ .

Figure 7 compares three strategies for approximating the survival function of R(10). The solid line is $\big\{1-\big(\frac23\big)^{-10}\big\}$ times the single-saddlepoint approximation applied using the CGF $\ln\mathcal{R}_{10+}(\mathrm{r})$ . The red dashed line is the single-saddlepoint approximation applied to $\ln\mathcal{\hat{R}}_{\mathrm{c},10}(\mathrm{r})$ . The dotted line is a Normal(16,20) approximation based on the moment approximations of Corollary 2. The normal approximation performs badly since the variance expansion value of 20 is quite far from the true variance, $37.15$ . The right panel plots the percentage relative error of the saddlepoint approximation using $\mathcal{\hat{R}}_{\mathrm{c},10}(\mathrm{r})$ and taking the saddlepoint approximation based on $\mathcal{R}_{10+}(\mathrm{r})$ as the standard. What this plot shows is that little is lost in applying a saddlepoint approximation to the residue approximation as opposed to the exact MGF $\mathcal{R}_{10+}(\mathrm{r})$ .

Figure 7. Left: Approximations for the survival function of R(10) using single-saddlepoint inversion of $\mathcal{R}_{10+}(\mathrm{r})$ (solid black line), single-saddlepoint inversion of $\mathcal{\hat{R}}_{\mathrm{c},10}(\mathrm{r})$ (dashed red line), and a Normal(16, 20) approximation (dotted line). Right: Plot of percentage relative errors respecting the tails for inversion using $\mathcal{\hat{R}}_{\mathrm{c},10}(\mathrm{r})$ as compared with inversion using $\mathcal{R}_{10+}(\mathrm{r})$ .

Example 14. (Reward as an interarrival random walk.) In the context of Example 13, we assess the approximation to the PGF for the first-passage time $P_{x}$ given in Theorem 5. With a threshold of $x=15$ , the moments of $P_{x}$ are approximated as in Corollary 9 as $\mathbb{E}\{P_{15}\}=10.25$ and $\mathbb{V}\{P_{15}\}=3.75$ . The dominant pole for $\mathcal{\hat{P}}_{x}(\mathrm{z})$ in (30) is the real solution in r to

(33) \begin{equation} 1 = \mathrm{z}\mathcal{G}_{0}(\mathrm{r}) = \mathrm{z}\exp(2-2\sqrt{(1+2\mathrm{r})}) \end{equation}

or $\hat{r}(\mathrm{z})=\frac{1}{2}\ln\mathrm{z}+\frac{1}{8}(\!\ln\mathrm{z})^{2}$ . Since the exponential is a multi-function, there are also pairs of complex-conjugate solutions to (33) given as

(34) \begin{equation} \hat{r}_{\pm k}(\mathrm{z}) = \bigg\{\frac{1}{2}\ln\mathrm{z} + \frac{1}{8}(\ln\mathrm{z})^{2} - \frac{\pi^{2}k^{2}}{2}\bigg\} \pm \mathrm{i}\pi k\bigg(1+\frac{\ln\mathrm{z}}{2}\bigg), \quad k\geq1. \end{equation}

These roots form additional simple poles which are ignored when using the approximation $\mathcal{\hat{P}}_{x}(\mathrm{z})$ but could be added to the expansion for greater accuracy, as described in [Reference Borovkov and Rogozin2, Section 2.3]. The dominant pole leads to the approximation

\begin{align*} \mathcal{\hat{P}}_{x}({\mathrm{e}}^{\mathrm{r}}) = {\mathrm{e}}^{x(\mathrm{r}/2+\mathrm{r}^{2}/8)}\frac{({\mathrm{e}}^{\mathrm{r}}-1)\big(1+\frac{1}{2}\mathrm{r}\big)} {(1-2/3{\mathrm{e}}^{\mathrm{r}})(3\mathrm{r})\big(1+\frac{1}{4}\mathrm{r}\big)}, \quad -2<\mathrm{r}<\ln\tfrac{2}{3}=0.405. \end{align*}

The value $\mathcal{\hat{P}}_{x}(1)=1$ holds by l’Hôpital’s rule. Moment computations from this expression give an approximate mean of $10.25$ and variance of $9.646$ , with the variance quite different from the expansion value of $3.75$ .

Figure 8 shows smooth plots for the two continuity-corrected single-saddlepoint survival approximations depicted in solid red and dashed blue which are visually coincident. A continuity-corrected Normal $(10.25,3.75)$ is shown as a dotted line.

For $\mathcal{\hat{P}}_{x}({\mathrm{e}}^{\mathrm{r}})$ to be a valid approximation to $\mathcal{P}_{x}({\mathrm{e}}^{\mathrm{r}})$ , the conditions of Theorem 5 must hold for this example. Condition $\mathcal{Z}$ is shown to hold in the Supplementary Material.

To assess the accuracy of the approximations in Figure 8, we consider the exact MGF $\mathcal{P}_{x}({\mathrm{e}}^{\mathrm{r}})$ . When the additional residue expansion terms from all the poles given in (34) are added up to approximate $\mathcal{P}_{x}({\mathrm{e}}^{\mathrm{r}})$ , then we show in the Supplementary Material using a non-trivial argument that this infinite sum converges pointwise to the exact MGF $\mathcal{P}_{x}({\mathrm{e}}^{\mathrm{r}})$ . This development is based on [Reference Borovkov and Rogozin2, Corollary 2.2], which specifies when such infinite residue expansions are pointwise convergent to the exact values they approximate.

Figure 8. Two continuity-corrected single-saddlepoint approximations (solid red and dashed blue lines) for the survival function of $P_{15}$ are compared with a continuity-corrected Normal $(10.25,3.75)$ approximation (dotted line).

The order of the various residue terms which contribute in this infinite expansion reveal that $\mathcal{\hat{P}}_{x}({\mathrm{e}}^{\mathrm{r}})$ and $\mathcal{P}_{x}({\mathrm{e}}^{\mathrm{r}})$ are computationally in agreement to 16-significant-digit accuracy or the accuracy capable when using 64-bit floating-point arithmetic. As a consequence, saddlepoint approximation using $\mathcal{\hat{P}}_{x}({\mathrm{e}}^{\mathrm{r}})$ is indistinguishable from the same computations using the exact MGF $\mathcal{P}_{x}({\mathrm{e}}^{\mathrm{r}})$ . Let $\mathcal{\hat{E}}_{k}({\mathrm{e}}^{\mathrm{r}})$ denote the additive contribution to the infinite residue expansion from the complex conjugate pair of poles $\hat{r}_{k}(\mathrm{z})$ and $\hat{r}_{-k}(\mathrm{z})$ . Figure 9 in the Supplementary Material plots $10^{32}\mathcal{\hat{E}}_{1}({\mathrm{e}}^{\mathrm{r}})/\mathcal{\hat{P}}_{x}({\mathrm{e}}^{\mathrm{r}})$ against r, and we see that over the range of saddlepoints $({-}1.9,0.36)$ used in the approximation,

\begin{align*} \max_{-1.9\leq\mathrm{r}\leq0.36}\frac{|\mathcal{\hat{E}}_{1}({\mathrm{e}}^{\mathrm{r}})|} {\mathcal{\hat{P}}_{x}({\mathrm{e}}^{\mathrm{r}})} < 5.7\times10^{-32}. \end{align*}

Likewise, the two derivatives of this ratio have similar size. This means that saddlepoint computations based on $\mathcal{\hat{P}}_{x}({\mathrm{e}}^{\mathrm{r}})$ and $\mathcal{\hat{P}}_{x}({\mathrm{e}}^{\mathrm{r}})+\mathcal{\hat{E}}_{1}({\mathrm{e}}^{\mathrm{r}})$ are computationally indistinguishable using 64-bit arithmetic. This is explained by the asymptotic orders of the various terms. For example, in the left tail where the error is the greatest at $\mathrm{r}=-1.9$ , the ratio of exponential factors in the residue terms for poles at $\hat{r}_{\pm k}$ versus the dominant pole $\mathrm{r}/2+\mathrm{r}^{2}/8$ are

\begin{align*} \exp\big[x\operatorname{Re}\{\hat{r}_{\pm1}({\mathrm{e}}^{\mathrm{r}})\} - x(\mathrm{r}/2+\mathrm{r}^{2}/8)\big]_{\mathrm{r}=-1.9} & = 7.1\times10^{-33}, \\ \exp\big[x\operatorname{Re}\{\hat{r}_{\pm2}({\mathrm{e}}^{\mathrm{r}})\} - x(\mathrm{r}/2+\mathrm{r}^{2}/8)\big]_{\mathrm{r}=-1.9} & = 2.6\times10^{-129}. \end{align*}

The contributions from $\mathcal{\hat{E}}_{3}({\mathrm{e}}^{\mathrm{r}})$ , etc. are even smaller. Thus we see that $\mathcal{\hat{P}}_{x}({\mathrm{e}}^{\mathrm{r}})$ and $\mathcal{P}_{x}({\mathrm{e}}^{\mathrm{r}})$ are indistinguishable in floating-point computation and the errors in the saddlepoint approximations of Figure 8 are only due to the saddlepoint methods themselves and not due to using an approximation for $\mathcal{P}_{x}({\mathrm{e}}^{\mathrm{r}})$ .

In Example 14, there are multiple poles involved in the residue expansion. However, the pole at $\hat{r}(\mathrm{z})$ is dominant enough that its single-term residue expansion $\mathcal{\hat{P}}_{x}(\mathrm{z})$ is computationally indistinguishable from the exact PGF $\mathcal{P}_{x}(\mathrm{z})$ . This is not always the case. Examples 15 and 16 in the Supplementary Material replace the inverse Gaussian reward distribution for $\mathcal{G}_{0}(\mathrm{r})$ with a Gamma(2,1) distribution. For the first passage in Example 16, there are two real solutions to $1=\mathrm{z}\mathcal{G}_{0}(\mathrm{r})$ in r which result in two poles. Using $\mathcal{\hat{P}}_{x}(\mathrm{z})$ based on only the dominant pole does not provide adequate approximation to the exact PGF, and in fact $\mathcal{\hat{P}}_{x}(\mathrm{z})$ is not even analytic at $\mathrm{z}=0$ as it needs to be if it is to be a PGF.

Examples 14 and 16 suggest that residue expansions such as $\mathcal{\hat{P}}_{x}(\mathrm{z})$ should not be used blindly as a surrogate for $\mathcal{P}_{x}(\mathrm{z})$ . The degree of dominance shown by the dominant pole depends on the spacing of the various poles for $\mathcal{P}_{\mathrm{c}}(\mathrm{r},\mathrm{z})$ , and without adequate dominance, additional residue terms may be required.