2.1. Expansions for f(t) using pole b

The inversion formula for f (t) is

(4) \begin{equation}f(t)=\frac{1}{2\pi \im}\lim_{N\rightarrow\infty}\int_{b^{-}-\im N}^{b^{-} +\im N}\mathcal{M}(s)\re^{-st}\sd s,\end{equation}

where b ⁻ = b − ε with 0 < ε < b, and the integration is performed along the vertical line {Re(s) = b ⁻} and within the convergence domain for $\mathcal{M}$. Consider deforming this contour to allow integration along the vertical line {Re(s) = b ⁺ = b + ε} which is outside the convergence domain of $\mathcal{M}$ but still within its analytic continuation. If singularity b is an m-pole of $\mathcal{M}$ and the only singularity in {b ≤ Re(s) ≤ b ⁺}, then we can use Cauchy’s residue theorem to relate the two integrals in terms of the residue of $\mathcal{M} (s)\re^{-st}$ at b > 0. If, apart from the pole at b, $\mathcal{M}$ possesses an analytic continuation to {Re(s) < b + ε ₀} for ε ₀ > ε, then Cauchy’s theorem can be applied to integration over the rectangular curve with the four corners b ⁻ ± iN and b ⁺ ± iN to give

(5) \begin{align} f(t)&=\frac{1}{2\pi \im}\lim_{N\rightarrow\infty}\int_{b^{-}-\im N}^{b^{-} +\im N}\mathcal{M}(s)\re^{-st}\sd s\nn\\ &=-\operatorname*{Res}\{\mathcal{M}(s)\re^{-st} ;\ b\} +\frac{1}{2\pi \im}\lim_{N\rightarrow\infty}\bigg( \int_{b^{+}-\im N}^{b^{+}+\im N}+\int_{b^{+}+\im N}^{b^{-}+\im N}+\int_{b^{-}-\im N} ^{b^{+}-\im N}\bigg) \mathcal{M}(s)\re^{-st}\sd s.\nonumber\\ \end{align}

Subject to condition $\mathcal{X}$ of Theorem 1 below, the last two integrals are negligible for large N so that (5) becomes

(6) \begin{equation}f(t)=-\operatorname*{Res}\{\mathcal{M}(s)\re^{-st};\ b\}+R_{1}(t),\end{equation}

with

(7) \begin{equation}R_{1}(t)=\re^{-b^{+}t}\frac{1}{2\pi}\int_{-\infty}^{+\infty}\mathcal{M} (b^{+}+\im y)\re^{-\im yt}\sd y,\end{equation}

where the improper integral is a principal-value integral. With the additional uniform integrability condition $\mathcal{UI}$ of Theorem 1, then the integral in (7), whose contour lies outside the convergence domain of $\mathcal{M}$ but within the analytic continuation, can be made as small as o(1) as t → ∞. Thus, R ₁(t) = o(e^−b+t) as t → ∞ so (6) becomes

(8) \begin{equation}f(t)=\re^{-bt}\sum_{k=1}^{\mathfrak{m}}t^{k-1}\frac{({-}1)^{k}\beta_{-k}} {(k-1)!}+o(\re^{-b^{+}t}), \qquad t\rightarrow\infty.\end{equation}

Here, the residue has been determined by multiplying the Taylor expansion of e^−st about s = b with the Laurent expansion of $\mathcal{M}$ about b whose principal part is denoted as $\sum_{k=1}^{\mathfrak{m}}\beta _{-k}\break (s-b)^{-k}$. Sufficient conditions for the expansion in (8) are formalised in Theorem 1.

The leading term in this expansion is e^−btt ^𝔪−1(−1)^𝔪β _−𝔪/(𝔪 − 1)! and was given in Butler (Reference Butler2017, Proposition 1) with error o(t ^𝔪−1e^−bt) as t → ∞. By contrast, the use of Cauchy’s theorem determines the expansion in (8) with smaller-order error o(e^{−b +t}), which should prove more accurate in numerical computations.

Before stating Theorem 1, it is worth emphasizing that the integral in (7) is not the Fourier inversion that leads to the tilted density f_b+ (t) = e^{b +t}f (t). For ε < 0, $\mathcal{M} (b+\varepsilon+\im y)$, as a function of y, is the characteristic function of f_{b +ε} (t) so that the inversion holds.

(9) \begin{equation} f_{b+\varepsilon}(t)=\re^{(b+\varepsilon)t}f(t)=\frac{1}{2\pi} \int_{-\infty }^{+\infty}\mathcal{M}(b+\varepsilon+\im y) \re^{-\im yt}\sd y, \qquad\varepsilon\in({-}b,0), \end{equation}

However, for ε ≥ 0, the identity in (9) does not hold because the integral is outside the convergence domain for $\mathcal{M}(s)$.

What is most notable about Theorem 1 is that its main assumptions $\mathcal{AC}, \mathcal{X}$, and $\mathcal{UI}$ all concern $\mathcal{M}$ and not f. In many stochastic modelling applications, $\mathcal{M}$ is available but f is not, so such conditions are much simpler to check than conditions on f as were used in Butler (Reference Butler2017, Proposition 1). The only assumption made about f in Theorem 1 is that it is locally of bounded variation, i.e. for each t > 0, there is a small δ_t > 0 for which f is of bounded variation within (t − δ_t, t + δ_t). This is the weakest assumption under which the inversion formula (4) holds and is satisfied, for example, over sets for which f is finitely piecewise continuous.

The uniform integrability condition $\mathcal{UI}$ can often be established by showing that M(b ⁺ + iy) is absolutely integrable in y. The next result follows directly from the Riemann–Lebesgue lemma (Feller (Reference Feller1971, p. 513)).

Example 2. (Interfering poles.) Suppose that $\mathcal{M}$ has 𝔪-poles at b and b ± yi, with respective mth order Laurent coefficients β _−𝔪 = α ₀, $\beta_{-\mathfrak{m}}^{+}= \alpha+\gamma \im,$, and $\beta_{-\mathfrak{m} }^{-}=\bar{\beta}_{-\mathfrak{m}}^{+}=\alpha-\gamma \im$, where the superscripts ± identify the associated pole b ± yi. Then expansion (53) in Corollary 6 of Butler (Reference Butler2019a) has the form

(10) \begin{equation}f(t)=\re^{-bt}\sum_{k=1}^{\mathfrak{m}}\frac{t^{k-1}(-1)^{k}c_{k} (t)}{(k-1)!}+o(\re^{-b^{+}t}),\end{equation}

where the term with the highest order has the coefficient

c_{\mathfrak{m}}(t)=\beta_{-\mathfrak{m}} +\re^{-\im yt} \beta_{-\mathfrak{m}}^{+}+\re^{\im yt}\bar{\beta}_{-\mathfrak{m}}^{+}=\alpha_{0}+2\alpha\cos(yt)+2\gamma\sin(yt).

The periodicity of c _𝔪(t) is due to the presence of the complex conjugate pair of poles. If, instead, the poles at b ± yi are of lower order 𝔪₁ < 𝔪 then c _𝔪(t) = β _−𝔪 and the periodicity does not appear in the leading O(e^−btt ^m−1) term, but rather in the lower-order O(e^−btt ^𝔪₁−1) term.

2.2. Expansions for S(t) using pole b

An asymptotic expansion for S(t) also follows by using Cauchy’s theorem in the same manner as used in Theorem 1 under slightly weaker conditions.

Theorem 2. (Survival expansions.) Suppose absolutely continuous X with support in (− ∞, ∞) has the survival function S(t) and MGF $\mathcal{M}$ with convergence boundary {s ∈ ℂ: Re(s) = b} with 0 < b < ∞. Subject to conditions $\mathcal{AC}$ and $\mathcal{X}$ of Theorem 1 as well as $\mathcal{UI}^{{S}}$ below on $\mathcal{M}$, then

(11) \begin{equation}S(t)=\sum_{k=1}^{\mathfrak{m}}S_{\text{G}(k,b)}(t)\frac{(-1)^{k}\beta_{-k} }{b^{k}}+R_{1}^{S}(t),\end{equation}

where $\sum_{k=1}^{\mathfrak{m}}\beta_{-k}(s-b)^{-k}$ is the principal part of the Laurent expansion for $\mathcal{M}(s)$ at b, and S_{G (k,b)} is the survival function of a gamma (k, b) distribution given by

(12) \begin{equation}\label{incomplete}S_{\text{G}(k,b)}(t)=\re^{-bt}\sum_{j=0}^{k-1}\frac{(bt)^{j}}{j!}.\end{equation}

with b ⁺ = b + ε, the error is

(13) \begin{equation}R_{1}^{S}(t)=\re^{-b^{+}t}\frac{1}{2\pi}\int_{-\infty}^{\infty}\frac {\mathcal{M}(b^{+}+\im y)}{b^{+}+\im y}\re^{-\im yt}\sd y=o(\re^{-b^{+}t}), \qquad t\rightarrow\infty.\end{equation}

$\mathcal{UI}_{\rm m}$. For some T > 0, the principal-value integral

\begin{equation*}\int_{-\infty}^{+\infty}\frac{\mathcal{M}(b^{+}+\im y)}{b^{+}+\im y}\re^{-\im yt}\sd y\end{equation*}

converges uniformly in t for t ≥ T.

Proof. The proof is the same as used in Theorem 1 but applied instead to the inversion formula of S(t). See Section 5.1.4 of Butler (Reference Butler2019a) for details.

If S ₁(t) and f ₁(t) are the expansions in (11) and (8), respectively, then $S_{1}(t)=\int_{t}^{\infty}f_{1}(u)\sd u$ also the error terms in (13) and (7) are related by $R_{1} ^{S}(t)=\int_{t}^{\infty}R_{1}(u)\sd u$. Thus, expansion (11) could have been derived directly from Theorem 1 through the integration of (8) and subsequently showing that the error of o(e^−b+t) is maintained after such integration. This is a legitimate argument but it also requires using the stronger condition $\mathcal{UI}$ of Theorem 1 rather than the weaker condition $\mathcal{UI}^{{\rm S}}$ of Theorem 2. Note also that no locally bounded variation assumption is needed for S(t), as there was for f (t), since S(t) has total variation 1.

Condition $\mathcal{UI}^{{\rm S}}$ can often be shown to hold by using the following result.

Lemma 2. If $\left| \mathcal{M} (b^{+}+\im y)\right| ^{p}$ is integrable for some p > 1 then $\mathcal{UI}^{{S}}$ of Theorem 2 holds. If δ > 0 and Y > 0 exist such that $\max_{b\leq x\leq b^{+}}\left| \mathcal{M}(x+\im y)\right| <c/|y|^{\delta}$ for |y| > Y and some c > 0, then conditions $\mathcal{X}$ and $\mathcal{UI}^{{S}}$ hold.

Proof. Use Hölder’s inequality so that

(14) \begin{equation} \int_{-\infty}^{+\infty}\frac{| \mathcal{M}(b^{+}+\im y)| }{| b^{+}+\im y| }\sd y\leq\bigg\{ \int_{-\infty}^{+\infty}| \mathcal{M}(b^{+}+\im y)| ^{p}\sd y\bigg\} ^{1/p}\bigg\{ \int_{-\infty }^{+\infty}\frac{1}{| b^{+}+\im y| ^{q}}\sd y\bigg\} ^{1/q}<\infty, \end{equation}

where 1/p + 1/q = 1. The absolute integrability in (14) implies $\mathcal{UI}^{{\rm S}}$.

Example 4. (Example 2 continued.) The survival function expansion given in Butler (Reference Butler2019a, Equation (56)) has the same general form as (10) but with different coefficients $\{c_{k}^{S}(t)\colon k=1,\ldots, \mathfrak{m\}}$ replacing {c_k(t)}. After using (12) to expand S _G(m,b) and S _G(m,b+iy) in (56), then the leading term of order O(e^−btt ^m−1) in S(t) is e^−btt ^m−1(−1)^mc^Sm(t)/(m − 1)!, with

\begin{align*} c_{\mathfrak{m}}^{S}(t) &=\frac{\beta_{-\mathfrak{m} }}{b}+2{\mathrm{Re}}\Big\{ \frac{\re^{-\im yt}\beta_{-\mathfrak{m}}^{+} }{b+\im y}\Big\} \\ &=\frac{\alpha_{0}}{b}+\frac{2}{b^{2}+y^{2}}\left\{ (b\alpha+y\gamma)\cos\left( yt\right) +(b\gamma-\alpha y) \sin(yt)\right\}\!. \end{align*}

A tedious computation shows that $c_{\mathfrak{m}}^{S}(t)=\int_{t}^{\infty }c_{\mathfrak{m}}(u)\sd u$, which is consistent with the fact that the survival expansion in Butler (Reference Butler2019a, Equation (56)) is the tail area for the density expansion in Butler (Reference Butler2019a, Equation (53)) when the latter expansion is valid. Note that the hazard function $f(t)/S(t)\sim c_{\mathfrak{m}}(t)/c_{\mathfrak{m}}^{S}(t)$ as t → ∞ and this ratio is a 2π/y-periodic function without a limit as discussed in Butler (Reference Butler2017, Section 3, Example 1).

2.3. Higher-order expansions using multiple poles

We now consider higher-order expansions for f (t) and S(t) related to poles of $\mathcal{M}$ that lie further into the analytic continuation {Re(s) ≥ b}. The following alternative assumption is needed.

• $\mathcal{AC}_{{\rm m}}$. Either $\mathcal{M}$ has a monotone increasing sequence of distinct real poles b = b ₁ < · · · < b_m, where b_j is an 𝔪_j-pole, or it has a monotonic increasing sequence of complex conjugate pairs of poles denoted by b = b ₁ < Re(b _2±) < · · · < Re(b_{m ±}). Apart from these poles, $\mathcal{M}$ admits an analytic continuation into {b ≤ Re(s) < b_m + ε ₀}, with ε ₀ > 0.

Nonreal poles must occur in complex conjugate pairs b_{j +} and b_{j −} = b̄_{j +} due to the property that $\mathcal{\bar{M}}(s)=\mathcal{M} (\bar{s})$ for all s. Denote the principal part of the Laurent expansion of $\mathcal{M}$ about pole b_j as $\sum_{k=1}^{\mathfrak{m}_{j}}\beta _{-k;j}(s-b_{j})^{-k}$.

Cauchy’s theorem is now used as in Theorems 1 and 2 to give an expansion in terms of the first m residues at {b_k: k ≤ m} by displacing the contour integral from {Re(s) = b ⁻} to $\{{\mathrm{Re}} (s)=b_{m}^{+}\}$, where $b_{m}^{+}=b_{m}+\varepsilon_{m}$ with 0 < ε_m < ε ₀.

Corollary 1. (Higher-order expansions.) Suppose that the density f (t) satisfies the conditions of Theorem 1 with conditions $\mathcal{AC}_{{m}}, \mathcal{X}_{m},$, and $\mathcal{UI}_{m}$ on $\mathcal{M}$ replacing $\mathcal{AC}, \mathcal{X}$, and $\mathcal{UI}$.

• $\mathcal{X}_{{\rm m}}$. For some ε_m ∈ (0, ε ₀), $\max_{b\,\leq\,x\,\leq b_{m}^{+}}\left| \mathcal{M}(x+\im N)\right| \rightarrow0$ as N → ∞.

• $\mathcal{UI}_{\rm m}$. For some T > 0, the principal-value integral $\int_{-\infty}^{+\infty}\mathcal{M}(b_{m} ^{+}+\im y)\re^{-\im yt}\sd y$ converges uniformly in t for t ≥ T.

Then,

(15) \begin{equation}f(t)=f_{m}(t)+R_{m}(t)\coloneqq\sum_{j=1}^{m}\re^{-b_{j}t} \sum_{k=1}^{\mathfrak{m}_{j} }t^{k-1}\frac{(-1)^{k} \beta_{-k;j}}{(k-1)!}+R_{m}(t),\end{equation}

where

(16) \begin{equation}R_{m}(t)=\re^{-b_{m}^{+}t}\frac{1}{2\pi}\int_{-\infty}^{+\infty} \mathcal{M} (b_{m}^{+}+\im y)\re^{-\im yt}\sd y=o(\re^{-b_{m}^{+}t}),\qquad t\rightarrow\infty.\end{equation}

Suppose that S(t) satisfies the conditions of Theorem 2 with conditions $\mathcal{AC}_{m}$ and $\mathcal{X}_{m}$ on $\mathcal{M}$ replacing $\mathcal{AC}$ and $\mathcal{X}$, and the following assumption $\mathcal{UI}_{m}^{S}$ replacing $\mathcal{UI}^{S}$.

• $\mathcal{UI}_{{\rm m}}^{{\rm S}}$. For some T > 0, the principal-value integral $\int_{-\infty}^{+\infty} \mathcal{M} (b_{m}^{+}+\im y)(b_{m}^{+}+\im y)^{-1}\re^{-\im yt}\sd y$ converges uniformly in t for t

Then,

(17) \begin{equation}S(t)=S_{m}(t)+R_{m}^{S}(t)\coloneqq\sum_{j=1}^{m}\sum_{k=1}^{\mathfrak{m}_{j} }S_{\text{G}(k,b_{j})}(t)\frac{({-}1)^{k}\beta_{-k;j}}{b_{j}^{k}} + R_{m}^{S}(t),\end{equation}

where

(18) \begin{equation}R_{m}^{S}(t)=\re^{-b_{m}^{+}t}\frac{1}{2\pi}\int_{-\infty}^{+\infty} \frac{\mathcal{M}(b_{m}^{+}+\im y)}{b_{m}^{+}+\im y}\re^{-\im yt}\sd y = o(\re^{-b_{m}^{+} t}), \qquad t\rightarrow\infty,\end{equation}

and S _G(k,bj) denotes the survival function of a gamma (k, b_j) distribution as in (12).

If b_j is rather a complex conjugate pair b_{j +} and b_{j −}, then (15) and (17) continue to hold but with the jth term summed over both j + and j −.

2.4. Expansions using essential singularities

If b and other finite singularities are isolated essential then similar residue expansions hold under the identical conditions used in Theorems 1 and 2 and Corollary 1. The formal statements of the corresponding results are identical apart from pole orders {m_j} taking value ∞ so that residue sums are infinite. We provide two examples.

Example 6. (Noncentral χ ²(2m, λ).) The MGF is

(19) \begin{equation}\mathcal{M}(s)=\frac{1}{(1-2s)^{m}}\exp\Big\{ \frac{\lambda} {2}\Big( \frac{1}{1-2s}-1\Big) \Big\}\end{equation}

and $b=\tfrac12$ is an essential singularity. From the form in (19) we see that this distribution is the convolution of a central χ ²(2m), which contributes the factor (1 − 2s)^−m, and a compound Poisson with Poisson rate λ/2 and Y ₁ ~ Exponential $(\tfrac12)$.

For m ≥ 2, the $\mathcal{UI}$ condition holds due to Lemma 1; when m = 1, $\mathcal{UI}$ holds by an integration-by-parts argument as given in Butler (Reference Butler2019a, Section 5.1.6). The expansion errors in (7) and (13) are R ₁(t) ≡ 0 ≡ R _1S(t) so the residue expansions are exact and take the form of infinite Poisson mixtures of well-known central χ ² distributions of the form

(20) \begin{align} f(t) & = -\operatorname*{Res}\Big\{\mathcal{M}(s)\re^{-st};\ \frac12 \Big\}=\sum _{k=0}^{\infty}\frac{\re^{-\lambda/2}}{k!}\Big( \frac{\lambda}{2}\Big) ^{k}f_{\chi^{2}(2m+2k)}(t),\label{esssing}\end{align}

(21) \begin{align}S(t) & = -\operatorname*{Res}\Big\{ \frac{\mathcal{M}(s)} {s}\re^{-st};\ \frac{1}{2}\Big\} =\sum_{k=0}^{\infty}\frac{\re^{-\lambda/2}} {k!}\Big( \frac{\lambda}{2}\Big) ^{k}S_{\chi^{2}(2m+2k)}(t). \label{esssing2} \end{align}

Here f_{χ ²(2m+2k)} and S_{χ ²(2m+2k)} are the density and survival functions of a central χ ²(2m + 2k) distribution. For details, see Butler (Reference Butler2019a, Section 5.1.6).

Alternatively, these expansions can be shown to hold for any real m > 0 by taking a Taylor expansion of the exponential term in $\mathcal{M}$ about 0 and inverting the infinite sum term-by-term. Since the sum is infinite, justification for this must be based on Doetsch (Reference Doetsch1974, Theorem 30.1) which allows the inversion integral to be interchanged with the infinite sum. From this, we conclude that the identities in (20) and (21) hold for any real m > 0 and for a.e. t > 0.

2.5. Saddlepoint approximations in the analytic continuation

The expansion errors in (16) and (18) can be approximated by using either saddlepoint approximations or through numerical integration. In both integrals, if the integrand factor $\mathcal{M}(s)\re^{-st}$ has another real pole at b_{m +1} > b_m then a real saddlepoint ŝ_m ∈ (b_m, b_{m +1}) may exist to approximate the integrals. In such cases, saddlepoint approximations for R_m(t) and $R_{m}^{S}(t)$ are

(22) \begin{equation}R_{m}(t) \simeq\hat{R}_{m}(t)=\frac{1}{\sqrt{2\pi\mathcal{K}^{\prime \prime}(\hat{s}_{m})}\,}\exp\{\mathcal{K}(\hat{s}_{m})-\hat{s}_{m} t\},\end{equation}

(23) \begin{equation}R_{m}^{S}(t) \simeq\hat{R}_{m}^{S}(t)=\frac {1}{\sqrt{2\pi\{\mathcal{K}^{\prime\prime}(\mathfrak{\hat{s}}_{m} )+\mathfrak{\hat{s}}_{m}^{-2}\}}\,}\frac{1}{\mathfrak{\hat{s}}_{m}}\exp\{\mathcal{K}(\mathfrak{\hat{s}}_{m})-\mathfrak{\hat{s}}_{m}t\},\end{equation}

where $\mathcal{K}(s)=\ln\mathcal{M}(s)$ and ŝ_m and $\mathfrak{\hat{s}}_{m}$ are the respective saddlepoint solutions to and in (b_m, b_{m +1}).

An example of the setting for Proposition 1 is given in Figure 1 of Section 3.2.1. In such cases, $\mathcal{M}(s)\re^{-st}$ achieves a minimum (maximum) value in (b_m, b_{m +1}) when $\mathcal{M}(s) 0 (<0)$. Also, in such cases, $\mathcal{K}(s)=\ln\mathcal{M}(s)$ is analytic on s ∈ (b_m, b_{m +1}) × (− δ, δ) ⊂ ℂ for sufficiently small δ > 0, so that $\mathcal{K}^{\prime}(s)-t$ must have a real zero for at least one ŝ_m ∈ (b_m, b_{m +1}). Such a point ŝ_m must therefore be a saddlepoint by the Cauchy–Riemann equations. This argument does not apply when $0\in\{\mathcal{M}(s)\colon s\in(b_{m},b_{m+1})\}$ This is because $\mathcal{K}(s)$ cannot be made analytic on (b_m, b_{m +1}) × (− δ, δ) since 0 is not in the domain of the ln function and an analytic branch of the multifunction $\ln\mathcal{M}(s)$ cannot be defined on (b_m, b_{m +1}) × (− δ, δ).

With this setup, the saddlepoint formulae in (22) and (23) can be derived as in Copson (Reference Copson1965, pp. 92–93) or Murray (Reference Murray1984, Chapter 3) based upon the steepest descents method. Both (22) and (23) presume that the saddlepoints are simple meaning that $\mathcal{K}^{\prime\prime}(\hat{s}_{m})\neq0$ so the expressions are well defined as complex computations.

In many settings there is a unique real saddlepoint for $\mathcal{M} (s)\re^{-st}$ in (b_m, b_{m +1}). This is the case for some simpler distributions such as the minus log-beta distribution (Examples 1 and 3), Wilks’ likelihood ratio statistic (Examples 7 and 8), the extreme value, logistic, and hyperbolic secant distributions (Example 10). This is also the case in a relatively complicated Cramér–Lundberg example (Example 19).

In settings where b_{m +1} is not real and is rather a complex conjugate pair of poles b_{m +1} and b̄_{m +1}, then there may be multiple real or complex conjugate pairs of saddlepoints in {b_m < Re(s) < Re(b_{m +1})}. This is illustrated in both the right and left tails of the compound distribution of Example 17 (see also Butler (Reference Butler2019a, Sections 5.3 and 5.3.1)) and in a Cramér–Lundberg example (Example 18). With multiple real or multiple complex conjugate pairs of saddlepoints, then a choice must be made between them when using (22) and (23) to approximate the inversion integrals for R_m(t) and $R_{m}^{S}(t)$. We recommend the closest (pair) in terms of its real component. If one must deal with a complex conjugate pair of saddlepoints then such saddlepoints typically lie somewhere in the vicinity between b_m and b_{m +1}, and b_m and b̄_{m +1}, but not necessarily on straight lines as would happen when b_m and b_{m +1} are both real poles. Also with a complex conjugate pair of saddlepoints, both saddlepoints are used in the approximation and the contour of the inversion integral must be deformable to pass along the path of steepest descent from each saddlepoint in the pair. Thus, in this case, (22) and (23) are applied by adding the contributions from both saddlepoints.

The saddlepoint approximation in (23) will be accurate without the need to apply the method of Bleistein (Reference Bleistein1966) to deal with the integrand factor $(b_{m}^{+}+\im y)^{-1}$, since $b_{m}^{+}$ is not in the convergence domain of $\mathcal{M}$ and therefore well away from 0. If a value for $b_{0}^{+} \in(0,b)$ in (18) were chosen so it is in the convergence domain, then $R_{0}^{S}(t)=S(t)$ and the approximation in (23) using a saddlepoint in (0, b) would not be as accurate as the Lugannani and Rice (Reference Lugannani and Rice1980) saddlepoint approximation for S(t), which is based upon the Bleistein method.

If $\mathcal{M}$ is analytic on {Re(s) > b_m} then deforming the contour integral further to the right cannot change the errors R_m(t) and $R_{m}^{S}(t)$. Also, there may not be any saddlepoints in this region to approximate these errors. An example is the random variable −ExV * Exp (1), where ExV is an extreme value distribution with MGF Γ(1 − s). The convolution distribution has MGF Γ(1 + s)/(1 − s) with poles {…, −2, −1, 1} and f ₁(t) = e^−t = S ₁(t) have errors R ₁(t) = −e^−et (1 + e^−t) and $R_{1}^{S}(t)=-\re^{-\re^{t} -t}$. In examples for which $\mathcal{M}$ is meromorphic with only a finite number of poles outside of {a < Re(s) < b}, then $\mathcal{M}$ is rational and such errors are 0. For example, with poles only at a_n < · · · < a ₁ < 0 < b < · · ·< b_m, then f (t) = f_m(t) and S(t) = S_m(t) for t > 0.

2.6. Examples

The expansions and saddlepoint approximations together can provide more accurate approximations to f (t) and S(t) than the individual methods used separately as seen in the next example.

Example 7. (Wilks’ likelihood ratio statistic, $\mathcal{M}$ rational.) The likelihood ratio statistic Λ_k,p,n in MANOVA has null distribution, with a_i = n − i+ 1 and all beta variates independent. Here, k denotes dimension, p is the degrees of freedom for hypothesis, and n is the degrees of freedom for error. Taking k = 3, p = 6, and n = 7, then the MGF for − ln Λ_3,6,7 is

(24) \begin{align} \mathcal{M}(s) & = c\,\prod_{\im=1}^{3}\frac{\Gamma\{(8-\im)/2-s\} }{\Gamma\{(14-\im)/2-s\} }\nonumber\\ & = c\big\{ \big( \tfrac52-s\big) \big( 3-s\big) \big(\tfrac72-s\big) ^{2}\big( 4-s\big) \big( \tfrac92-s\big)^{2}\big( 5-s\big) \big(\tfrac{11}{2}-s\big) \big\} ^{-1} \end{align}

where $c=\tfrac{3\,274\,425}{16}$ and $s<a_{3}/2=\tfrac52=b$. The MGF has seven poles ranging from $\tfrac52$ to $\tfrac{11}{2}$ in increments of $\tfrac12$ with 1-poles at $\tfrac52$, 3, 4, 5, and $\tfrac{11}{2}$ and 2-poles at $\tfrac72$ and $\tfrac92$. This rational form for $\mathcal{M}$ occurs whenever either p is even or k is even with the latter setting following from the fact that Λ_k,p,n and Λ_p,k,p+n−k have the same distribution. From the rational form in (24), the exact density/survival functions can be easily computed through partial fraction inversion. This approach can be shown to be equivalent to that used by Schatzoff (Reference Schatzoff1966) who worked in the time domain by successively convoluting log-beta densities.

Table 1:

Various density and survival approximations, f_m(3) + R̂_m(3) and $S_{m}(3)+\hat{R}_{m}^{S}(3)$, for f (3) = 0.297 316 and S(3) = 0.206 673 47 where m denotes the number of residues included. The boldface 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. Also shown are saddlepoint approximations R̂_m(3) and $\hat{R}_{m}^{S}(3)$ for errors R_m(3) and $R_{m}^{S}(3)$, respectively, along with exact (numerically integrated) computations for these errors.

^a Indicates that R̂_m or $\hat{R}_{m}^{S}$ is undefined.

^b Indicates that the seven residue sum is exact.

In Table 1 we show the approximations of Corollary 1 for f (3) = 0.297 316 and S(3) = 0.206 673. The entries demonstrate that a sufficient number m of residue terms are needed in order for the asymptotic expansions to achieve an accurate approximation. This should have been expected in this example from the nature of the asymptotic expansion since the poles are spaced relatively close together with increments of $\tfrac12$. The entries for saddlepoint error approximations R̂_m(3) and $\hat{R}_{m}^{S}(3)$ provide guidance about the expansion error and appear quite useful in determining the value of m needed to ensure sufficient accuracy. Both approximations need to sum residues to at least the 2-pole at m = 5 in order to achieve any kind of accuracy. Inaccuracy for m ≤ 4 occurs because successively added residue terms fluctuate in sign and magnitude. This is caused by quite large residues that are not adequately dampened by the gradually diminishing exponential values {exp (− b_j t)} using poles {b_j} in increments of $\tfrac12$ which are not well spaced.

Saddlepoint density and Lugannani and Rice (Reference Lugannani and Rice1980) survival approximations applied in the convergence strip $\{s<\tfrac52\}$ of the MGF lead to the values 0.2984 and 0.2061, respectively, where the boldface digits indicate the last ‘accurate’ digit or the last digit in agreement with the exact result when both are rounded to the same number of digits. These approximations attain remarkable accuracy which is slightly better than expansions to m = 5 and exceeded only by using m = 6 residues. When term m = 7 is included, the expansions are exact.

Example 8. (Wilks’ likelihood ratio statistic, $\mathcal{M}$ is not rational.) Consider the case in which both p and k are odd. This is the setting in which there are no exact computational procedures (however, see Section 2.7 below). We take k = 3, p = 5, and n = 7 so the MGF for − ln Λ_3,5,7 is

(25) \begin{equation}\mathcal{M}(s)=c_{1}\,\frac{1}{\left( 5/2-s\right) \left( 3-s\right)\left( 7/2-s\right) ^{2}\left( 4-s\right) \left( 9/2-s\right) ^{2}}\times\frac{\Gamma(11/2-s)}{\Gamma(6-s)},\end{equation}

where $c_{1}=30240/\sqrt{\pi}$ There are poles at $\tfrac52$, 3, $\tfrac72$, 4, $\tfrac92, \tfrac{11}{2}$ with saddlepoints in between; however, above $\tfrac{11}{2}$ there are no saddlepoints. At $\tfrac{11}{2}$ and above there are simple poles at $\tfrac{11}{2},\tfrac{13}{2},\ldots$, but $\mathcal{M}(s)$ is monotone decreasing in between each pair of poles with $\infty=\mathcal{M}\{(k-\tfrac12)^{+}\}> \mathcal{M}\{(k+\tfrac12)^{-}\}=-\infty$ for k ≥ 6, so Proposition 1 does not hold. Since $\Gamma(\tfrac{11}{2}-s)/\Gamma (6-s)\sim1/\sqrt{-s}$ when Re(s) is fixed and Im (s) → ∞, then $|\mathcal{M}(s)|= O(|s|^{-{15}/{2}} )$; thus, conditions $\mathcal{AC}_{\rm m}, \mathcal{X}_{\rm m}$, and $\mathcal{UI}_{\rm m}$ of Corollary 1 hold for all m ≥ 1.

Table 2:

Various density and survival approximations for f (2.7) = 0.301 133 and S(2.7) = 0.201 400 8 where m denotes the number of residues included. See Table 1 for further description.

In Table 2 we show the accuracy that can be achieved when using the expansions and saddlepoint approximations to determine f (2.7) = 0.301 133 and S(2.7) = 0.201 400 8. Five terms (m = 5) are needed to stabilise the expansion but five terms result in 5 digit accuracy (after rounding) for both density and survival approximations. For comparison, the saddlepoint density and Lugannani and Rice (Reference Lugannani and Rice1980) survival approximations are 0.3028 and 0.2009, respectively, with 2 and 3 digit accuracy (after rounding); this is considerably better than the four term expansion but much less accurate than using five residue terms.

In both examples, the saddlepoint error approximations R̂_m and $\hat{R}_{m}^{S}$ make it very clear which values of m provide adequate and inadequate expansions. This is an important aspect of the methods if they are to ensure a specified accuracy particularly when poles are not well spaced as occurs in these two examples.

2.7. Expansions using an infinite number of poles

If $\mathcal{M}$ admits an infinite sequence of increasing poles b = b ₁ < b ₂ < · · · and Corollary 1 can be assumed to hold for all m = 1, 2, …, then we obtain infinite expansions for f (t) and S(t) which we can denote as

\begin{equation}%\label{dwig} f(t)\approx f_{\infty}(t),\qquad S(t)\approx S_{\infty}(t).\nonumber \end{equation}

The proper meaning for these expansions is that Corollary 1 holds for each m < ∞. The infinite sums f _∞ (t) and S _∞(t) may or may not be pointwise convergent and if so, may or may not converge to f (t) and S(t). We formalise the conditions for equality in the next result whose proof is based on Cauchy’s theorem.

The condition R_m(t) → 0 as m → ∞ is an existence condition for a particular sequence of vertical contours $\{{\mathrm{Re}}(s)=b_{m}^{+}=b_{m}+\varepsilon_{m}\}$ such that R_m(t), the integral over this vertical line, converges to 0 as m → ∞. It is not necessary that this holds for all possible sequences. Typically $b_{m}^{+}\rightarrow\infty$, so in relatively straightforward cases R_m(t) → 0 if it can be shown that $|\mathcal{M}\mathbf{(}b_{m}^{+}+\im y)|$ is uniformly integrable for sufficiently large values of m. Examples where these series correctly converge to f (t) and S(t) are now given.

Example 11. (Example 8 continued.) From (25), the MGF has the form $\mathcal{M}(s)=\mathcal{M}_{5}(s)\times\Gamma(11/2-s)/\Gamma(6-s)$, where $\mathcal{M}_{5}(s)$ contributes the first 5 poles which lead to the term f ₅(t) based on the associated residues. The latter factor $\Gamma(\frac{11}{2}-s)/\Gamma(6-s)$ has an infinite sequence of simple poles at $\{\frac{11}{2},\frac{13}{2},\ldots\}$, which contributes an infinite sum of residues. All conditions of Corollary 2 are shown to hold in Section 5.1.7 of Butler (Reference Butler2019a). After simplification, the residue summation becomes the convergent hypergeometric expansion in

\begin{equation}f(t)=f_{\infty}(t)=f_{5}(t)-\frac{672}{\pi}\re^{-(11/2)t}\,_{4}F_{3}\Big(\frac {1}{2},1,1,\frac{3}{2};\ 3,\frac{7}{2},4;\ \re^{-t}\Big).\nonumber\end{equation}

The exact survival function has S(t) = S _∞(t), which is also expressible in terms of ₄F ₃ hypergeometric functions. Thus, Corollary 2 provides exact expressions for the density and survival functions of Wilks’ likelihood ratio statistic in the intractable setting where k and m are odd-valued, which is the setting in which no exact results have previously been derived. General expressions can be derived for arbitrary odd values of k and m and a computational algorithm is under development to make such computations simple in the general case. Similar infinite residue expansions can be developed for all the other log-likelihood ratio test statistics in normal theory MANOVA as they involve sums of independent log-beta distributions with structure similar to this example; see Anderson (Reference Anderson2003, Chapters 9 and 10).

2.8. Convolutions and finite mixtures

For convolution distributions, the infinite expansions in Corollary 2 may or may not converge as seen in the next three examples. The first two examples convolve a density that has a convergent expansion with a density whose MGF is an entire function.

3.1. Higher-order expansions

If $\mathcal{M}_{R}$ is a meromorphic function in its analytic continuation {Re(s) ≥ b}, then higher-order expansions can be computed subject to the conditions below. In order to simplify the resulting expansions, we assume simple poles in $\mathcal{AC}_{\rm m}^{\rm CD}$ below, which result by assuming that $\mathcal{P}$ admits only one simple pole at r.

• $\mathcal{AC}_{\rm m}^{\rm CD}$. The function $\mathcal{M}_{R}$ admits a sequence of simple real poles b = b ₁ < · · · < b_m or a complex conjugate sequence of simple poles for which b = b ₁ < Re(b _2±) < · · · < Re(b_{m ±}). Apart from these poles, $\mathcal{M}_{R}$ can be continued analytically to {b ≤ Re(s) < b_m + ε ₀} for some ε ₀ > 0.

• $\mathcal{AI}_{\rm mp}^{\rm CD}$. There exists ε_m ∈ (0, ε ₀) and integer p ≥ 1 such that $|\mathcal{M}_{X}(b_{m}^{+}+\im y)|^{p}$ is integrable in y for $b_{m}^{+}=b_{m}+\varepsilon_{m}$. If $b_{m}^{+}>c$ then $\max_{c\,\leq\,x\,\leq b_{m}^{+}}| \mathcal{M} _{X}(x+iN)| \rightarrow0$ as N → ∞.

• $\mathcal{UI}_{\rm m}^{\rm CD}$. If p ≥ 2 in condition $\mathcal{AI}_{\rm mp}^{\rm CD}$ then, for k = 1, …, p − 1, there exists T_k > 0 such that the principal-value integral $\int_{-\infty}^{\infty}\mathcal{M}_{X}(b_{m}^{+} +\im y)^{k}\re^{-\im ty}\sd y$ converges uniformly in t for t > T_k.

Corollary 3. (Higher-order expansions for compound distributions.) Let $\mathcal{P}(z)$ and $\mathcal{M}_{X}(s)$ be as described in Theorem 3.

If f_X satisfies the bounded variation conditions of $\mathcal{BV}^{CD}$, and $\mathcal{M}_{R}$ and $\mathcal{M}_{X}$ satisfy conditions $\mathcal{AC}_{m}^{CD}$, $\mathcal{AI}_{mp}^{CD}$, and $\mathcal{UI}_{m}^{CD}$ above with simple real poles, then

(30) \begin{equation}f_{R}(t) = f_{m}(t)+R_{m}(t)\coloneqq\sum_{j=1}^{m}\re^{-b_{j}t} \frac{-\rho_{-1}}{\mathcal{M}_{X}^{\prime}(b_{j})}+R_{m}(t),\qquad t\rightarrow\infty,\end{equation}

(31) \begin{equation}R_{m}(t) = \frac{1-p(0)}{2\pi}\re^{-b_{m}^{+}t}\int_{-\infty}^{\infty}\mathcal{M}_{R^{+}}(b_{m}^{+}+\im y)\re^{-\im ty}\sd y=o(\re^{-b_{m}^{+}t}),\end{equation}

where ρ ₋₁ is the residue for $\mathcal{P}$ at its single simple pole r and $\mathcal{M}_{R^{+}}$ is given in (29).

Assuming that conditions $\mathcal{AC}_{m}^{CD}$ and $\mathcal{AI}_{mp}^{CD}$ hold, then

(32) \begin{equation}S_{R}(t)=S_{m}(t)+R_{m}^{S}(t)\coloneqq\sum_{j=1}^{m}\re^{-b_{j}t}\frac{-\rho_{-1} }{b_{j}\mathcal{M}_{X}^{\prime}(b_{j})}+o(\re^{-b_{m}^{+}t}),\qquad t\rightarrow\infty,\end{equation}

where $R_{m}^{S}(t)$ is the integral expression in (31) with the additional integrand factor $(b_{m}^{+}+\im y)^{-1}$. With complex conjugate pairs of poles, the same expansions apply when summed over all complex conjugate pairs.

Proof. The proof is given in Section 5.2.3 of Butler (Reference Butler2019a) and broadly follows that of Theorem 3.

If r is rather an m-pole of $\mathcal{P}$ then b ₁, · · ·, b_m are also m-poles, and a more general and complicated version of the corollary can be stated in which the jth terms of (30) and (32) are replaced with terms of the form given in (26) and (28), respectively. Furthermore, if $\mathcal{P}$ admits more than one pole then the residue ρ ₋₁ in (30) and (32) becomes dependent on j and is the residue for the pole of $\mathcal{P}$ at the value $\mathcal{M}_{X}(b_{j})$.

3.2. Negative binomial weights

Suppose that the weight distribution for N is negative binomial (𝔪, p) with q = 1 − p = 1/r and $\mathcal{P}(z)=(p/q)^{\mathfrak{m}} (1/q-z)^{-\mathfrak{m}}$. Then

(33) \begin{equation}\mathcal{M}_{R^{+}}(s)=\frac{p^{\mathfrak{m}}}{1-p^{\mathfrak{m}}}\left[\frac{1}{\{1-q\mathcal{M}_{X}(s)\}^{\mathfrak{m}}}-1\right]\!.\end{equation}

From (33) we see that $\mathcal{M}_{R^{+}}$ can be analytically extended past the value c if c is a pole of $\mathcal{M}_{X}$ but not when it is a branch point. If c is a pole then $\mathcal{M}_{R^{+} }(s)\rightarrow-p^{\mathfrak{m}}/(1-p^{\mathfrak{m}})$ as s → c, so c is a removable singularity for $\mathcal{M}_{R^{+}}$. Additional poles of $\mathcal{M}_{X}$ above c are also removable singularities of $\mathcal{M} _{R^{+}}$, while the zeros of $1-q\mathcal{M}_{X}(s)$ determine 𝔪-poles of $\mathcal{M}_{R^{+}}$ for use in the higher-order expansions of Corollary 3.

Table 3:

Various density and survival function approximations for the true density f (t) and survival S(t) in three settings: using a negative binomial (𝔪, $\tfrac12$) weight distribution with 𝔪 = 3, 2, and 1 along with an extreme value claim distribution. For various approximations, the last accurate digit is boldface.

^a,c Denotes a saddlepoint density or Lugannani and Rice (Reference Lugannani and Rice1980) survival approximation applied directly to $\mathcal{P}\{\mathcal{M} _{X}(s)\}$ and retaining the point mass at t = 0.

^b,d Denotes a saddlepoint density or Lugannani and Rice (Reference Lugannani and Rice1980) survival approximation applied to $\mathcal{M}_{R^{+}}(s)$ and without the point mass and subsequently correcting by using the factor 1 − p(0).

^e Path of steepest descent for computation lies in the directions θ = 0 and 180°.

The Laurent coefficients from the PGF of N are ρ _−𝔪 = (− p/q)^𝔪 and ρ _−k = 0 for k = 1, …, 𝔪 − 1. Values in (27) are

\beta_{-k}=\frac{\rho_{-\mathfrak{m}}}{(\mathfrak{m}-k)!}\frac{d^{\mathfrak{m} -k}}{ds^{\mathfrak{m}-k}}\Big(\frac{1}{\{\mathcal{N}(s)\}^{\mathfrak{m}} }\Big) _{s=b},\qquadk=1,\ldots,\mathfrak{m}-1.

A comparison of f (t) with f ₁(t) and S(t) with S ₁(t) leads to the following observations. For m fixed, accuracy increases with t as expected. If the percentile rank of t is held fixed for R|R ≠ 0, so the accuracy comparison involves columns 2, 4, and 6, then accuracy decreases as the pole order decreases from 𝔪 = 3 to 1.

All 𝔪 terms in expansions (26) and (28) are needed in order to attain the high accuracy seen in the table. For density approximation, if only the leading term in (26) of order O(e^−btt ^𝔪−1) is used in the 𝔪 = 3 and 𝔪 = 2 cases, then the resulting approximate density values are not accurate. This leading term attains the values 0.0²674, 0.0168, and 0.0150 respectively over columns 2–4. For the survival function, the leading term of order O(S _G(𝔪,b)(t)) is not accurate either and attains the respective values 0.0215, 0.0417, and 0.0365.

First-order expansions are least accurate for 𝔪 = 1; however, saddlepoint error corrections also substantially improve accuracy in this case. When t = 3.010 (4.015) then f ₁ + R̂ ₁ = 0.040 904 (0.023 443) and each error correction adds one further digit of accuracy; also S ₁ +$\hat{R}_{1}^{S}=0.073\,6\mathbf{3}9$ $(0.042\,0\mathbf{9}9)$ and error correction adds two (one) additional digits of accuracy.

Saddlepoint density approximations f _SP(t) and f _SP\(t) in Table 3 are based on the approximation (22) applied with and without the point mass p(0) included in the MGF used for the approximation. Thus, we approximate f (t) either directly from $\mathcal{M}_{R}$ to get f _SP(t), or indirectly from $\mathcal{M}_{R^{+}}$ to get f _SP\(t), through f (t) = {1 − p(0)}f_R+ (t) 1_{t>0}. Likewise, S _SP(t) and S _SP\(t) use the Lugannani-Rice approximation (Butler (Reference Butler2007, Section 1.2.1)) to approximate S(t) either directly or else indirectly through S(t) = {1 − p(0)}S_R+(t) 1_{t>0}. Both approximations provide adequate accuracy for all these applications but do not achieve the extraordinary accuracy of f ₁(t) and S ₁(t). Approximations f _SP\(t) and S _SP\(t) have comparable accuracy to their counterparts f _SP(t) and S _SP(t) for m = 3 and 2 and show slightly greater accuracy with 𝔪 = 1. This is explained by the fact that p(0) increases as m decreases so more mass is placed at 0 and the methods that do not remove this point mass encounter greater difficulty and, hence, lose accuracy.

3.2.1. Inversion integration and the analytic continuation of $\mathcal{P}\{\mathcal{M}_{X}(s)\}$

The expansions of Theorems 1 and 2 do not apply directly to the density and survival of R, but rather indirectly by way of the absolutely continuous distribution for R ⁺ = R| R = 0 with MGF $\mathcal{M}_{R^{+} }(s)=[\mathcal{P}\{\mathcal{M}_{X}(s)\}-p(0)]/\{1-p(0)\}$, density f_R+ (t), and survival S_R+ (t); thus, the relations f_R(t) = {1 − p(0)}f_R+(t) and S_R(t) = {1 − p(0)}S_R+ (t) for t > 0 lead to the expansions in Theorem 3 for the distribution of R. In fact, if the density inversion integral in (4) were used with MGF $\mathcal{P} \{\mathcal{M}_{X}(s)\}$, then it would not converge due to the leading term p(0) > 0 in its expansion. Using $\mathcal{P}\{\mathcal{M}_{X}(s)\}$ in the survival inversion integral (54) of Butler (Reference Butler2019a, Section 5.1.4) results in a conditionally convergent principal-value integral but not an absolutely convergent integral; thus, numerical inversion requires integration from ε − Ni to ε + Ni for quite large N. By comparison, both of these integrals are absolutely convergent using $\mathcal{M}_{R^{+}}(s)$ if $|\mathcal{M}_{X}(\varepsilon+\im y)|=O(|y|^{-1-\alpha})$ as y → ∞ for some α > 0. From a practical perspective, however, absolute convergence is typically not enough and α must be large if the numerical inversion is to be accurate using reasonable values of N.

For our example, the ‘exact’ computations for f (t) and S(t) in Table 3 were based on the numerical inversion of $\mathcal{M}_{R^{+}}(s)$ to determine f_R+ (t) and S_R+ (t). The integration was in the convergence strip along the contour line {Re(s) = ŝ ₀} determined by saddlepoint ŝ ₀ ∈ (0, b) which solves $\rd\ln\,\{\mathcal{M}_{R^{+} }(s)\}/\sd s=t$. This leads to an inversion integral expression

(34) \begin{equation}f_{R^{+}}(t)\approx \re^{-\hat{s}_{0}t}\frac{1}{\pi}\int_{0}^{N}{\mathrm{Re}}\{\mathcal{M}_{R^{+}}(\hat{s}_{0}+\im y)\re^{-\im ty}\}\sd y\end{equation}

for sufficiently large N. Survival S_R+ (t) is determined as in (34) by including the additional factor (ŝ ₀ + iy)⁻¹ in the argument of Re{·} in the integrand. We used N = 20 to obtain a 17 digit accuracy based on an inversion integrand which, from (33), has the order

|\mathcal{M}_{R^{+}}(\hat{s}_{0}+\im y)|\sim\frac{p^{\mathfrak{m}}\mathfrak{m} q}{1-p^{\mathfrak{m}}}|\Gamma(1-\hat{s}_{0}-\im y)|=O(|y|^{1/2-\hat{s}_{0} }\re^{-\pi|y|/2}),\qquad |y|\rightarrow\infty,

with ŝ ₀ ∈ (0, 0.55); see 5.11.9 of NIST DLMF (Reference Nist2014).

Figure 1:

Plot of $\mathcal{M}_{R^{+}}(s)\re^{-6s}$ versus s ∈ [0, 5.5] for the 𝔪=3 setting showing the real poles as well as the potential for saddlepoints at critical values in between the integers 1–5.

The analytic continuation of $\mathcal{P}\{\mathcal{M} _{X}(s)\}$. The m-pole b = 0.5571 applies to all five examples in Table 3 since it is a pole for the factor {2 − Г(1 − s)}^−𝔪, which appears in $\mathcal{P}\{\mathcal{M}_{X}(s)\}$ and $\mathcal{M}_{R^{+}}(s)$. The exact remainder terms R ₁(t) and $R_{1}^{S}(t)$ in Table 3 were computed using numerical inversion of $\mathcal{M}_{R^{+}}(s)$ as in (34) but rather passing the vertical contour through ŝ ₁, the closest real saddlepoint above b. In each of the five examples, Cauchy’s theorem, as reflected in the identities f (t) = f ₁(t) + R ₁(t) and $S(t)=S_{1}(t)+R_{1}^{S}(t)$, was confirmed to 17 significant digits through numerical computation of f (t), R ₁(t), S(t), and S ₁(t) using MAPLE®.

In the analytic continuation of $\mathcal{M}_{R^{+}}(s)$ in which Re(s) > b, there are poles at b _2± = 2.581 ± 0.1723i, b ₃ = 4.076, 4.978, 6.004, 6.999, … arising as zeros of 2 − Г(1 − s). For the 𝔪 = 3 example, Figure 1 shows a plot of $\mathcal{M}_{R^{+}}(s)\re^{-6s}$ which captures these real poles as well as potential saddlepoints spaced in between the integers 1 − 5. Saddlepoints are located at 0.2687, 2.196, 2.546, 4.019, and 4.838. For all saddlepoints except 2.546, Figure 1 suggests that $|\mathcal{M}_{R^{+}}(s)\re^{-6s}|$ has a local minimum at each of these saddlepoints along the real line. This is indeed the case as $\mathcal{K}_{R^{+}}^{\prime\prime}>0$ at each of these saddlepoints so the steepest descent (ascent) directions are ±90° (0°, 180°) and the saddlepoint is ‘vertically oriented’. By contrast, the saddlepoint at 2.546 occurs at a local maximum for $|\mathcal{M}_{R^{+}}(s)\re^{-6s}|$ along the real line. This is confirmed by computing $\mathcal{K}_{R^{+}} ^{\prime\prime}(2.546)<0$ so its steepest descent (ascent) directions are θ = 0°, 180° (± 90°) and are horizontally oriented. Thus, the geometry at 2.546 shows a real saddlepoint roughly in between the complex conjugate pair of poles b _2± which have the bearings ±1.370 = ±78.5°.

The extraordinary accuracy of f ₁(t) and S ₁(t) in the examples of Table 3 is partly due to having the poles b _2± lying well above b = 0.5571, so the asymptotic order of error is O(e^{−Re(b ₂₊)t}).

The computations for R ₁(t) and $R_{1}^{S}(t)$ integrate along the respective contours {Re(s) = ŝ ₁} and $\{{\mathrm{Re}}(s)=\mathfrak{\hat{s}}_{1}\}$. In Figure 2 we plot the modulus of the integrand along the former contour leaving saddlepoint ŝ ₁ = 2.196. While the steepest descent curve is tangent at ŝ ₁, it clearly deviates from this vertical line quite quickly since, for y > 0.14, the modulus increases along the vertical line due to the presence of the pole at 2.581 + 0.1723i.

In all the cases for which entries of R̂ ₁(t) and $\hat{R}_{1}^{S}(t)$ are given in Table 3, the next highest saddlepoint above b is vertically oriented and the saddlepoint approximation at this saddlepoint gives a reasonable assessment for the error magnitude in Cauchy’s theorem. In the 𝔪 = 2 case, these saddlepoints were horizontally oriented and this leads to inaccuracy when approximating R ₁(t) and $R_{1}^{S}(t)$. For example, in the case t = 4.658, R̂ ₁(t) = 0.0³146i and $\hat{R} _{1}^{S}(t)=0.0^{3}149\im$ and the imaginary factor i results from the contour integral passing horizontally through the saddlepoints at ŝ ₁ = 0.861 and $\mathfrak{\hat{s}}_{1}=0.851$, respectively, to the right of pole b. Apart from the obviously incorrect imaginary factors, their magnitudes are also quite far from the true values of R ₁(t) and $R_{1}^{S}(t)$. Judging from this example and others, the accuracy of saddlepoint computations as error approximations are in doubt if the saddlepoints are horizontally oriented.

Figure 2:

Plot of ${\mathrm{Re}}\{\mathcal{K}_{R^{+}} (\hat{s}_{1}+\im y)-6(\hat{s}_{1}+\im y)\}$ versus y ∈ [0, 1] for the 𝔪 = 3 setting, where ŝ ₁ = 2.196. The path of steepest descent is tangent at ŝ ₁ (y = 0).

4.1. Higher-order expansions for the Cramér–Lundberg model

Increasingly accurate approximations can be obtained by adding in additional residue terms when $\mathcal{M}_{R^{+}}$ in (39) admits an increasing sequence of simple real poles or complex conjugate pairs. The next result is new and follows from Corollary 3 as shown in Section 5.4.4 of Butler (Reference Butler2019a).

Corollary 5. (Cramér–Lundberg model.) Assume that the conditions of Corollary 4 hold. Suppose, furthermore, that $\mathcal{M}_{R^{+}}$ can be analytically continued as in $\mathcal{AC}_{m}^{CD}$ and that $\mathcal{M}_{X}$ satisfies $\mathcal{X}_{m}^{CL}$.

• $\mathcal{X}_{\rm m}^{\rm CL}$. If b_m > c, then ε_m ∈ (0, ε ₀) exists such that $\max_{c\leq x\leq b_{m}^{+}}|\mathcal{M}_{X}(x+iN)|\rightarrow0$ as N → ∞, where $b_{m}^{+}=b_{m}+\varepsilon_{m}$ and ε ₀ is specified in condition $\mathcal{AC}_{m}^{CD}$.

For real-valued poles, then

(40) \begin{align}f_{R}(t) = f_{m}(t)+R_{m}(t)\coloneqq(1-\rho)\sum_{k=1}^{m}\re^{-b_{k}t}\frac{b_{k}}{\lambda\mathcal{M}_{X}^{\prime}(b_{k})/\sigma-1}+o(\re^{-b_{m}^{+}t}),\end{align}

(41) \begin{align}S_{R}(t) = S_{m}(t)+R_{m}^{S}(t)\coloneqq(1-\rho)\sum_{k=1}^{m}\re^{-b_{k}t}\frac {1}{\lambda\mathcal{M}_{X}^{\prime}(b_{k})/\sigma-1}+o(\re^{-b_{m}^{+}t}),\end{align}

as t → ∞. Here,

(42) \begin{equation}R_{m}(t)=\frac{1-\rho}{2\pi}\re^{-b_{m}^{+}t}\int_{-\infty}^{\infty} \mathcal{M}_{R^{+}}(b_{m}^{+}+\im y)\re^{-\im ty}\sd y,\end{equation}

and $R_{m}^{S}(t)$ is expression (42) with the additional integrand term $(b_{m}^{+}+\im y)^{-1}$. The same expansions hold with complex conjugate pairs of poles when summed over all pairs of poles.

4.2. Heavy-traffic diffusion approximations

Siegmund (Reference Siegmund1979) and Blanchet and Glynn (Reference Blanchet and Glynn2006) provided diffusion or heavy-traffic approximations for the distribution of R in the Sparre Andersen model which are directly related to the expansions above. In particular, the complete asymptotic expansion of Blanchet and Glynn (Reference Blanchet and Glynn2006, Equation 5) is exactly the survival expansion (36) in Theorem 4, and the corrected diffusion approximation of Siegmund (Reference Siegmund1979, Equation 28) is a lower-order approximation to this expansion as discussed further below. The importance of these papers is that they address the computation of the residue coefficient $\beta_{-1}=\operatorname*{Res} [\mathcal{P}\{\mathcal{M}_{L^{+}}(s)\},b],$, a task made quite difficult due to its dependence on $\mathcal{M}_{L^{+}}$. These papers are less relevant in the Cramér–Lundberg setting since explicit values for this residue are available and so the complete expansion of Blanchet and Glynn (Reference Blanchet and Glynn2006, Equation 5) is given in (38) of Corollary 4 with even higher-order terms in (41) of Corollary 5.

Heavy traffic is related to properties of the convex function $\mathcal{K} _{Y}(s)=\ln\mathcal{M}_{Y}(s)$ for the random step size Y. From Siegmund (Reference Siegmund1979), the adjustment coefficient Δ = b − 0 since 0 and b are the unique zeros for $\mathcal{K}_{Y}(s)$. Heavy traffic occurs when the negative drift of the random walk $\mathcal{K}_{Y}^{\prime}(0)\uparrow0$ and Δ = b ↓ 0. In terms of convex $\mathcal{K}_{Y}$, this means that the value and location of its minimum approach 0. The asymptotic regime for heavy traffic has t → ∞, b ↓ 0, and bt → ξ > 0; see Siegmund (Reference Siegmund1979). The expansion in (36) holds under asymptotic regime t → ∞ so it also holds under the more restrictive heavy-traffic regime. Hence, equating the asymptotic expansion in (36) with that in Blanchet and Glynn (Reference Blanchet and Glynn2006, Equation 5) shows their leading terms must be the same within the heavy-traffic asymptotic regime.

This discussion suggests that improvement over the expansions in both Siegmund (Reference Siegmund1979) and Blanchet and Glynn (Reference Blanchet and Glynn2006, Equation 5) is possible in the Sparre Andersen setting by using higher-order expansions comparable to those in Corollary 5. The impediment to this, however, is in determining the additional residues for poles of $\mathcal{P}\{\mathcal{M}_{L^{+}}(s)\}$ further into its analytic continuation.

4.3. Applications and Cramér–Lundberg examples

Wiener–Hopf factorization in the Sparre Andersen model leads to an explicit rational form for $\mathcal{M}_{L^{+}}$ when $\mathcal{M}_{Y}(s)=\mathcal{M} _{X}(s)\mathcal{M}_{T}(-\sigma s)$ is a rational function and B < ∞. Factorization of $1-\mathcal{M}_{Y}(s)$ into monomials in its numerator and denominator fully identifies those monomials belonging to the rational factor that determines $\mathcal{M}_{L^{+}}(s)$, so that $\mathcal{M}_{L^{+}}$ is rational; see Butler (Reference Butler2017, Section B.4.2). Since $\mathcal{P}$ is rational, then $\mathcal{P}\{\mathcal{M}_{L^{+}}(s)\}$ is rational and its partial fraction expansion leads to exact computation of the density and survival function of R ⁺.

By contrast, when $\mathcal{M}_{Y}$ is not rational, then extraction of $\mathcal{M}_{L^{+}}(s)$ is quite difficult except in the Cramér–Lundberg setting where $\mathcal{M}_{L^{+}}$ is the excess life MGF constructed from $\mathcal{M}_{X}$. Tijms (Reference Tijms2003, Section 9.2.2) considered two Cramér–Lundberg examples which use hyperexponential and Erlang claim distributions. In both instances, however, $\mathcal{M}_{X}$ is rational which leads to rational expressions for $\mathcal{M}_{L^{+}}(s)=\{1-\mathcal{M}_{X}(s)\}/(-\mu s)$ and an exact analysis is possible through partial fraction expansion. We avoid these simpler situations by presenting two examples in which $\mathcal{M} _{L^{+}}$ is not rational. In the first example, $\mathcal{M}_{X}$ is an entire function while in the second example, $\mathcal{M}_{X}$ admits a sequence of increasing poles.

In Table 4 we show the accuracy for asymptotic expansions in (40) and (41) of orders m = 1 to m = 7 for the values f_R(2.5) and S_R(2.5). First-order expansion f ₁(2.5) = 0.056 68 attains 4 (3 significant) digit accuracy and this exceeds the accuracy of saddlepoint approximation f _SP\ (2.5) = 0.060 46. First-order expansion S ₁(2.5) = 0.062 071 4 attains similar accuracy and exceeds that of saddlepoint approximation S _SP\ (2.5) = 0.062 17.

Figure 3:

A plot showing the alternating pattern of poles (◊) and saddlepoints (°) for $\mathcal{M}_{R^{+}}(s)\re^{-2.5s}$ for s = x + iy ∈ ℂ with x, y ≥ 0. Lines through saddlepoints indicate the directions of steepest descent.

Table 4:

(Cramér–Lundberg, Raleigh claims). Errors of various residue approximations for density and survival function values f_R(2.5) and S_R(2.5) using a Raleigh (1) claim distribution with mean 1. For a description of other entries, see Table 2.

Exact values for f_R and S_R are displayed to a 30 decimal point accuracy and were computed using numerical integration which roughly follows the path of steepest descent from ŝ ₀ as described in Section 4.4 below.

There are some subtleties concerning how the saddlepoint computations were made. The saddlepoint density value f _SP\ (2.5) = 0.060 46 results if (22) is used only with ŝ ₀ = 0.4753. However, upon viewing the pattern of saddlepoints and poles in Figure 3, additional saddlepoint terms are needed to approximate f (2.5) when based on the method of steepest descents. The contour for this path must be deformed to pass through the saddlepoint ŝ ₀ but also the saddlepoints ŝ _1± = 2.703 ± 3.807i with steepest descent directions ±0.723 ≈ ±41.4°. The sum of all three contributions is 0.060 44 and about the same value. Likewise, in the computation of R̂ ₁, saddlepoint computation requires summing four contributions over ŝ _1± and ŝ _2± = 2.817 ± 2.821i with the latter having steepest descent directions ±2.23 ≈ ±128.0°. This yields the value for R̂ ₁ in Table 4. The value S _SP \(2.5) = 0.062 17 was computed using the Lugannani–Rice approximation, which accommodates for the pole at s = 0 using the Bleistein (Reference Bleistein1966) method, but which is not based upon following the path of steepest descent. Saddlepoint approximation using (23) and summing over saddlepoints $\mathfrak{\hat{s}}_{0}=0.6504$ and $\mathfrak{\hat{s}}_{1\pm}=2.679\pm3.794\im$ leads to 0.058 17 which, as expected, is less accurate due to not accommodating the pole at s = 0.

Table 5:

(Cramér–Lundberg, truncated extreme value claims). Errors of various residue approximations for density and survival function values f_R(4.5) and S_R(4.5) using a truncated extreme value claim distribution. For a description of other entries, see Table 2.

Example 19. (Cramér–Lundberg, truncated extreme value claims.) Assume that the claim distribution has the density (1 − e⁻¹) exp (− t − e^−t) 1_{t>0} with mean μ = 1.260. Its MGF is derived in Section 5.4.8 of Butler (Reference Butler2019a) as

(43) \begin{equation}\mathcal{M}_{X}(s)=\frac{1}{1-\re^{-1}}\frac{1}{1-s}\,_{1}F_{1}(1-s;2-s;-1),\end{equation}

where ₁F ₁ denotes the confluent hypergeometric function. Expression (43) admits simple poles at s 𢘈 {1, 2, . . .} at which $\mathcal{M}_{R^{+}}(s)$ is analytic. Interspersed between these poles are simple zeros for the denominator of $\mathcal{M}_{R^{+}}(s)$ which occur at

\begin{align*}b=b_{1}=0.4632<b_{2}=2.226<2.896<4.023<4.995<6.001=b_{6}.\end{align*}

We choose λ = 1/μ and σ = 2 so that $\rho=\tfrac12$. In between the poles, there are real saddlepoints for $\mathcal{M}_{R^{+}}(s)\re^{-4.5s}$ which are used to approximate error terms for f_R(4.5). They occur at

\hat{s}_{1}=2.049<2.735<3.906<4.312(-)<5.958<6.055(-)=\hat{s}_{6},

where the negative sign (no sign) indicates a saddlepoint that is oriented horizontally (vertically) so the path of steepest descent is in directions 0, π (± π/2). Slightly different saddlepoints apply when dealing with the survival inversion.

In Table 5 we show that increasingly accurate approximations can be obtained by adding in additional residue terms from these simple poles. The leading residue terms f ₁(4.5) = 0.030 175 8 and S ₁(4.5) = 0.065 154 3 achieve a 5 and 4 decimal point accuracy, and this exceeds that of the saddlepoint approximations f _SP\ (4.5) = 0.032 54 and S _SP\(4.5) = 0.065 21, which achieve a 2 and 4 decimal point accuracy. Each additional residue term contributes 2 or 3 more accurate digits in the expansions. The saddlepoint computation of errors R̂_m(4.5) and $\hat{R}_{m}^{S}(4.5)$ provides a reliable guide for expansion accuracy by capturing the correct order for error except, as previously seen, when horizontally oriented saddlepoints are involved.

Exact computations of f_R(4.5) and S_R(4.5) were computed by roughly following the path of steepest descent as described below.

Figure 4:

(Raleigh claim distribution). A plot of ${\mathrm{Re}} \{\mathcal{M}_{R^{+}}(\hat{s}_{0}+\im y)\re^{-\im ty}\}$ versus y (solid line) and its asymptotic expression μ ⁻¹(1 − ρ)y ⁻¹ sin (ty) versus y (dashed line) shows that the former plot resembles a slowly dampened sine function. Here, t = 2.5, μ = 1, and $\rho=\tfrac12$.

4.4. Steepest descent inversion required

In both examples of Section 4.3, numerical confirmation for the accuracy of various approximations for f_R(t) and S_R(t) required choosing inversion contours that are more specific to the examples than those in (34). In fact, for all Cramér–Lundberg examples, the integrand for the inversion in (34) is of order O{y ⁻¹ sin(ty)} as y → ∞ and not absolutely integrable. A proof is given in Section 5.4.5 of Butler (Reference Butler2019a).

Theorem 5. (Cramér–Lundberg inversion contour.) Let the inversion integral (34) for f_R+ (t) use an arbitrary contour {Re(s) = x} for x ∈ [0, b). Then, the integrand in (34) for the resulting inversion integral which computes f_R+ (t) has the asymptotic expansion

(44) \begin{equation}\re^{-xt}\frac{1}{\pi}{\mathrm{Re}}\{\mathcal{M}_{R^{+}}(x+\im y)\re^{-\im ty} \}\sim\frac{(1-\rho)\re^{-xt}}{\pi\mu}\frac{\sin(ty)}{y},\qquad y\rightarrow \infty.\end{equation}

Thus, the left-hand side of (44) is conditionally integrable in y but not absolutely integrable. The integrand for inversion of S_R+ (t) if x ≠ 0 has the expansion

(45) \begin{equation}\re^{-xt}\frac{1}{\pi}{\mathrm{Re}}\Big\{ \frac{\mathcal{M}_{R^{+}}(x+\im y)\re^{-\im ty}}{x+\im y}\Big\} \sim\frac{(1-\rho)\re^{-xt}}{\pi\mu}\,\frac {\cos(ty)}{y^{2}},\qquad y\rightarrow\infty,\end{equation}

so the left-hand side of (45) is absolutely integrable but still converges very slowly.

In Figure 4 we show how quickly the inversion integrand for the density conforms to its oscillatory asymptotic expansion in (44) when x = ŝ ₀. Such slowly dampened oscillations are a major impediment to accurate density and survival function inversion. To avoid such difficulties in our numerical computations, we deform the inversion contour so it becomes a contour ray which ultimately approximates the path of steepest descent in direction θ ∈ [0, π/2]. Instead of integrating from ŝ ₀ − i∞ to ŝ ₀ + i∞ as in (34), we choose ray angle θ above the real axis and ray angle −θ below the real axis. Thus, our deformation of the integration is along three lines: from ŝ ₀ − iπ + e^−iθ∞ to ŝ ₀ − iπ to ŝ ₀ + iπ to ŝ ₀ + iπ + e^iθ∞. As shown in Section 5.4.6 of Butler (Reference Butler2019a), this provides the following inversion expression.

Theorem 6. (Cramér–Lundberg deformed contour.) Let absolutely continuous variable R ⁺ have density f_R+ (t) on (− ∞, ∞) which is locally of bounded variation at t. Suppose that its MGF $\mathcal{M}_{R^{+}}(s)$ converges in the right half-plane on {0 ≤ Re(s) < b} for b > 0. Furthermore, assume that

\begin{equation}\max_{\theta\leq\phi\leq\pi/2}|\mathcal{M}_{R^{+}}(\hat{s}_{0}+\im\pi+N\re^{\im\phi })|\rightarrow0, \qquad N\rightarrow\infty,\nonumber\end{equation}

and also assume that $\mathcal{M}_{R^{+}}$ is analytic in the sector {ŝ ₀ + iπ + re^iϕ : r ≥ 0, ϕ ∈ [θ, π/2]}. Then, the inversion integral (4) for f_R+ (t) which passes through the saddlepoint ŝ ₀ may be deformed so that

(46) \begin{align}f_{R^{+}}(t) & = \re^{-\hat{s}_{0}t}\frac{1}{\pi}\int_{0}^{\pi}{\mathrm{Re}}\{\mathcal{M}_{R^{+}}(\hat{s}_{0}+\im y)\re^{-\im ty}\}\sd y\label{invertray}\end{align}

(47) \begin{align}&\jump +\re^{-\hat{s}_{0}t}\frac{1}{\pi}\int_{0}^{\infty}\operatorname{Im} [\mathcal{M}_{R^{+}}(\hat{s}_{0}+\im\pi+r\re^{\im\theta})\exp\{ \im(\theta -t\pi)-tr\re^{\im\theta}\} ] \sd r. \label{invertray2}\end{align}

The inversion for S_R + (t) is the same as in (46) and (47) but includes the additional integrand factors (ŝ ₀ + iy)⁻¹ and (ŝ ₀ + iπ + re^iθ)⁻¹ inside the curly and square braces of (46) and (47), respectively.

Contour deformation, as in Theorem 6, leads to extraordinary gains in computational efficiency over (34) for numerical inversion. This is demonstrated in Example 20, which provides a setting in which the ultimate steepest descent direction is θ ≈ 50.6°, and also in Example 21 where the ultimate direction is θ = 0.

The problem with conventional inversion algorithms is their insistence on using vertical contours that stay in the convergence domain, which ensures that convergence is slow, rather than attempting to follow the path of steepest descent which is ultimately in some direction θ ∈ [0, π/2 − ε) and which takes the contour outside the convergence domain. If we deform the contour to crudely follow this path, as given in the integration formula of Theorem 6, then Example 21 shows that numerical integration in MAPLE® achieves a 30 digit accuracy using a contour that is $12\%=100(\frac{60}{500})\%$ of the length used in following the vertical contours. Thus, our integration after deformation far outperforms the truncated integration along the vertical contour and is, furthermore, more stable and faster per unit of contour length. There seems to be no sensible reason for applying sophisticated convergence acceleration methods to slowly converging integrals in the convergence domain in the first place when much faster convergence can be achieved by deforming the contour to pass into the analytic continuation and to roughly follow the direction of steepest descent.

From a practical point of view, however, numerical integration is not even needed in this example since the first three residue terms suffice as a substitute for exact computation. We know this without making exact computations since the suggested saddlepoint approximations in the analytic continuation accurately provide the correct orders for these errors.