\section{Introduction}
\vskip.1truecm{\bf Motivation:}
The  problem of
approximating  distributions   based on  empirical data like moments
is one of the bread and butter problems of applied probability.

In
  risk theory and related {\bf first passage
  applications}  (queueing,
mathematical finance, ...) there is special interest in
approximating
densities of nonnegative random variables by {\bf affine combinations of exponentials}, also
called {\bf GHE densities (generalized hyper-exponential)}, or, more generally, by  {\bf matrix exponential distributions} (which allow for "collision of exponents"). One reason for that is that this class
captures the asymptotic behavior in the important "light tails case".

There exists  already a quite extensive literature, based on inverting
the explicit Pollaczek-Khinchine formula for the Laplace transform, which assumes
  however  complete knowledge of the model, for example knowledge of all its moments.

  Since  input data are never certain, it is interesting to  develop {\bf approximations based on finitely many moments} (the  coefficients of the power series expansion of the Laplace transform
around $0$).

Some well-known
such approximations in risk theory and queueing   are  the    {\bf Renyi, De Vylder, Gamma},  and {\bf Whitt -Ramsay} approximations,
obtained by {\bf fitting one, two or three moments}, and the minimal (but arbitrarily high) order {\bf three moments fitting formulas of Bobbio, Horvath and Telek}
\cite{bobbio2005matching}.
Also useful    is the {\bf Cram\'er-Lundberg}
approximation \begin{eqnarray}} \def\eeq{\end{eqnarray} \label} \def\fr{\frac} \def\im{\item{CLa}\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q(x) \sim C e^{- \gamma}  \def\d{\delta} \def\z{\zeta}  \def\th{\theta x}\eeq
 where $\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q(x)$ denotes the ruin probability \eqref{eq:ru} and $-\gamma}  \def\d{\delta} \def\z{\zeta}  \def\th{\theta$ is the so called adjustment coefficient (i.e. the
largest negative root of the Cram\'er-Lundberg equation \eqref{eq:CL}), which captures  the asymptotic
behavior in the case of light tail claims.
All these approximations may be derived from the explicit Laplace transform provided by the Pollaczek Khinchine  formula \eqref{PK}.

Producing higher order approximations in the  "intermediate regime" when a finite,  but  larger number of moments is known, seems a  very challenging problem.


\begin{Qu} In view of the scarcity of
approximations fitting  more than three moments, it is natural to ask {\bf what are the difficulties blocking the development of high  order  moments based approximations}?
\end{Qu}


Below, we examine
this question in the context of  first passage theory for {\bf spectrally negative Levy processes}.

 \vskip.1truecm{\bf Pad\'e, two-point Pad\'e and other rational   approximations}. It turns out  that most of the approximations currently used are Pad\'e   approximations of Laplace transforms, and that higher order approximations
 are quite easy to
obtain (for example using   the Mathematica command {\bf PadeApproximant}), since the conversion from moments to a {rational Laplace approximation}
 requires  only solving a linear system.
\begin{comment}
 \begin{ex}} \def\eeXa{\end{ex} \label} \def\fr{\frac} \def\im{\item{ex:22}  For the classic Cram\'er-Lundberg  model, the
$(1,2)$ Pad\'e approximation for $\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q^*(s)$,   in terms of the stationary excess  moments, is:
\begin{eqnarray*}&& \Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q^*(s)\approx\frac{N}{D},\\ &&N=s \left(-3 \tilde} \def\B{\bar} \def\H{\hat m_2^3+4 \tilde} \def\B{\bar} \def\H{\hat m_1 \tilde} \def\B{\bar} \def\H{\hat m_3 \tilde} \def\B{\bar} \def\H{\hat m_2-\tilde} \def\B{\bar} \def\H{\hat m_1^2 \tilde} \def\B{\bar} \def\H{\hat m_4\
\right)-4 \tilde} \def\B{\bar} \def\H{\hat m_3 \tilde} \def\B{\bar} \def\H{\hat m_1^2+6 \tilde} \def\B{\bar} \def\H{\hat m_2^2 \tilde} \def\B{\bar} \def\H{\hat m_1\\&&D= s^2\left(-3 \tilde} \def\B{\bar} \def\H{\hat m_2^3+4 \tilde} \def\B{\bar} \def\H{\hat m_1 \tilde} \def\B{\bar} \def\H{\hat m_3 \tilde} \def\B{\bar} \def\H{\hat m_2-\frac{1}{2} \th  \tilde} \def\B{\bar} \def\H{\hat m_4 \tilde} \def\B{\bar} \def\H{\hat m_2+\frac{2 \th  \tilde} \def\B{\bar} \def\H{\hat m_3^2}{3}-\tilde} \def\B{\bar} \def\H{\hat m_1^2 \tilde} \def\B{\bar} \def\H{\hat m_4\right)\\&&+s\left(-4 \tilde} \def\B{\bar} \def\H{\hat m_3 \tilde} \def\B{\bar} \def\H{\hat m_1^2+6 \tilde} \def\B{\bar} \def\H{\hat m_2^2 \tilde} \def\B{\bar} \def\H{\hat m_1-\th  \tilde} \def\B{\bar} \def\H{\hat m_4 \tilde} \def\B{\bar} \def\H{\hat m_1+2 \th  \tilde} \def\B{\bar} \def\H{\hat m_2 \tilde} \def\B{\bar} \def\H{\hat m_3\right)+6 \th  \tilde} \def\B{\bar} \def\H{\hat m_2^2-4 \th  \tilde} \def\B{\bar} \def\H{\hat m_1 \tilde} \def\B{\bar} \def\H{\hat m_3\\&&= s N +\frac{1}{6} \th   \left(36 \tilde} \def\B{\bar} \def\H{\hat m_2^2-3 s \left(s \tilde} \def\B{\bar} \def\H{\hat m_4-4 \tilde} \def\B{\bar} \def\H{\hat m_3\right) \tilde} \def\B{\bar} \def\H{\hat m_2+4 \
s^2 \tilde} \def\B{\bar} \def\H{\hat m_3^2-6 \tilde} \def\B{\bar} \def\H{\hat m_1 \left(4 \tilde} \def\B{\bar} \def\H{\hat m_3+s \tilde} \def\B{\bar} \def\H{\hat m_4\right)\right)\end{eqnarray*}
or, in terms of the claims moments,
\begin{eqnarray*}&& N=2
   \lambda  \left(45 m_4 m_2^2-60 m_3^2 m_2+s \left(20 m_3^3-30 m_2
   m_4 m_3+9 m_2^2 m_5\right)\right)
   \\&&D=\left(40 \lambda  m_3^3-60 \lambda  m_2 m_4 m_3+12
   p m_5 m_3-15 p m_4^2+18 \lambda  m_2^2
   m_5\right) s^2\\&&+\left(90 \lambda  m_4 m_2^2-120
   \lambda  m_3^2 m_2+36 p m_5 m_2-60 p m_3
   m_4\right) s-240 p m_3^2+180 p m_2 m_4\\&&= s N +
   3 p \left(-80 m_3^2+4 s \left(s m_5-5 m_4\right)
   m_3-5 s^2 m_4^2+12 m_2 \left(5 m_4+s
   m_5\right)\right)\end{eqnarray*}
   \eeXa

\end{comment}
It is also quite easy to produce Pad\'e approximations with specified limiting
   behavior of the Laplace transform at $\infty$, so called {\bf two-point Pad\'e approximations} -- see Example \ref{ex:exp}.

   Other
rational approximations of interest
 are those
   minimizing in least squares
 sense the sum of the coefficients -- see Beylkin and Monz\'on \cite{beylkin2005approximation}.
 Another interesting class are
 "{\bf Tijms approximations}", which    try to  incorporate moments fitting with including the exact Cram\'er-Lundberg  asymptotics  \eqref{CLa}
as dominant term --see \cite{willmot1998class}.
These may also be obtained by using  a Pad\'e  approximation with a  prescribed pole.

   \vskip.1truecm{\bf The   admissibility of  Pad\'e    approximations in probability}.
While in principle a great tool due to their easiness of implementation, and their convergence for large
$n,$  Pad\'e   approximations (and variations, like two point Pad\'e) applied to Laplace transforms in probability
have the
drawback of the difficulty to  guarantee   {\bf "admissible inverses"}, i.e. {\bf nonnegative densities}
   and {\bf non-increasing survival functions},   when fitting three moments or more is desired
   (note that for fitting two moments $m_1,m_2$ of a nonnegative random  variable, the admissible Gamma approximation
   $$f(x)\sim \tilde} \def\B{\bar} \def\H{\hat f(x)=\fr{(\mu x)^{\alpha} \def\b{\beta -1}}{\G(\alpha} \def\b{\beta)} \mu e^{- \mu x}, \; \alpha} \def\b{\beta=\frac{m_1^2}{m_2-m_1^2}, \mu=\frac{m_1}{m
   _2-m_1^2}$$   provides an easy solution).









  Even  ensuring the
{\bf nonnegativity  of combinations of exponentials with fixed given rates} is quite a difficult
question (since this
involves an infinite number  of constraints), still open nowadays,  except for two exponentials
(when nonnegativity of $\tilde} \def\B{\bar} \def\H{\hat f(0)$ and nonnegativity of the coefficient of the asymptotically dominant exponent are clearly necessary and sufficient), and for three exponentials \cite{dehon1982geometric}.



\begin{ex} \label} \def\fr{\frac} \def\im{\item{e:Bot}
Consider  the example due to Harris \cite{harris1992note} (see also
\cite[Ch. 5.4]{fackrell2003characterization}) \begin{eqnarray*} f(t)= 2 e^{-t} -
6 e^{-2 t} + 6 e^{-3t}=2 e^{-t} - 3( 2 e^{-2 t}) +  2 ( 3
e^{-3t})\\
\Leftrightarrow \bar F(t)= 2 e^{-t} - 3 e^{-2 t} +  2
e^{-3t}\end{eqnarray*} with {\bf canonical coordinates} $(2,-3,2),$
and Laplace transform


{\small \begin{eqnarray*} f^*(s)=\frac{2 \left(s^2+2 s+3\right)}{(s+1)
   (s+2) (s+3)} =\fr 23 \fr{6}{(s+1)(s+2)(s+3)}-\fr 13 \fr{6}{(s+2)(s+3)}+\fr 23 \fr{3}{s+3}\end{eqnarray*}}




Since the "Coxian coordinates" $(2/3,-1/3,2/3),$ produced by the partial fractions decomposition above  
are not nonnegative, this is not a phase-type distribution of order $3$.

However, by an admissibility criteria for  {\bf combinations of three negative
exponentials} due to
\cite{dehon1982geometric}\symbolfootnote[4]{with exponents equal to $i=1,2,3$,
this is
$-\w_2 \leq 2 \sqrt{\w_1 \w_3},$ or, after normalization
$$\w_1 + 1 \leq \w_3 \leq \w_1 + 1 + 2 \sqrt{\w_1 }.$$}
(see also
\cite{fackrell2003characterization,bean2008characterization} for a
criterion which allows colliding exponents), we know  this is a proper density. In fact, it is a
phase-type density of order $4$ -- see \eqref{e:4}
 \eeXa

\vskip.1truecm{\bf Admissibility by phase-type
representations}.
 One  "lucky case" in which admissibility is automatic is when one has obtained somehow
      {\bf any phase-type representation} $PH(\bff \alpha} \def\b{\beta, A)$
      with $\bff \alpha} \def\b{\beta$
      a {\bf probability vector} and $A$ a {\bf subgenerator} matrix
      (satisfying $A_{ij} \geq 0$ for $i \neq j$ and $A \bff 1 \leq \bff 0$).
       In that case, nonnegativity follows from the
probabilistic interpretation of $f(t)$ as the density of the
absorbtion time of the corresponding Markovian semigroup.

      However,  determining when  a phase-type representation exists is
        again notoriously difficult, the so
      called {\bf positive realization}  problem of systems theory
      \symbolfootnote[4]{Finding the minimum possible order of such a representation is even harder, and known currently only for three moments fitting representations).
Attesting further to the difficulty of providing admissible phase-type approximations are several interesting recent
approaches, like the recursive {\bf minimal order three moments
fitting} formulas of \cite{bobbio2005matching},
and the {\bf EM algorithm} approach
 \cite{willmot2011risk}.}.

In example \ref{e:Bot}, it is possible to show that a phase-type of minimal order $4$ is available,
  by  a  recursive approach  of
decomposing the   admissibility domain as union of higher order
admissibility polytopes associated to {phase-type representations}
-- see  \cite{colm1993triangular}.%
One  phase-type representation is:
\begin{eqnarray}} \def\eeq{\end{eqnarray} \label} \def\fr{\frac} \def\im{\item{e:4} A=\left(\begin{array}{llll}
 -1 & 1 & 0 & 0 \\
 0 & -2 & 2 & 0 \\
 0 & 0 & -3 & 3 \\
 0 & 0 & 0 & -4
\end{array}
\right), \q \alpha} \def\b{\beta=(1/2,0,0,1/2)\eeq

BUTools http://webspn.hit.bme.hu/~telek/tools/butools/butools.html
only obtains a representation of order $5$,
further illustrating the difficulty of this problem.



\begin{Qu} Currently,   no algorithmic approach for  testing admissibility of combinations of more than four given negative exponentials (or four terms M\"untz polynomials) is available \symbolfootnote[4]{a solution might  however be possible by the approach of Faybusovich \cite[Thm 4,5]{faybusovich2002self}, who offers a general representation of the
Koszul-Vinberg  characteristic function  \cite{seeger2012epigraphical} of the positivity cone generated by any
Chebyshev system, as Pfaffian of a matrix of multiple integrals (and taking logarithm
yields a  "self-concordant barrier" function).}.

 \end{Qu}

The question of providing admissible  approximations
(with
non-specified rates) fitting more than four given values of a Laplace transform (for example moments) is a priori even more difficult. However, two remarkable exceptions in which this the challenging admissibility problem was solved are the {\bf minimal order three moments fitting admissible
    approximations} of Bobbio, Horvath and Telek
\cite{bobbio2005matching}, obtained by a recursive approach on the order, and the non-minimal ones of Johnson and Taaffe \cite{johnson1989matching}, an outcome of the classical moments theory, which  work for any number of moments.
















    \begin{comment}(their "admissibility" and asymptotic correctness are however hard to enforce.

  The basic question of rational approximation is how to chose  the {\bf interpolation nodes} where we fit our function and whether to impose constraints on the {\bf poles of the denominator}   -- see \cite{avdis2007power} for a nice review.




  A further
generalization of the De Vylder and Ramsay's approximations is to look
for


The simplest possibility is the classic Hermite-Pad\'e approach, which provides a  multiple fit of the function of interest  at $0$, and,  when applied to  a Laplace transform,  is  clearly equivalent to  equating {\bf integer moments}.
 This approach
  may be traced  back to   the Chebyshev-Markov-Stieltjes moments
 problem, and may be efficiently implemented nowadays using Rutishauser's {\emph  qd quotient-difference} scheme
  -- see
 \cite{henrici1993applied,wenocur1989approximating,van1990moment} for  superb accounts.


After obtaining any {\bf RLT} (rational
 approximation of the
 Laplace transform),
     partial fractions and Laplace  inversion will result
in a
 combination of exponentials, which may also be
 represented as a {\bf matrix exponential function}\symbolfootnote[4]{This is  more general than moments fitting, which corresponds to the particular case  of the {\bf Hermite-Pad\'e}
   method of fitting the MacLaurin expansion of the Laplace transform around $0$.}.
   Thus, RLT/matrix exponential approximations may be viewed as a method
   of Laplace inversion.
\end{comment}



   \vskip.1truecm{\bf Contents and contributions}.
   This problem was motivated by the desire to provide new admissible matrix exponential approximations in ruin theory.

Necessary ruin theory background, including the Pollaczek Khinchine  transform, is reviewed in Section \ref{s:PK}. Section \ref{s:L} presents the "key characters" in our ruin application, the aggregate loss and
 its moments  -- see \eqref{expexp}, \eqref{Lmts}.

First  order Pad\'e approximations of the Pollaczek Khinchine  transform due to Renyi and DeVylder
are reviewed in Section \ref{s:Ap}. We also provide  here a new approximation, for the perturbed Cram\'er-Lundberg  model -- see Theorem \ref{t:per}.
The remarkable Johnson-Taaffe approximations are reviewed in Section \ref{s:JT}, Theorem \ref{t:JT}. Using their approach, we may "update" a second order approximation due to
 Ramsay to make it work for arbitrary  claims having three moments -- see Theorem \ref{t:Ram}.

In the case of random sums,    it is possible to apply the Johnson-Taaffe approach
both to the individual summands -- we call this a "Ramsay-type approximation", and directly to the sum -- yielding  "Beekman-Bowers-type
approximations".
Theorem \ref{t:new} in Section \ref{s:comp} provides a  comparison between these
two methods, by comparing their explicit {\bf 3 moments
Johnson-Taaffe orders}.

Section \ref{s:prop}  discusses  two-point Pad\'e approximations --see
Theorem \ref{t:Ram2}, whose admissibility is an open problem for the
moment.  Finally, in Section~\ref{sec:num} we provide some numerical
examples.












\begin{comment}

















\end{comment}




\sec{Ruin theory background \label} \def\fr{\frac} \def\im{\item{s:PK}}
{\bf The  perturbed Cram\'er
Lundberg risk process} models the reserves of an insurance compan
by:

\begin{equation}\label{CLmod} X(t) = u
+ {c} \def\w{\omega} \def\q{\qquad} \def\D{\Delta}\def\Si{\Sigma \,} t- S(t) + \s W(t), \q S(t)=\sum_{k=1}^{N_\l(t)} Z_k
\end{equation}
 used in collective risk theory to describe
the surplus $X=\{X(t), t\ge0\}$ of an insurance company. Here, \begin{enumerate}}  \def\BI{\begin{itemize}
\im $u$ is the initial capital, \im $c} \def\w{\omega} \def\q{\qquad} \def\D{\Delta}\def\Si{\Sigma \, t$ represents the premium
income up to time $t$, \im $Z_k$ are i.i.d. positive random
variables representing the claims made, with cumulative distribution
function  and density denoted by $F(x)$ and $f(x),$ and (some)
moments denoted by $m_i, i=1,2,...$,  \im $ N=\{N_\l(t), t\ge0\}$ is
an independent Poisson process with intensity $\l$ modeling the
times at which the claims occur, and \im $ W(t)$ is an infinite variation  spectrally negative perturbation, for example a standard Wiener  motion, and $\sigma>0$  is  a scale parameter.\end{enumerate}}   \def\EI{\end{itemize}




Since the  jumps of  $X$ are all negative, the moment generating
function ${\mathbb E}} \def\Z{{\mathbb Z}} \def\R{{\mathbb R}[\te{s X(t)}]$ exists for all $s\ge0$ and $t\ge 0$,
and is log-linear in $t.$  The symbol/Laplace exponent/{cumulant generating function}
$\k(s)$ is defined by \begin{eqnarray}} \def\eeq{\end{eqnarray} \label{eq:psi} \k(s)=\log
\left({\mathbb E}} \def\Z{{\mathbb Z}} \def\R{{\mathbb R}_0 \te{s X(1)}\right) = \boxed{s\left} \def\ri{\right} \def\i{\infty( c} \def\w{\omega} \def\q{\qquad} \def\D{\Delta}\def\Si{\Sigma - \l
\bar{F}^*(s)+\tilde} \def\B{\bar} \def\H{\hat \k(s) \ri)},\eeq where $\bar{F}^*(s)$ denotes the Laplace
transform of the survival function of the claims, and $s \tilde} \def\B{\bar} \def\H{\hat \k(s) $ is the symbol of the perturbation $\s W(t)$. For example, in the Wiener case $\tilde} \def\B{\bar} \def\H{\hat \k(s) = \frac{\s^2}{2} s$.




Let $T $ be the first passage time of a stochastic process $X(t)$
below $0$: { \begin{equation*} T : = \inf\{t\ge0: X(t) <
0\
.\end{equation*} The objects of interest in classical ruin theory are the
 "finite-time" and "eventual" ruin probabilities

\begin{eqnarray*} &&\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q(t,u)= P_u [T
 \leq t],
 \q \boxed{\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q(u)= P_u [T < \infty]=P_0[L > u]},  \label} \def\fr{\frac} \def\im{\item{eq:ru}\end{eqnarray*}
 where $L=-\underline{X},$ also called maximal aggregate loss, is the negative of all-time infimum of the process \eqref{CLmod},
 started from $0$.




 \begin{comment}


{\bf The integro-differential equation.} The Markov property of the
process \eqref{CLmod} implies that $\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q(x)$ is a homogeneous
solution of the semigroup generator:
\begin{eqnarray} \label} \def\fr{\frac} \def\im{\item{IDE} && \frac{\s^2}{2} \,\frac{\partial^2 \Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q
(x)}{\partial x^2}+c} \def\w{\omega} \def\q{\qquad} \def\D{\Delta}\def\Si{\Sigma \,\frac{\partial \Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q (x)}{\partial x}  +
 \l \int_{0}^{\infty} \Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q(x-z) F_Z (dz) -\l  \Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q (x) =0,
 \\&& {\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q}(x)=1, \text{ if} \; {x < 0}, \q  \Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q (0)=1, \text{ if} \; {\s > 0}
 \nonumber
\end{eqnarray}


\end{comment}

The ultimate ruin probability $\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q(u)$ for \eqref{CLmod}  is not identically $1$ iff
the Levy drift/profit rate
 \begin{eqnarray}} \def\eeq{\end{eqnarray} \boxed{p:= c} \def\w{\omega} \def\q{\qquad} \def\D{\Delta}\def\Si{\Sigma - \l m_1 >0:=\l m_1 \th} \eeq
 is positive, in which case  adding the condition $\lim_{x \to
 \infty} \Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q(x)=0$  determines it uniquely.




\vskip.1truecm{\bf The Pollaczek-Khinchine formula.} Taking
Laplace transform of the Kolmogorov equation for the {\bf
 ultimate ruin probabilities} of the perturbed
Cram\'er-Lundberg model  yields the
Pollaczek-Khinchine formula:
 \begin{eqnarray}\label} \def\fr{\frac} \def\im{\item{PK}\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q^*(s)=\frac{1}{s}-
\frac{\k'(0)}{\k(s)}=\frac{\rho} \def\s{\sigma(1-f_e^*(s))+ \tilde} \def\B{\bar} \def\H{\hat \k(s)/c} \def\w{\omega} \def\q{\qquad} \def\D{\Delta}\def\Si{\Sigma}{s(1- \rho} \def\s{\sigma f_e^*(s)+\tilde} \def\B{\bar} \def\H{\hat \k(s)/c} \def\w{\omega} \def\q{\qquad} \def\D{\Delta}\def\Si{\Sigma)}:=\frac{1}{s}(1-
\f(s))\end{eqnarray} The first expression, in terms of the symbol
$\k(s)$ of the Levy process involved, is valid  for all spectrally negative L\'{e}vy processes,
-- see for example \cite{kyprianou2006introductory}). The second
emphasizes the fact that the result in the Cram\'er-Lundberg case
depends only on the "equilibrium density" of the claims, defined
by \begin{eqnarray}} \def\eeq{\end{eqnarray} \label} \def\fr{\frac} \def\im{\item{eq:sex} f_e(x):=\bar{F}(x)/m_1,\eeq the estimation of which may be a
convenient starting point. The third expression is equivalent to $ \label} \def\fr{\frac} \def\im{\item{eq:rud}\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q(x)=\int_x^\infty \tilde} \def\B{\bar} \def\H{\hat \p(u) du, \q  x>0$ where  $\tilde} \def\B{\bar} \def\H{\hat \p(u)$ is the {\bf density of the continuous part of the distribution of the aggregate loss}
$L$
 (while $\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q(x) $ is the  survival function of $L$).
The Pollaczek Khinchine  formula for $\f(s)$
 \begin{eqnarray}} \def\eeq{\end{eqnarray} \label} \def\fr{\frac} \def\im{\item{eq:L} \boxed{ \f(s)={\mathbb E}} \def\Z{{\mathbb Z}} \def\R{{\mathbb R} e^{-s L}=1- s \Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q^*(s)=\frac{\k'(0)}{\k(s)/s}

 \frac{1-\rho} \def\s{\sigma}{1-\rho} \def\s{\sigma f_e^*(s)+ \tilde} \def\B{\bar} \def\H{\hat \k(s)/c} \def\w{\omega} \def\q{\qquad} \def\D{\Delta}\def\Si{\Sigma}}. \eeq




In the case $\s=0$,
a beautiful probabilistic interpretation of the (nonperturbed) Pollaczek Khinchine  formula \eqref{eq:L} was
 discovered (independently)
by Benes, Kendall and Dubordieu, by expanding the denominator into a geometric series:
\begin{eqnarray}\label} \def\fr{\frac} \def\im{\item{PKf}&&
\boxed{\f(s)=\fr{1-\rho} \def\s{\sigma}{1- \rho} \def\s{\sigma
\widehat f_e(s)}=\sum_{k=0}^\infty (1-\rho} \def\s{\sigma) \rho} \def\s{\sigma^k \widehat f_e^{k}(s)}\; \; \; \end{eqnarray} where $f_e(x)=\bar F(x)/m_1$ is the
equilibrium distribution of $Z_i$. This reveals that  $\f(s)$  is the Laplace
transform of
a geometric sum of convolutions of the equilibrium/stationary distribution,
i.e. $\p(x)=(1-\rho} \def\s{\sigma) \sum_{n=0}^\infty \rho} \def\s{\sigma^k f_e^{(*,k)}(x)$\symbolfootnote[4]{Another interpretation of
$\f(s)$ is  as
   Laplace transform of the stationary waiting time of the M/G/1 queue --see for example \cite[Thm VIII.5.7]{asmussen2003applied}.}, which may be visualized by examining the "ladders" of the paths (the amounts by which the process $-c} \def\w{\omega} \def\q{\qquad} \def\D{\Delta}\def\Si{\Sigma t + S(t)$ jumps to new maxima).




This is the so called Pollaczek Khinchine  "ladder decomposition"    $L=\sum_{n=1}^N L_n$ of the
 "maximal aggregate loss random variable" $L$, with
 $N$  a geometric r.v.
$\Pr[N=k]=(1-\rho} \def\s{\sigma) \rho} \def\s{\sigma^{k}, k=0,1,..., \rho} \def\s{\sigma=\frac{\l m_1}{c} \def\w{\omega} \def\q{\qquad} \def\D{\Delta}\def\Si{\Sigma}$,
representing the number of ladders.




\begin{rem}} \def\eeR{\end{rem} {\bf The ladder decomposition, $\s>0$}.
In the Brownian perturbed case, a  beautiful probabilistic interpretation of the  Pollaczek Khinchine  formula
  \eqref{eq:L} was recently discovered by Dufresne-Gerber and rederived  in an elementary way by  Kella, by rewriting \eqref{eq:L} as:
 \begin{eqnarray}} \def\eeq{\end{eqnarray} \label} \def\fr{\frac} \def\im{\item{eq:Kel} \f(s)=\fr{1-\rho} \def\s{\sigma}{1- \rho} \def\s{\sigma f_e^*(s)+ \tilde} \def\B{\bar} \def\H{\hat \k(s)/c} \def\w{\omega} \def\q{\qquad} \def\D{\Delta}\def\Si{\Sigma}=\fr{1}{1+ \tilde} \def\B{\bar} \def\H{\hat \k(s)/c} \def\w{\omega} \def\q{\qquad} \def\D{\Delta}\def\Si{\Sigma} \fr{1-\rho} \def\s{\sigma}{1- \rho} \def\s{\sigma f_e^*(s) \fr{1}{1+ \tilde} \def\B{\bar} \def\H{\hat \k(s)/c} \def\w{\omega} \def\q{\qquad} \def\D{\Delta}\def\Si{\Sigma}}\eeq
  reflects  the fact that $L$ is an independent sum of a "first creep at the current infimum" (which in the diffusion case is an exponential of rate $\frac{2 c} \def\w{\omega} \def\q{\qquad} \def\D{\Delta}\def\Si{\Sigma}{\s^2}$)  and of  an alternating geometric sum
 of "compound Poisson ladders and further creeps" when $\s>0$
 -- see \cite[Fig 2]{dufresne1991risk}, \cite{kella2011class}.


In the case that the Brownian perturbation $\sigma W(t)$ is replaced
by a {\bf general spectrally negative perturbation} $Y$ with non-zero expectation $E[Y(1)]\ge 0$,
a similar ladder decomposition holds true. Let
$$\kappa^*(s) = s \tilde} \def\B{\bar} \def\H{\hat \k(s)$$
 denote the Laplace exponent of \, $Y$. Note that
by the Wiener-Hopf factorisation--see e.g. \cite{kyprianou2006introductory}, $\tilde} \def\B{\bar} \def\H{\hat \k(s)$ is Laplace exponent of the possibly killed downward ladder process
(and $s$ is the Laplace exponent of the up-crossing ladder process),  that is, $\tilde} \def\B{\bar} \def\H{\hat \k$ takes the form
$$\tilde} \def\B{\bar} \def\H{\hat \k(s) = c + \delta s + \int_{(0,\infty)} (1-\te{-s x}) \nu(dx),$$
where the killing rate $c$ and the drift $\delta$ are non negative constants and the L\'{e}vy measure $\nu$ satisfies
the integrability condition
$$\int_{(0,\infty)}[ 1\wedge x] \nu(dx) < \infty.$$
Then, the formula \eqref{eq:Kel} still holds,
 providing a decomposition of $L$ as an independent sum of the {\bf increment of the ladder process of the perturbation $Y$ at an independent exp$(c)$- random time}
and the geometric sum of further  such increments and ``compound Poisson ladder height increments''.
In particular, if the perturbation is given by a completely asymmetric stable process, i.e.
$\kappa^*(s) = s^\alpha$, $\alpha\in(1,2)$, then $\kappa_-^*(s)=s^{\alpha-1}$
and we identify $\frac{1}{1 + \kappa_-^*(s)/c}$ as the Laplace transform of
the non-negative random variable $\tilde Y^{\alpha-1}_{e(c)}$,
where $Y^{\alpha-1}$ is a stable subordinator with parameter $\alpha-1$.
\eeR



\begin{rem}} \def\eeR{\end{rem}
 Note that the aggregate loss $L$ is the mixture of a discrete mass of
$\f(\infty)=1-\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q(0)=1-\rho} \def\s{\sigma$ at $0$, and of a  continuous random
variable. Letting $\tilde} \def\B{\bar} \def\H{\hat \f(s)$ denote the Laplace transform of the
density $\tilde} \def\B{\bar} \def\H{\hat \p(u)$, note the decomposition \begin{eqnarray}} \def\eeq{\end{eqnarray} \label} \def\fr{\frac} \def\im{\item{e:fn}
&&\f(s)=1-\rho} \def\s{\sigma+ \rho} \def\s{\sigma \tilde} \def\B{\bar} \def\H{\hat \f(s) \; \Leftrightarrow \; \p(x) = (1-\rho} \def\s{\sigma) \d_0(x) + \rho} \def\s{\sigma \tilde} \def\B{\bar} \def\H{\hat
\p(x)\nonumber} \def\bs{\bigskip  \\&& \tilde} \def\B{\bar} \def\H{\hat \f(s):=\fr{\f(s)-\f(\infty)}
{\f(0)-\f(\infty)}=\fr{\f(s)-(1-\rho} \def\s{\sigma)} {\rho} \def\s{\sigma}=(1-\rho} \def\s{\sigma)\fr{ f_e^*(s)}{1- \rho} \def\s{\sigma
 f_e^*(s)}, \eeq
 where we denoted by $\p(u)$ the inverse Laplace transform of $\f(s)$, given by the generalized function
 $\p(u)= \tilde} \def\B{\bar} \def\H{\hat \p(u) + (1-\rho} \def\s{\sigma) \d_0(u).$




The behavior at $\infty$  of $\f(s)$ distinguishes between
the nonperturbed ($\s=0$) and perturbed case ($\sigma>0$):
\begin{eqnarray}\label} \def\fr{\frac} \def\im{\item{lim} \lim_{s \to \infty}
  \f(s)=\lim_{s \to \infty}
  1- s \Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q^*(s)=1- \Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q(0)=\begin{cases}
  1-\frac{\l m_1}{c} \def\w{\omega} \def\q{\qquad} \def\D{\Delta}\def\Si{\Sigma}=1-\rho} \def\s{\sigma=\frac{p}{c} \def\w{\omega} \def\q{\qquad} \def\D{\Delta}\def\Si{\Sigma}, \; & \s=0\\
  0, \; & \s>0\end{cases}\end{eqnarray}

  The Laplace transform
   $\f(s)$
is an essential quantity in the
theory of L\'{e}vy processe
 



\eeR

\begin{rem}} \def\eeR{\end{rem} The roots of the denominator $\k(s)=0$ in the Pollaczek Khinchine  formula determine the asymptotic behavior of ultimate ruin probabilities. More generally, an important role is played by the roots of the Cram\'er Lundberg equation \begin{eqnarray}} \def\eeq{\end{eqnarray}\label} \def\fr{\frac} \def\im{\item{eq:CL}
\k(s)=q, q >0.\eeq \eeR

\begin{rem}} \def\eeR{\end{rem} In
practice, the true distribution of the claims (and interarrival
times) is of course unknown, and since the Pollaczek Khinchine  formula requires
this knowledge,  it should be viewed more as a theoretical than a practical
tool. It may be argued that  the most reliable information
available in insurance data is  contained in the first few integer
moments, and thus it seems natural to start building approximations by
fitting moments, or, equivalently, by Pad\'e interpolation of the
Laplace transform at the origin.
\eeR





\section
  The moments of the aggregate loss \label} \def\fr{\frac} \def\im{\item{s:L}}

     From now on, we will assume a classic Brownian perturbation. As an alternative  to Ramsay's approximation of the ladder distribution, we  may approximate   directly the aggregate loss distribution. One advantage is that the second moment of the aggregate loss satisfies automatically the second order representability constraint $c_v \geq 1/2$. This allows focusing on the third moment
     constraint.




  Consider   the
expansion
 \begin{eqnarray}\label} \def\fr{\frac} \def\im{\item{mtsexp} &&\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q^*(s) =\frac{
 \n_{2,\s} /2-  \n_3 s/6+...}{p+  \n_{2,\s} s/2- \n_3
s^2/6+...} \Leftrightarrow \\&&\f(s) =\frac{p}{p+  \n_{2,\s} s/2- \n_3
s^2/6+...}\nonumber} \def\bs{\bigskip \\&&=\frac{\th}{1 + \th- \left} \def\ri{\right} \def\i{\infty(1-(\frac{m_2}{m_1}+ \frac{\s^2}{2 \l m_1})s/2+ \frac{m_3}{m_1} s^2/6-\frac{m_4}{m_1}
s^3/4!+...\ri)}\nonumber} \def\bs{\bigskip\\&&=\frac{1-\rho} \def\s{\sigma}{1-\rho} \def\s{\sigma\left} \def\ri{\right} \def\i{\infty(1-(\frac{m_2}{m_1}+
\frac{\s^2}{2 \l m_1})s/2+ \frac{m_3}{m_1} s^2/6-\frac{m_4}{m_1}
s^3/4!+...\ri)}\nonumber} \def\bs{\bigskip
\end{eqnarray}where $\n_i=\l m_i, i=0,1,2,...$ are
the moments of the L\'{e}vy measure, and $ \n_{2,\s}= \n_2+
{\s^2}.$




\begin{rem}} \def\eeR{\end{rem} When $\s=0,$ the  expression in the last parenthesis of
\eqref{mtsexp} $$\widehat f(s):=1-\frac{m_2}{2 m_1}s+ \frac{m_3}{3 m_1}
s^2/2-\frac{m_4}{4 m_1} s^3/3!+...={\mathbb E}} \def\Z{{\mathbb Z}} \def\R{{\mathbb R} [e^{-s L_i} |\{L_i
>0\}]$$  has moments $\tilde} \def\B{\bar} \def\H{\hat
m_i=\frac{m_{i+1}}{(i+1) m_1},$  identifying $L_i$ as  the famous
equilibrium/stationary excess/ladder variable generated by $Z_i$,  with
density $f_e(x)=\bar{F}(x)/m_1$, and {\bf stationary excess moments} $\tilde} \def\B{\bar} \def\H{\hat
m_i=\frac{m_{i+1}}{(i+1) m_1}$. \eeR




\begin{rem}} \def\eeR{\end{rem} Let $l_i= {\mathbb E}} \def\Z{{\mathbb Z}} \def\R{{\mathbb R} L^i, \l_i= l_i/i!, i \geq 1$ denote the
  moments and
  "factorially reduced moments" obtained by "normalizing" with respect
  to the
exponential moments of the maximal aggregate loss. These may be easily obtained, either by the  recursion equivalent of the Pollaczek Khinchine  formula:
\begin{eqnarray*} \th \l_n=\tilde} \def\B{\bar} \def\H{\hat \mu_n +\sum_{k=1}^{n-1} \tilde} \def\B{\bar} \def\H{\hat \mu_k \l_{n-k}, \q n \geq 1, \tilde} \def\B{\bar} \def\H{\hat \mu_k:=\fr{\tilde} \def\B{\bar} \def\H{\hat m_{k}}{k!}=\fr{\mu_{k+1}}{m_1}\end{eqnarray*}
or by expanding \eqref{mtsexp} in power series.
The first factorially reduced moments are:
\begin{eqnarray}\label} \def\fr{\frac} \def\im{\item{Lmts} && \l_1={\mathbb E}} \def\Z{{\mathbb Z}} \def\R{{\mathbb R} L= \frac{\tilde} \def\B{\bar} \def\H{\hat{m}_1}{\th}+ \fr{\s^2}{ p}= \frac{\l {m}_2+\s^2}{2 p},  \\&&
 \l_2={\mathbb E}} \def\Z{{\mathbb Z}} \def\R{{\mathbb R} L^2/2!= \frac{\tilde} \def\B{\bar} \def\H{\hat{m}_2}{2 \th}+  \l_1^2,
 \;
   \l_3={\mathbb E}} \def\Z{{\mathbb Z}} \def\R{{\mathbb R} L^3/3!= \frac{\tilde} \def\B{\bar} \def\H{\hat{m}_3}{3!\th}
  + 2  \frac{\tilde} \def\B{\bar} \def\H{\hat{m}_1 \tilde} \def\B{\bar} \def\H{\hat{m}_2/2}{ \th^2}  +   \frac{\tilde} \def\B{\bar} \def\H{\hat{m}_1^3}{\th^3}
  \nonumber} \def\bs{\bigskip \\&&
   \l_4={\mathbb E}} \def\Z{{\mathbb Z}} \def\R{{\mathbb R} L^4/4!= \frac{\tilde} \def\B{\bar} \def\H{\hat{m}_4}{4! \th}
  +  \frac{ 2\tilde} \def\B{\bar} \def\H{\hat{m}_1 \tilde} \def\B{\bar} \def\H{\hat{m}_3/3!  +  (\tilde} \def\B{\bar} \def\H{\hat{m}_2/2)^2
  }{\th^2}  +  \frac{3 \tilde} \def\B{\bar} \def\H{\hat{m}_1^2 \tilde} \def\B{\bar} \def\H{\hat{m}_2/2}{\th^3} +
  \frac{\tilde} \def\B{\bar} \def\H{\hat{m}_1^4  }{\th^4},\nonumber} \def\bs{\bigskip
  \end{eqnarray}
and the mass of the continuous part is $\l_0:=\Pr [L>0]=\rho} \def\s{\sigma=\frac{1}{1+\th}.$


\eeR

\begin{rem}} \def\eeR{\end{rem}

When $\s=0$, the factorially reduced moments $\l_k$ admit an
interesting decomposition:
\begin{eqnarray*}\label} \def\fr{\frac} \def\im{\item{expexp} &&\f(s)=(1-\rho} \def\s{\sigma)\frac{1}{1-\rho} \def\s{\sigma \widehat f_e(s)}=\sum_{k=0}^\infty (1-\rho} \def\s{\sigma) \rho} \def\s{\sigma^k  (f_e^{*}(s))^k=
\\&&1+
 \sum_{k=1}^\infty \l_k (-s)^k =
1+\sum_{k=1}^\infty (-s)^k \sum'_{i_1+i_2+...=k} \p_{i_1}\p_{i_2}...
\end{eqnarray*}where $\sum'$ denotes sum over all
decompositions of $k$ as a sum, and  \begin{eqnarray*} \p_i=\frac{\l m_{i+1}}{p
(i+1)!}=\th^{-1}\frac{\tilde} \def\B{\bar} \def\H{\hat m_{i}}{ i!}=\th^{-1}\widetilde {\mu}_i>0,\end{eqnarray*}
where $\tilde} \def\B{\bar} \def\H{\hat m_{i} :=\frac{m_{i+1}}{(i+1)m_1}, \mu_i:=\frac{m_i}{i!}$
and $\tilde} \def\B{\bar} \def\H{\hat \mu_i:=\frac{\tilde} \def\B{\bar} \def\H{\hat m_i}{i!}$. \eeR



 \begin{rem}} \def\eeR{\end{rem} The moments of the conditioned continuous r.v. $L|L>0$,  necessary for applying  certain results from the literature, may be obtained by dividing by $\rho} \def\s{\sigma=\fr{1}{1+ \th}$. \eeR

  \begin{rem}} \def\eeR{\end{rem} Note that the corresponding moments $l_i$ may also be viewed as moments of the
  ruin function:
  \begin{eqnarray*}  \Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q_0&:=&\l_1=\int_0^\infty \Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q(x) dx=\frac{\tilde} \def\B{\bar} \def\H{\hat{m}_1}{\th}, \; \Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q_1:=\l_2=
   \int_0^\infty x \Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q(x) dx=\frac{\tilde} \def\B{\bar} \def\H{\hat{m}_2}{2 \th}+   (\frac{\tilde} \def\B{\bar} \def\H{\hat{m}_1}{\th})^2,\\
  \Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q_{k}&:=& \fr{\int_0^\infty x^{k}\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q(x) dx}{k!}=  \l_{k+1},...\end{eqnarray*}
   \eeR






 \section{Renyi, De Vylder, and a new simple approximation for ruin probabilities

 \label} \def\fr{\frac} \def\im{\item{s:Ap}}

{\bf Approximations  of ultimate ruin probabilities}.  The problem of approximating ultimate ruin probabilities $\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q(u)$
for the Cram\'er Lundberg model \eqref{CLmod} using data on the
distribution $F(u)$ of the claims is a classic of applied
probability, dating back before 1900.  Its roots may be traced back
to the Danish mathematician TN. Thiele, who
founded  the first insurance company, Hafnia, who
is also  the {\bf inventor of
cumulants}  and of {\bf Thiele continued fractions} useful for
Laplace transform inversion \cite{cohen2007numerical}.

In this section, we make  the
observation  that the exponential mixture approximations recalled in the introduction are particular  cases of Pad\'e approximations of  Laplace transforms (a theme already present  in Thiele's research preoccupations). We also provide a new simple
approximation in this vein for the perturbed Cram\'er-Lundberg  model in Theorem \ref{t:per}.
  



    \begin{rem}} \def\eeR{\end{rem} Concerning  our application which involves random sums, we make the observation that   it is possible to use  Pad\'e, Johnson-Taaffe and any other rational approximations of Laplace transforms at three levels: \begin{enumerate}}  \def\BI{\begin{itemize} \im for the density of the claims, based on the estimates of $m_i, i \geq 1$
\im for the equilibrium density of the claims, using the
estimates of the equilibrium moments $\tilde} \def\B{\bar} \def\H{\hat m_i, i \geq 1,$ and  the profit rate $p $.
The second level is intuitively  superior to the first, since the
equilibrium density is monotonically decreasing, even when the claims
density isn't \symbolfootnote[4]{For example, this  gives rise to a smaller JT index, as illustrated  in Example \ref{ex:uni}}. We will call this Ramsay type approximation.
  \im for
the aggregate loss density transform $\f(s)$ (or,
equivalently, the ccdf  transform $\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q^*(s)$), using directly the moments
$\l_{i+1}, i \geq 0,$ given in \eqref{Lmts}. This amounts
to {\bf working directly with the Pollaczek-Khinchine formula},
instead of approximating its denominator, and so intuitively,
should be better,   at least under certain conditions.\symbolfootnote[4]{More precisely, the third method is expected to be  better in the case of light tails claims and heavy traffic, while the second is expected to be  better in the case of heavy tails claims and light traffic. The fact that both methods are
 better sometimes is illustrated in Theorem \ref{t:new}, from the point of view of yielding a smaller JT index. Identifying the domains within which   methods two and three are preferable in "boundary cases"  is not an easy task. Note also that a mixture
 of the two has been also proposed \cite{sakurai2004approximating}.}. We will call this Beekman-Bowers type approximation.
 \end{enumerate}}   \def\EI{\end{itemize}
  \eeR





 One  Pad\'e
  approximation we consider here is: \begin{eqnarray}\label} \def\fr{\frac} \def\im{\item{Pad}
  (\f)^{(n)}(s)={\mathcal Pade}_{(m,n)}(\f(s))=
  {\mathcal Pade}_{(m,n)}([\frac{p}{\k(s)/s}]_{N}
 
 
\end{eqnarray} where ${\mathcal Pade}_{(m,n)}$ denotes the "classic Pad\'e
approximation" based on the Taylor series around $0$, where $m=n$
or $m=n+1$ for the classic/perturbed Cram\'er Lundberg  process,
respectively, and where $[f(s)]_N$ denotes truncation of a power
series to its first $N$ terms, with $N=m+n+1$  (for "theoretical
models" where an expression for $\f(s)$ is available, we may also
take $N=\infty$). The first case is applicable to the classic, and the
second to the perturbed Cram\'er Lundberg  model.
 As mentioned, the motivation of \eqref{Pad} is that
   the classic
 De Vylder approximation, 
 is  precisely
 the one point Pad\'e approximation of $\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q^*(s),$ around $0$, of orders $(n-1,n)$,
 with $n=1$.

 A second class of
Pad\'e
  approximations we experiment with  is: \begin{eqnarray}\label} \def\fr{\frac} \def\im{\item{Pad1}
(\f)^{(n)}(s)={\mathcal Pade}_{((m,n,m_1,n_1))}(\f(s))= {\mathcal
Pade}_{((m,n,m_1,n_1))}([\frac{p}{\k(s)/s}]_N)
\end{eqnarray}where ${\mathcal Pade}_{((m,n,m_1,n_1))}$ denotes a
two point Pad\'e approximation based on the power series around
$0$ and $\infty$. These are indispensable when dealing with the perturbed model.













\begin{ex} The simplest Ramsay-type approximation is   the {\bf one moment Renyi  exponential
approximation of the equilibrium density}  (which may also be viewed as a Pad\'e $(0,1)$ approximation of the aggregate loss density, which imposes also  the correct limiting behavior $\lim_{s
\to \infty}
      s \Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q^*(s)=\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q(0)=\rho} \def\s{\sigma=1/(1+\th)$ of the Laplace transform  at $\infty$).
 This  amounts to
    looking for an approximation of the form
    $$\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q^*(s)\approx
    \frac{\rho} \def\s{\sigma}{s+b_0} \Leftrightarrow f_e^*(s)\approx
    \frac{b_0/(1-\rho} \def\s{\sigma)}{s+b_0/(1-\rho} \def\s{\sigma)},$$
    where $ f_e^*(s)$ denotes  the Laplace transform of the stationary excess density of the claims \eqref{eq:sex} (the two being related by the Pollaczek Khinchine  formula \eqref{PK}).

      Fitting one moment yields $\frac{ 1-\rho} \def\s{\sigma}{b_0}={\tilde} \def\B{\bar} \def\H{\hat m_1},$ where $ \tilde} \def\B{\bar} \def\H{\hat m_1=\frac {m_2}{2  m_1} $ (the first moment
    of the stationary excess density), and  
  
    \begin{eqnarray}} \def\eeq{\end{eqnarray} \label} \def\fr{\frac} \def\im{\item{e:Re} \Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q^*(s)&\approx&
   \frac{\rho} \def\s{\sigma}{\frac{ 1-\rho} \def\s{\sigma}{\tilde} \def\B{\bar} \def\H{\hat m_1}+ s}, \q \Leftrightarrow \q
     \Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q(x)
    \approx \rho} \def\s{\sigma  e^{-  x (1-\rho} \def\s{\sigma)/\tilde} \def\B{\bar} \def\H{\hat m_1},\eeq
    which is also correct in the limit $\rho} \def\s{\sigma \to 1$ when $\tilde} \def\B{\bar} \def\H{\hat m_1$ exists see \cite{kalashnikov1997geometric}, \cite[(31)]{grandell2000simple}.



\label} \def\fr{\frac} \def\im{\item{ex:exp}
\end{ex}

\vskip.1truecm{\bf De Vylder's  exponential approximation}.
 The simplest Beekman-Bowers type approximation is De Vylder's \eqref{e:DVL}, one of the most popular
   approximations for ruin probabilities, due to its simplicity and  asymptotic correctness \cite{grandell2000simple}, despite of  its
   its being expressed  in terms  of only the {  first three moments} of the claims.





\vskip.1truecm{\bf Derivation of De Vylder as a two  moments Pad\'e
approximation}. De Vylder's approximation  was obtained originally by
equating the first three moments  of   the original Cram\'er
Lundberg process with those of   a new process
 with
 { exponential claims}, and with different arrival intensity and premium rate.
We check now, using the moments \eqref{Lmts} of the aggregate loss density,
that De Vylder's
formula coincides with the Pad\'e $(0,1)$ approximation around $0$ of
the  Laplace transform  $\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q^*(s)$ of the ruin probabilities.   We start  by expanding in power series the numerator and denominator of the
Pollaczek-Khinchine formula \eqref{PK}:

\begin{eqnarray*} &&\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q^*(s)=\fr{ \l m_2 /2- \l m_3 s/6+...}{p+ \l m_2 s/2- \l m_3
s^2/6+...}\approx \fr{a}{s + \alpha} \def\b{\beta} \\&& \Leftrightarrow a s (p+ \l m_2 s/2- \l m_3
s^2/6+...) \\&&\approx (s+ \alpha} \def\b{\beta) (\l m_2 s/2- \l m_3 s^2/6+...)\\&&\Leftrightarrow
 \begin{cases} a p= \alpha} \def\b{\beta \l m_2 /2,\\  a  m_2 /2= m_2 /2-\alpha} \def\b{\beta m_3 /6 \end{cases} \Leftrightarrow
 \begin{cases}
 a= \fr{3 \l m_2^2}{3 \l  m_2^2+2 p  m_3},\\
 \alpha} \def\b{\beta=\fr{6 p m_2}{ 3 \l
m_2^2+2 p  m_3} \end{cases}\end{eqnarray*}
and so
\begin{eqnarray}} \def\eeq{\end{eqnarray} \label} \def\fr{\frac} \def\im{\item{e:DVL}
\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q_{DV}^*(s)=\fr{a}{s + \alpha} \def\b{\beta} =
\fr{ \fr{3 \l m_2^2}{3 \l  m_2^2+2 p  m_3}}{s+\fr{6 p m_2}{ 3 \l
m_2^2+2 p  m_3}}=\frac{3 \l m_2^2}{s(3 \l  m_2^2+2 p  m_3)+6 p m_2}\eeq

\begin{comment}


 \vskip.1truecm{\bf Derivation of De Vylder by equating
moments.}
 Alternatively, De Vylder's approximation could  be obtained by
 fitting the moments of order $0$ and $1$ of $\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q(x)$ with those of a defective exponential approximation
$$\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q(x)
    \approx  a e^{- x /\b }$$
   or, equivalently, {\bf by equating the moments of order
   $1$ and $2$ of the density $\p(x)$ of the aggregate loss} $L$ \eqref{Lmts}.

  Solving $ \begin{cases}a {\b} =\frac{\tilde} \def\B{\bar} \def\H{\hat{m}_1}{\th}\\
   a {\b^2}=\frac{\tilde} \def\B{\bar} \def\H{\hat{m}_2}{2 \th}+   (\frac{\tilde} \def\B{\bar} \def\H{\hat{m}_1}{\th})^2\end{cases}$
   yields:
   \begin{eqnarray}} \def\eeq{\end{eqnarray} \label} \def\fr{\frac} \def\im{\item{e:DV} && a=\fr{1}{1+\th \frac{\tilde} \def\B{\bar} \def\H{\hat m_2 }{2 \tilde} \def\B{\bar} \def\H{\hat m_1^2}}, \q \b=\fr{\l_2}{\l_1}=\frac{\tilde} \def\B{\bar} \def\H{\hat{m}_2}{2 \tilde} \def\B{\bar} \def\H{\hat m_1 }+   \frac{\tilde} \def\B{\bar} \def\H{\hat{m}_1}{\th}=\frac{\tilde} \def\B{\bar} \def\H{\hat{m}_2}{2 \tilde} \def\B{\bar} \def\H{\hat m_1 } \frac{1}{1- a }
   \\&&\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q(x)
    \approx {a}
   e^{- x (1- a) /\widetilde{\tilde} \def\B{\bar} \def\H{\hat m_1}  }=\frac{1}{1+\tilde} \def\B{\bar} \def\H{\hat {\th}}
   e^{- x {\tilde} \def\B{\bar} \def\H{\hat \th} /(1+ \tilde} \def\B{\bar} \def\H{\hat \th)\widetilde{\tilde} \def\B{\bar} \def\H{\hat m_1} },
   \eeq
   where $ a=(1+ \tilde} \def\B{\bar} \def\H{\hat
\th)^{-1}, \; \tilde} \def\B{\bar} \def\H{\hat \th=\th \frac{\tilde} \def\B{\bar} \def\H{\hat m_2 }{2 \tilde} \def\B{\bar} \def\H{\hat m_1^2}  =\th \frac{2 m_1 m_3}{3  m_2^2}$, and
$$\tilde} \def\B{\bar} \def\H{\hat m_i= \fr{m_{i+1}} {(i+1) m_1}, i \geq 1$$ are the
{\bf stationary excess/equilibrium moments}.

\begin{rem}} \def\eeR{\end{rem} In conclusion, the Pad\'e $K=1, L=1$ approximation of
the ruin density \eqref{eq:rud} is

\begin{eqnarray}} \def\eeq{\end{eqnarray} \label} \def\fr{\frac} \def\im{\item{e:DVf}
\f_{DV}(s)=1- s \Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q_{DV}^*= \frac{2 p \left(3 m_2+s m_3\right)}{3 s \lambda
   m_2^2+6 p m_2+2 p s m_3}\eeq

Note that $\f_{DV}(\infty)=\frac{p }{p+\frac{3 \lambda  m_2^2}{2 m_3}}$
is a lower(upper) bound of the survival probability $\overline{\rui}(0)=\fr{p}{p + \l m_1},$
according to whether $\fr{m_3 m_1}{m_2^2} <(>) \fr 32$. \eeR
\end{comment}

\begin{comment}

The Laplace transform  is simpler  in terms of equilibrium moments:
   \begin{eqnarray*}&& \fr{a \b}{s \b +1}=\frac{\frac{\tilde} \def\B{\bar} \def\H{\hat{m}_1}{\th}}{s(\frac{\tilde} \def\B{\bar} \def\H{\hat{m}_2}{2 \tilde} \def\B{\bar} \def\H{\hat m_1 }+   \frac{\tilde} \def\B{\bar} \def\H{\hat{m}_1}{\th}) +1}= \\
   &&=\frac{1}{\frac{2 {\th} m_1}{m_2} + s (\frac{2
   {\th} m_1 m_3}{3 m_2^2}+1)}=\frac{1}{\tilde} \def\B{\bar} \def\H{\hat \th + 1}\frac{1}{{\tilde} \def\B{\bar} \def\H{\hat \th} /(1+ \tilde} \def\B{\bar} \def\H{\hat \th)\widetilde{\tilde} \def\B{\bar} \def\H{\hat m_1}+ s}
   =\frac{1}{{\tilde} \def\B{\bar} \def\H{\hat \th} /(1+ \tilde} \def\B{\bar} \def\H{\hat \th)\widetilde{\tilde} \def\B{\bar} \def\H{\hat m_1}+ s(\th \frac{\tilde} \def\B{\bar} \def\H{\hat m_2 }{2 \tilde} \def\B{\bar} \def\H{\hat m_1^2} + 1)}
   \end{eqnarray*}



 \end{comment}

 \begin{Qu}


  Identify  domains within which   methods two and three  are "preferable" in some sense, in the simplest case of exponential approximations, i.e. {\bf compare the Renyi and De Vylder approximations}. \end{Qu}





\begin{comment}
 \vskip.1truecm{\bf A  well behaved Laplace inversion method, and
Erlang mixtures with common order}. One inversion benchmark  is
the Post-Widder method:
 \begin{eqnarray*} f_n(t)&=&P_n (f)(t)=\fr{(-1)^{n-1} (\fr{n}{t})^n }{(n-1)!}  \H f^{(n-1)}(\fr{n}{t})\\&=&\fr{(\fr{n}{t})^n}{(n-1)!}  \int_0^\infty e^{- \fr{n}{t} u} u^{n-1}  f(u) du, t \geq 0, \; \;  n \geq 1\end{eqnarray*}
 which was shown in
 \cite{jagerman1978inversion} to have  many ideal "shape preserving"
 properties (the approximants share with the approximated function positivity, bounds, monotonicity, diminishing of variation, ...).


  Note however that   this "ideal method"
 
   requires complete knowledge of the Laplace
  transform.


 To turn this into a practical method  with empirical information

 one could use empirical density approximations in the second formula\symbolfootnote[4]{other statistical estimation methods, using
 empirical Laplace transform \cite{horvathexplicit},
  mixture estimation
 via the
 EM algorithm \cite{lee2010modeling}, or  Bayesian estimation \cite{ausin2004bayesian}
 have been shown to work
quite well with large data sets. Here, we   investigate however
the analytical moments
 fitting approach, leaving for future work
   comparisons to the
 statistically inspired approaches.} leading to approximations via
 Erlang mixtures with common
order.



\end{comment}




\begin{comment}
 For example, at order $2$, the Tijms approximation is
\begin{eqnarray*} \bar{F}^*(x)= C e^{-\gamma}  \def\d{\delta} \def\z{\zeta}  \def\th{\theta x} + (C-\rho} \def\s{\sigma) e^{- \gamma}  \def\d{\delta} \def\z{\zeta}  \def\th{\theta_1 x},\end{eqnarray*}
 combining the exact Cram\'er-Lundberg asymptotics $C e^{-\gamma}  \def\d{\delta} \def\z{\zeta}  \def\th{\theta x}$  of light tailed geometric sums with one more exponential term with $\gamma}  \def\d{\delta} \def\z{\zeta}  \def\th{\theta_1=(\rho} \def\s{\sigma -C)(\fr{\rho} \def\s{\sigma \tilde} \def\B{\bar} \def\H{\hat m_1}{1-\rho} \def\s{\sigma}- \fr{C}{\gamma}  \def\d{\delta} \def\z{\zeta}  \def\th{\theta})^{-1}$, obtained by fitting the first moment.  Asymptotic correctness require assuring that
 $\gamma}  \def\d{\delta} \def\z{\zeta}  \def\th{\theta_1 > \gamma}  \def\d{\delta} \def\z{\zeta}  \def\th{\theta \;  \Leftrightarrow \; 0 < (1-\rho} \def\s{\sigma- \gamma}  \def\d{\delta} \def\z{\zeta}  \def\th{\theta \tilde} \def\B{\bar} \def\H{\hat m_1)({\rho} \def\s{\sigma {\gamma}  \def\d{\delta} \def\z{\zeta}  \def\th{\theta} \tilde} \def\B{\bar} \def\H{\hat m_1}- {C}(1-\rho} \def\s{\sigma)),$
 and admissibility  requires  assuring that
 $C \gamma}  \def\d{\delta} \def\z{\zeta}  \def\th{\theta+ (C-\rho} \def\s{\sigma) \gamma}  \def\d{\delta} \def\z{\zeta}  \def\th{\theta_1 \geq 0$.


 Note that both the nonnegativity constraint and the boundedness by $e^{-\gamma}  \def\d{\delta} \def\z{\zeta}  \def\th{\theta x}$ for large $x$
 make this a problem of combining "Laplace rational approximation" with "shape constraints on the inverse", two tasks which do not fit easily together.
\end{comment}





\begin{ex} \label} \def\fr{\frac} \def\im{\item{ex:per} {\bf The Cram\'er Lundberg
process  with Brownian perturbation $\s B(t)$}. The
Laplace
   exponent  is
\begin{eqnarray}} \def\eeq{\end{eqnarray} \label} \def\fr{\frac} \def\im{\item{eq:ps} \kappa(s)=p s + \frac{\k_2 }{2} s^2 + \sum_{k=3}^\infty (-s)^{k} \fr{\l m_k}{k!},\eeq
where $\l m_k$ are the moments of the L\'{e}vy measure, and $\k_2 =
\l m_2 + \s^2$.

Now we have two further unknowns of interest, the
probability of "creeping ruin" and that of "ruin by jump". The respective Laplace transforms   satisfy \begin{eqnarray}} \def\eeq{\end{eqnarray} \label} \def\fr{\frac} \def\im{\item{eq:per}
\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q_d^*(s)=\frac{\s^2/2}{\k'(0)} \f(s)=\frac{\s^2/2}{p} \f(s), \q \Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q_j^*(s)=\fr 1s- \frac{\f(s)}{s} -\frac{\s^2/2}{p} \f(s) \eeq

The Pad\'e $K=0, L=1$ approximations
\begin{comment}of
the transforms of the
Wiener-Hopf factor:
$
\f_{DV}(s)=\frac{2 p}{s (\sigma ^2+ \lambda  m_2)+2 p}$$
yields
\begin{eqnarray*} && \Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q_d^*(s)=\frac{\sigma ^2}{s (\sigma ^2+ \lambda  m_2)+2 p}, \q \Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q_j^*(s)=\frac{\lambda  m_2}{s (\sigma ^2+ \lambda  m_2)+2 p} \; \Leftrightarrow\\&&
\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q_d(x)=\frac{
   \sigma ^2}{\sigma ^2+\lambda  m_2} e^{-\frac{2 p x}{\sigma ^2+\lambda  m_2}}, \q \Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q_j(x)=\frac{
   \lambda  m_2}{\sigma ^2+\lambda  m_2} e^{-\frac{2 p x}{\sigma ^2+\lambda  m_2}}
   \end{eqnarray*}
   \end{comment}
 are  unreasonable,  since they cannot satisfy  the boundary conditions $\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q_d(0)=1, \Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q_j(0)=0.$
To satisfy those as well as the equation
\begin{eqnarray}} \def\eeq{\end{eqnarray} \label} \def\fr{\frac} \def\im{\item{eq:mper} \Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q_d^*(s)=\frac{\s^2/2}{p} \f(s) \Leftrightarrow 1-s \Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q_j^*(s)= (s+ \fr{2 p}{\s^2}) \Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q_d^*(s)\eeq
which follows from \eqref{eq:per},
 we must use at least a  $K=1, L=2$ approximation.

 \begin{thm \label} \def\fr{\frac} \def\im{\item{t:per} Consider the exponential approximation
 \begin{eqnarray} \label} \def\fr{\frac} \def\im{\item{eq:pa}
\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q_d(x)
&=&  \frac{a_d - \m_1}{\m_2-\m_1} e^{-\m_1x} + \frac{ \m_2-a_d }{\m_2-\m_1} e^{-\m_2x} \\
\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q_j(x) &=& \frac{a_j}{ \m_2-\m_1 } [e^{-\m_1 x}   - e^{-\m_2 x}],
\end{eqnarray}
for the "creeping ruin" and  "ruin by jump", which satisfy $\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q_d(0)=1, \Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q_j(0)=0.$ Then, by fitting the first two moments of the aggregate loss $L$,
one is led to the following admissible approximation:
\begin{eqnarray}} \def\eeq{\end{eqnarray} \label} \def\fr{\frac} \def\im{\item{eq:a} a_d=\fr{3 m_2}{m_3}, \; a_j=a_d \fr{ \l m_2}{\s^2 }=\fr{3 \l m_2^2}{\s^2 m_3}
  
   
 
\eeq
and $-\m_1$ and $-\m_2$  the negative roots  of $s^2 + (a_d + a_j + \fr{2 p}{\s^2}) s +  a_d \fr{2 p}{\s^2}=0$ (whose  discriminant is non-negative).

\end{thm}

 {\bf Proof:}  $\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q^*(s),$   $\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q_d^*(s)=\fr{s+ a_d}{s^2 + b_1 s+ b_0}$ are quotients of  monic polynomials (to satisfy $\lim_{s \to \infty} s \Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q^*_d(s)=\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q_d(0)=1$), with three free coefficients, but the second condition in \eqref{eq:mper} imposes one more condition $b_0=\fr{2 p}{\s^2} a_d$ (so that $\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q_j^*(s)= \fr{a_j}{s^2 + b_1 s + a_d \fr{2 p}{\s^2} }$), leaving only two free coefficients.
 Finally,  fitting the first two coefficients around $0$ of
\begin{eqnarray}} \def\eeq{\end{eqnarray} \label} \def\fr{\frac} \def\im{\item{eq:dt} \f(s)=1- s\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q^*(s)=\fr{2 p}{\s^2} \Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q_d^*(s)\approx \fr{2 p}{\s^2} \fr{s+a_d}{s^2 +(a_d + a_j+ \fr{2 p}{\s^2}) s + a_d \fr{2 p}{\s^2} } \eeq
 yields:
 \begin{eqnarray*} \begin{cases}\fr {\s^2}{2 p} (1+\fr{a_j}{ a_d})=\fr{\l m_2 +\s^2}{2 p} \Leftrightarrow
 \fr{a_j}{ a_d}=\fr{\l m_2}{ {\s^2}}\\
 \fr {\s^2}{2 p} \fr{a_j}{ a_d^2}=\fr{\l m_3 }{6 p}\end{cases} \end{eqnarray*}
 with solution \eqref{eq:a}. Then,   assuming
 $\m_1 < \m_2$, Laplace inversion yields \eqref{eq:pa}.



Moreover, both $\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q_j $ and $\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q_d$ are  admissible. Indeed,
this is obvious for $\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q_j $, since its initial value  and its dominant coefficient $\frac{a_j}{ \m_2-\m_1 }$ are non-negative.

The same is true for $\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q_d $; indeed, we may    check  that its dominant coefficient is non-negative, i.e.  that $\m_1 \leq a_d$, by noting that $s^2 + (a_d + a_j + 2p/\sigma^2) s + a_d 2p/\sigma^2$
is negative at $s=-a_d$ (since
$
(a_d)^2 + (a_d + a_j + 2p/\sigma^2) (-a_d) + a_d 2p/\sigma^2 = -a_da_j < 0
$,
with $a_d$ and $a_j$ being positive).
Therefore, $\m_1 < a_d < \m_2$.

\begin{rem}} \def\eeR{\end{rem} It is easy to check  that this approximation is exact for exponential claims. Indeed, in that case  the density transform in the parametrization \eqref{eq:dt} is:
 \begin{eqnarray*}
 &&\f(s)=\fr{p}{\k(s)} = \fr{1}{1+ \fr{\sigma^2}{ 2 p} s +  \fr{\l}{p} s \fr{m_1}{s + m_1^{-1}}}
=  \fr{{s + m_1^{-1}}}{({s + m_1^{-1}})(1+ \fr{\sigma^2}{ 2 p} s) +  \fr{\l}{p} s {m_1}}\\&&=  \fr{s + m_1^{-1}}{\fr{\sigma^2}{ 2 p} s^2 +  s(\fr{\l}{p}  {m_1}+ \fr{\sigma^2}{ 2 p} m_1^{-1}+ 1) + m_1^{-1}}=  \fr{ 2 p} {\sigma^2} \fr{s + m_1^{-1}}{ s^2 +  s(\fr{2 \l {m_1}}    {\sigma^2} + m_1^{-1}+ \fr{ 2 p} {\sigma^2} ) + m_1^{-1} \fr{ 2 p} {\sigma^2} }
\end{eqnarray*} with $a_d=1/m_1, a_j=2 \l m_1/\s^2$, and
it is  easy to check  that this coincides with our approximation \eqref{eq:a},  in the case of exponential claims. \eeR

\begin{rem}} \def\eeR{\end{rem} \label} \def\fr{\frac} \def\im{\item{r:per} Another admissible approximation exact in the  exponential case, but fitting now only one moment is
$a_j=2 \l m_1/\s^2,$
$a_d=1/ \tilde} \def\B{\bar} \def\H{\hat m_1=\fr{2 m_1}{m_2}$ (with similar admissibility proof). \eeR

 \begin{comment}

 $$a_j=\frac{4 p \tilde} \def\B{\bar} \def\H{\hat m_1-\sigma ^2 \th }{2 p
   \tilde} \def\B{\bar} \def\H{\hat m_2}, \; a_d=\frac{\left(4 p \tilde} \def\B{\bar} \def\H{\hat m_1-\sigma ^2 \th
   \right){}^2}{4 p \sigma ^2 \th \tilde} \def\B{\bar} \def\H{\hat m_2},$$

\begin{eqnarray*} && \fr{ \k_2/2- ...}{p+ \k_2 s/2-...}\approx \fr{s+a_d+ a_j}{s^2 +(a_d + a_j+ \fr{2 p}{\s^2}) s + \fr{2 p}{\s^2} a_d} \Leftrightarrow \\&&(\k_2/2)(s^2 +(a_d + a_j+ \fr{2 p}{\s^2}) s + \fr{2 p}{\s^2} a_d) \\&&=(s+a_d+ a_j)(p+ \k_2 s/2) \Leftrightarrow
\begin{cases} \end{cases}\end{eqnarray*}



For the Brownian perturbed model,  we find the Pad\'e $K=1, L=2$ approximation of
the Wiener-Hopf factor is: \begin{eqnarray}\label} \def\fr{\frac} \def\im{\item{PaDV}&&
\f(s)\approx \frac{6 p \left(s m_4+4 m_3\right)}{6 p (m_4 s+4
    m_3)+  \l s(s(3  \k_2 m_4 -4 m_3^2
   )+12
 \k_2 m_3 )}\end{eqnarray}   and the "perturbed De Vylder approximation" is:

\begin{eqnarray*} \Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q^*(s)&\approx& \frac{\l (s(3  \k_2 m_4 -4 m_3^2
   )+12
 \k_2 m_3)}{6 p (m_4
s+4
    m_3)+  \l s(s(3  \k_2 m_4 -4 m_3^2
   )+12
 \k_2 m_3)}\end{eqnarray*}

This corresponds to a combination of two real exponentials with
negative exponents
$$\gamma}  \def\d{\delta} \def\z{\zeta}  \def\th{\theta_{1,2}=\frac{-12 \l \k_2 m_3-6 p m_4 \pm \sqrt{\left(12 \l \k_2
   m_3+6 p m_4\right)^2-96 \l p m_3 \left(3
    \k_2 m_4-4  m_3^2\right)}}{2 \l \left(3  \k_2 m_4-4
   m_3^2\right)}.$$







\end{ex}














\begin{comment}
Just as De Vylder's, Ramsay's  approximation -- see below -- may
 also be obtained as a   Pad\'e $(1,2)$   approximation
 (of the equilibrium distribution).

  A natural generalization
of De Vylder's and Ramsay's work is to look for admissible solutions
of  higher order moments/Pad\'e systems, or
   two-point Pad\'e,

 or  other  rational approximations,
 minimizing for example in least squares
 sense the moment-fitting errors -- see Beylkin and Monz\'on
 \cite{beylkin2005approximation}.
 \end{comment}





\sec{Johnson-Taaffe approximations \label} \def\fr{\frac} \def\im{\item{s:JT}}
This elegant  moment fitting method approximates  by {\bf common order Erlang mixtures approximations}\symbolfootnote[4]{One hint which may explain this choice comes from
discretizing the well-behaved Post-Widder Laplace inversion method \cite{jagerman1978inversion}, which leads
to Erlang mixtures of common order. Note also that pure Erlang distributions
are the { unique extremal points} minimizing the coefficient of variability,
within phase-type distributions representable at order $n$, and that this is a  dense subclass of phase-type distributions.}:
\begin{eqnarray}} \def\eeq{\end{eqnarray} \label} \def\fr{\frac} \def\im{\item{JTrep} \boxed{f(x) \sim \sum _{i=1}^K w_i e_n(\rho} \def\s{\sigma_i,x):=\int_0^\infty e_n(\rho} \def\s{\sigma,x)  \nu_K(d \rho} \def\s{\sigma)}\eeq
where $e_n(\rho} \def\s{\sigma,x)$ is the Erlang density. Finding a discrete measure $\nu_K(d \rho} \def\s{\sigma)$ with positive coefficients $w_i$ requires   solving
a classic De Prony  system
with respect to   {\bf Erlang reduced moments} $$\mu_1=\fr{m_1}{n}, \mu_2=\fr{m_1}{n(n+1)},..., \mu_k=\fr{m_k}{n_{(k)}}, \; \;  n_{(k)}=n(n+1)...(n+k-1)$$ obtained dividing by the moments $n_{(k)}$ of the scaled Erlang density of shape $n$.

 \begin{ex} When $n=1,$ one finds the {\bf factorially reduced moments}
 $$\mu_1=\fr{m_1}{1}, \mu_2=\fr{m_1}{2!}, \mu_3=\fr{m_1}{3!}, ...$$
 which intervene in constructing GHE approximations (which have positive weights, if the moments $m_k$ are far enough from the boundaries of the Stieltjes space of moments $\fr{m_2}{m_1^1}\geq 1, \fr{m_3 m_1}{m_2^1}\geq 1$). \end{ex}

 It turns out that   Erlang reducing leads to "positive moment sequences" $\mu_k$ satisfying $\fr{\mu_2}{\mu_1^1}\geq 1, \fr{\mu_3 \mu_1}{\mu_2^1}\geq 1, ....$ (i.e. having positive Hankel matrices of moments up to any desired degree), for $n$ big enough.   Then, the classic Stieltjes method of creating discrete moment fitting measures with {\bf positive weights} may be applied.


 \begin{thm \label} \def\fr{\frac} \def\im{\item{t:JT} {\bf Johnson-Taaffe three moments
approximation}. Let  $\bff m=(m_1, m_2, m_3)$ denote the first three moments of a nonnegative r.v. Let \begin{eqnarray*}
n^*&=&\lceil \max\{ (m_0 m_2/m_1^2-1)^{-1},
 (m_1 m_3 /m_2^2-1)^{-1} -1 \}\rceil \\&
 =&\lceil \max\{ \fr{1}{\H m_2-1},
 \fr{2 \H m_2-\H m_3}{\H m_3-\H m_2}\}\rceil\end{eqnarray*}
 denote the {\bf explicit "JT index of degree" $3$} giving the smallest "JT sector" -- see  \cite[Fig 3-7]{bobbio2005matching} -- containing our target
 moments.
  Let $\bff \mu =(\mu_i, i=1,...3)$ denote the Erlang reduced moments of order
 $n^*$, and let
 $$ \boxed{b_0= \mu_2 - \mu_1^2, b_1=\mu_3 - \mu_1 \mu_2, b_2= \mu_3 \mu_1 - \mu_2^2}$$
   denote the coefficients of the denominator of the classic second order Pad\'e approximant, and
   let
   $$ \m_{1,2}=\fr{2 b_0}{b_1 \mp \sqrt{b_1^2 - 4 b_0
 b_2}}, \; \w_1= \fr{1- \m_1 \mu_1}{\m_2-\m_1},
 \w_2=1-\w_1
 $$
 denote its roots and (non-negative) partial fractions coefficients.

 Then, the Erlang reduced moments satisfy $\bff \mu \in \{ b_0 > 0, b_2> 0 \}$, and
 $$f(x)= \sum_{i=1}^2 \w_i e_{n^*}(\m_i, x), \;
 $$
 is a nonnegative, decreasing "common order
 Erlang $n^*$" 
  density fitting  our  three moments.

 \end{thm}

Using this result, we can "fix" Ramsay's approximation by sums of two exponentials \cite{ramsay1992practical} to make it work for any valid three moments, at the price of using higher (not minimal) order approximations -- see Theorem \ref{t:Ram}.

\begin{rem}} \def\eeR{\end{rem} {\bf In principle, $2 k+ 1$ moments fitting Johnson-Taaffe Pad\'e
approximations} could be obtained by determining numerically higher order JT
indices ensuring the positivity of  corresponding Hankel
determinants. The Hankel determinants that must be taken into
consideration are: \begin{eqnarray}} \def\eeq{\end{eqnarray} &&\mM_k({\mathcal P})=\{
m_4(m_2 - m_1^2) + m_2 (m_3 m_1 -m_2^2)
+ m_3
(m_2 m_1 -m_3^2)\geq 0,  \label} \def\fr{\frac} \def\im{\item{e:St} \nonumber} \def\bs{\bigskip\\
&&m_5 m_4\left(m_1 m_3-m_2^2\right)+m_4\left(-m_3^3+2 m_2
   m_4 m_3-m_1 m_4^2\right)\geq 0,...,\} \eeq

   The inequalities above  define implicitly the {\bf JT index of degree $5$}.
\begin{comment}
   The JT approach could thus be used to obtain admissible approximations fitting  an arbitrary number of moments
   of either the Pollaczek-Khinchine transform, or
   of the equilibrium distribution.

   An interesting  extension would be to develop admissible two-point Pad\'e
   approximations
   \end{comment}
\eeR


\begin{comment}
Unfortunately, the approximations obtained this way are not of minimal order, and may display considerable oscillation, which brings us to a  second major difficulty, that of not  straying too much from the minimal order. In fact, the order in the set  $JT(n)-JT(n-1)$  fittable by Er (n), but not by Er(n-1) may be always reduced from $2 n$ to at most $n+1$, as shown in \cite{osogami2006closed}.

Recently, an exact formula  for the  minimal order three moments
Coxian representation was obtained in \cite[Fig
3-7]{bobbio2005matching}. It turns out that throughout $JT(n)-JT(n-1)$, the minimal order is either $n+1$ or $n$,  depending on the position of the pair of normalized moments
  $ (\widehat{m}_2, \widehat{m}_3) $ on the
   "Bobbio map"  \cite[Fig 3-7]{bobbio2005matching}.
Once the first Johnson-Taafe sector $JT(n)$ containing  a target
point has been determined, three cases are possible -- see
\cite[Fig 3-7]{bobbio2005matching}: the point can be in the
"spitz" of "Telek's boot of order $n$" in which case an $APH(n)$
fit can be given explicitly, or it can be done above or below the
spitz, in which cases different types of $APH(n+1)$ fits are found
-- see \cite[Fig 5]{bobbio2005matching}.

\end{comment}


\begin{comment}
\vskip.1truecm{\bf Ramsay's approximation.} Ramsay
\cite{ramsay1992practical}  developped a {\bf second order}
 approximation for the classical Cram\'er-Lundberg  model,
  by fitting
  both the correct value $\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q(0)=\fr{1}{1+\th}= \fr{\l m_1}{c} \def\w{\omega} \def\q{\qquad} \def\D{\Delta}\def\Si{\Sigma}$  of the ruin
 probability with $0$ initial reserves (where $\th$  is the  loading factor)
 and  the first
 three moments ${\tilde} \def\B{\bar} \def\H{\hat m_i}, i=1,2,3$ of the stationary
 excess/overshoot/ladder/equilibrium distribution
 \eqref{eq:L} (this was accomplished by a "home made" solution of the de Prony system
 of order $4$.


Ramsay's paper drew attention to the somewhat neglected fact
 that since the  Pollaczek-Khinchine formula only
depends on the nonincreasing  equilibrium density of the claims
$$f_e(x) \sim  \bar F(x) $$
 approximating directly this function
is likely to
lead to better results than approximating the claims
density.

It also provided conditions   ensuring  that the approximation obtained
are valid nonnegative densities (the nonincreasingness of $f_e(x)$
was not ensured).


However, outside the second order Telek nonnegativity domain, \cite{ramsay1992practical}
only proposed using some "experience suggested fix-ups". Nowadays,
one can replace these fix-ups by a systematic three moments Coxian
fitting,   cf Bobbio \& al \cite{bobbio2005matching}, or, for more moments, one may use the Johnson Taaffe approach of using Erlang reduced  moments of order big enough so that the classic Stieltjes constraints applied to the factorially reduced moments hold $\{{\mu_2}-{\mu_1^2} \geq 0, {\mu_3
\mu_1}-{\mu_2^2} \geq 0, ...\}$







\end{comment}

 \begin{thm \label} \def\fr{\frac} \def\im{\item{t:Ram} {\bf Ramsay updated}  Consider the ruin problem for the classic (nonperturbed) Cram\'er-Lundberg  model \eqref{CLmod}.
 \begin{enumerate}}  \def\BI{\begin{itemize}
      \im  Imposing the correct limiting behavior $\lim_{s \to \infty}
      s \Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q^*(s)=\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q(0)=\rho} \def\s{\sigma=1/(1+\th)$  at $\infty$, the
      Pad\'e $(1,2)$  approximations  for the equilibrium density and ruin  transforms are:
    \begin{eqnarray*}&& f_e^*(s)
    \approx \fr N D :=\frac{b_0 + a_1  s }{b_0 + b_1 s +b_2  s^2}, \\&& \Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q^*(s)
    \approx \rho} \def\s{\sigma \fr {D -N} {{D}-\rho} \def\s{\sigma N}: =\rho} \def\s{\sigma \frac{b_2 s +b_1-a_1
     }{b_2 s^2  +\tilde} \def\B{\bar} \def\H{\hat b_1 s+\tilde} \def\B{\bar} \def\H{\hat b_0 },\end{eqnarray*}
    \begin{eqnarray*} && b_0=\tilde} \def\B{\bar} \def\H{\hat \mu_2-\tilde} \def\B{\bar} \def\H{\hat \mu_1^2, \q
    b_1 = \tilde} \def\B{\bar} \def\H{\hat \mu_{3} -\tilde} \def\B{\bar} \def\H{\hat \mu_{2} \tilde} \def\B{\bar} \def\H{\hat \mu_{1}, \q b_2= \tilde} \def\B{\bar} \def\H{\hat \mu_{1} \tilde} \def\B{\bar} \def\H{\hat \mu_{3}-\tilde} \def\B{\bar} \def\H{\hat \mu_{2}^
    \\&& a_1= b_1 -\tilde} \def\B{\bar} \def\H{\hat \mu_1 b_0 = \tilde} \def\B{\bar} \def\H{\hat \mu_3  - 2 \tilde} \def\B{\bar} \def\H{\hat \mu_1 \tilde} \def\B{\bar} \def\H{\hat \mu_2 + \tilde} \def\B{\bar} \def\H{\hat \mu_1^3
     
    
   \\&& \tilde} \def\B{\bar} \def\H{\hat b_1 s = b_1 -\rho} \def\s{\sigma a _1, \q   \tilde} \def\B{\bar} \def\H{\hat b_0 =(1-\rho} \def\s{\sigma)  b_0,
   \end{eqnarray*}
    where $\tilde} \def\B{\bar} \def\H{\hat \mu_{2}$ are {\bf factorially reduced equilibrium
    moments}.

    \im  For the approximate inverse transform $f_e(x)$  to be hyperexponential (in particular, { nonnegative and completely monotone}), it is sufficient that $b_0 >0, b_2>0$.

 
 

  
  





   \im  Higher order moments fitting is possible by Erlang approximations of order $n^*$ big enough, where for three moments $n^*=\lceil \max\{\fr {m_1^2}{m_2-m_1^2},\fr
 {m_2^2}{m_1 m_3 -m_2^2} -1 \}\rceil $.
 \end{enumerate}}   \def\EI{\end{itemize}
   \end{thm}

   \begin{comment}

   {\bf Proof:} The equivalence of the two Laplace transform formulas
   follows from the identity
   \begin{eqnarray}} \def\eeq{\end{eqnarray} \label} \def\fr{\frac} \def\im{\item{eq:id} && s \Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q^*(s)= \rho} \def\s{\sigma \fr{1- f_e^*(s)}{1- \rho} \def\s{\sigma f_e^*(s)}=\rho} \def\s{\sigma \frac{b_2  s^2 + (b_1 -a_1)  s }{b_2  s^2 + (b_1 - \rho} \def\s{\sigma a_1)  s + b_0(1-\rho} \def\s{\sigma)} \eeq

   The rest is a slight extension of well known  and  easily checked results summarized
   in Theorem \ref{t:fea2} (for example $f_e^*(s) (b_0 + b_1 s +b_2  s^2)-(b_0 + a_1  s )= O(s^4)$, etc).
   \end{comment}




\sec{Comparison between the JT indices of  Ramsay and Beekman-Bowers
type approximations \label} \def\fr{\frac} \def\im{\item{s:comp}}
 In this section we present a simple comparison between these two approaches, based on comparing their  Johnson-Taaffe indices.

 The  moments of the aggregate loss
$L$  satisfy
\begin{eqnarray}\label} \def\fr{\frac} \def\im{\item{sp2} &&\widehat{l}_2=\fr{{l}_2}{{{l}_1^2}}=2+ \th
\frac{\tilde} \def\B{\bar} \def\H{\hat{m}_2}{\tilde} \def\B{\bar} \def\H{\hat{m}_1^2}=2+ \th \widehat{\tilde} \def\B{\bar} \def\H{\hat m}_2 \geq  2 + \th >2,\\
&&\widehat{l}_3=\fr{{l}_3}{{{l}_1 {l}_2}}=3+  \frac{ \frac{\tilde} \def\B{\bar} \def\H{\hat{m}_3}{\th}
  +  \frac{\tilde} \def\B{\bar} \def\H{\hat{m}_1 \tilde} \def\B{\bar} \def\H{\hat{m}_2}{ \th^2}
  }{
  \frac{\tilde} \def\B{\bar} \def\H{\hat{m}_2 \tilde} \def\B{\bar} \def\H{\hat{m}_1}{\th^2}+   2(\frac{\tilde} \def\B{\bar} \def\H{\hat{m}_1}{\th})^3}=
  3+  \frac{ \th \widehat{\tilde} \def\B{\bar} \def\H{\hat{m}}_2(1+ {\th} {\widehat{\tilde} \def\B{\bar} \def\H{\hat{m}}_3)}
  }{2+ \th \widehat{\tilde} \def\B{\bar} \def\H{\hat{m}}_2} \nonumber} \def\bs{\bigskip \end{eqnarray}




 The second moment   satisfies the necessary  inequality  ${\widehat{l}_2}
  \geq \fr{n+1}{n}$ \cite{bobbio2005matching} already with $n=1.$
 We only need to investigate  the  necessary inequality  for the third moment,
 which is:
 \begin{eqnarray*} n &\geq& \boxed{J(L):=\frac{2 \widehat{l}_2-\widehat{l}_3}
 { \widehat{l}_3-\widehat{l}_2}}=\frac{\left(\widehat{\tilde} \def\B{\bar} \def\H{\hat{m}}_2 \widehat{\tilde} \def\B{\bar} \def\H{\hat{m}}_3-2
   \widehat{\tilde} \def\B{\bar} \def\H{\hat{m}}_2^2\right) \theta ^2-4
   \widehat{\tilde} \def\B{\bar} \def\H{\hat{m}}_2 \theta -2}{\theta ^2
   \left(\widehat{\tilde} \def\B{\bar} \def\H{\hat{m}}_2^2-\widehat{\tilde} \def\B{\bar} \def\H{\hat{m}}_2
   \widehat{\tilde} \def\B{\bar} \def\H{\hat{m}}_3\right)-2} \Leftrightarrow \\&  \fr{n+2}{n+1}&
\leq \boxed{\n(L):=\fr{  \widehat{l}_3}{  \widehat{l}_2}
} =\frac{\theta  \widehat{\tilde} \def\B{\bar} \def\H{\hat{m}}_2 \left(\theta
\widehat{\tilde} \def\B{\bar} \def\H{\hat{m}}_3+4\right)+6}{\left(\theta  \widehat{\tilde} \def\B{\bar} \def\H{\hat{m}}_2+2\right){}^2}
\end{eqnarray*}

Here, $J(L)=\frac{2 \widehat{l}_2-\widehat{l}_3}
 { \widehat{l}_3-\widehat{l}_2}$ is a "partial JT index", based on the third moment
 admissibility condition, and the "normalized moment" $\n(L)=\fr{  \widehat{l}_3}{  \widehat{l}_2}$ is a monotone transformation (since $J=\fr{2-\n}{\n-1}$ is a decreasing function in the relevant range $\n \in (1,\infty)$), which has been already used in the literature.

 A rough indication of the performance of the Beekman-Bowers and Ramsay methods will be obtained now by checking which
 of $\n(L)$,  $\n(L_i)$ is higher.





 \begin{thm \label} \def\fr{\frac} \def\im{\item{t:new} a) The partial  J index of the aggregate loss is strictly  smaller than $J(L_i)=\fr{2\widehat{\tilde} \def\B{\bar} \def\H{\hat m}_2-\widehat{\tilde} \def\B{\bar} \def\H{\hat m}_3} {\widehat{\tilde} \def\B{\bar} \def\H{\hat m}_3-\widehat{\tilde} \def\B{\bar} \def\H{\hat m}_2}$  iff \begin{eqnarray*} \theta \widehat{\tilde} \def\B{\bar} \def\H{\hat m}_2 <\frac{3/2 \widehat{\tilde} \def\B{\bar} \def\H{\hat m}_2-
   \widehat{\tilde} \def\B{\bar} \def\H{\hat m}_3}{  \widehat{\tilde} \def\B{\bar} \def\H{\hat m}_3- \widehat{\tilde} \def\B{\bar} \def\H{\hat m}_2} =\frac{3/2 -\n(L_i)}{  \n(L_i)-1} \end{eqnarray*}

 In particular, if $\n(L_i)=\fr{\widehat{\tilde} \def\B{\bar} \def\H{\hat m}_3} {\widehat{\tilde} \def\B{\bar} \def\H{\hat m}_2} \geq  \fr 32,$
 in which case $\n(L) <  \fr 32$ holds as well,
 a  two terms exponential mixture distribution matching  the
first three moments exists only for the equilibrium ladders $L_i$.


b) More generally, let $n$ denote the unique integer such that  $\fr {n+2}{n+1} \leq \n(L_i) \leq  \fr {n+1}{n}, n \geq 2,$ is satisfied. If furthermore $ x_1 < \th \widehat{\tilde} \def\B{\bar} \def\H{\hat m}_2 <x_2,$ where $x_{1,2}=\frac{2 \mp \sqrt{2} \sqrt{n^2 +n-a(n^2-1)
   }}{(n +1)a-(n+2)}$, then
a  $n$'th order Erlang mixture distribution matching  the
first three moments exists only for the equilibrium ladders $L_i$.





c) For any $\th,$ there exists a mixture of Erlang $n$ distribution   matching the
first three moments of $L$ if
 $$n \geq \frac{\widehat{\tilde} \def\B{\bar} \def\H{\hat m}_3}{\widehat{\tilde} \def\B{\bar} \def\H{\hat m}_3-\widehat{\tilde} \def\B{\bar} \def\H{\hat m}_2}.$$
 \end{thm}

 {\bf Proof:} a) Consider the unimodal function
 $$\n(a,x)=\frac{\widehat{l}_3}{ \widehat{l}_2} =\fr{a x^2 + 4 x +6}{(x+2)^2},$$
 where $a=\frac{\widehat{\tilde} \def\B{\bar} \def\H{\hat m}_3}{\widehat{\tilde} \def\B{\bar} \def\H{\hat m}_2}, x= \th \widehat{\tilde} \def\B{\bar} \def\H{\hat m}_2$,
 which takes values  $$\fr{3 a -2}{2a -1} \leq  \n(a,x) \leq \max[\n(a,\infty)=a,\n(a,0)=3/2]$$
 (the lower minimal value is achieved for $x^*=\fr{1}{a-1} \Leftrightarrow \th^*=\frac{1}{ \widehat{\tilde} \def\B{\bar} \def\H{\hat m}_3-\widehat{\tilde} \def\B{\bar} \def\H{\hat m}_2}$).

b) Fitting $L$ is possible when
\begin{eqnarray*} &&\fr{a x^2 + 4 x +6}{(x+2)^2} > \fr{n+2}{n+1} \Leftrightarrow ((n+1) a -(n+2))x^2 - 4 x + 2(n-1) >0. \end{eqnarray*}

In the prescribed range of $a,$ the discriminant $\D=2(n^2 +n-a(n^2-1))$ is always positive, and the inequality holds  when $x < x_1$ or when $x > x_2,$ where $x_{1,2}=\frac{2 \mp \sqrt{2} \sqrt{n^2 +n-a(n^2-1)
   }}{(n +1)a-(n+2)}$.

 c) Follows  by minimizing $\n(\th)$
  in $\th,$ which yields $\th^*=\frac{1}{ \widehat{\tilde} \def\B{\bar} \def\H{\hat m}_3-\widehat{\tilde} \def\B{\bar} \def\H{\hat m}_2}$ and
 $\n(\th^*)=\frac{3 \widehat{\tilde} \def\B{\bar} \def\H{\hat m}_3-2 \widehat{\tilde} \def\B{\bar} \def\H{\hat m}_2}{2 \widehat{\tilde} \def\B{\bar} \def\H{\hat m}_3-\widehat{\tilde} \def\B{\bar} \def\H{\hat m}_2}, J(\th^*)=\frac{\widehat{\tilde} \def\B{\bar} \def\H{\hat m}_3}{\widehat{\tilde} \def\B{\bar} \def\H{\hat m}_3-\widehat{\tilde} \def\B{\bar} \def\H{\hat m}_2}$.

\begin{ex} \label} \def\fr{\frac} \def\im{\item{ex:uni} Consider   a three moments fitting example:  the $U[0,1]$ rv,
  with moments $m_1 =\fr 12, m_2 =\fr 13, m_3 =\fr 14, m_4 =\fr 15$.  The Johnson-Taaffe index is $7$, and so this procedure yields $14$ phases, while the Bobbio, Horvath and Telek \cite{bobbio2005matching} and He-Zhang methods
   \cite{he2006spectral}  yield only $9$ phases.

For the equilibrium distribution, the moments are $\tilde} \def\B{\bar} \def\H{\hat m_1  =\fr 13, \tilde} \def\B{\bar} \def\H{\hat m_2 =\fr 16, \tilde} \def\B{\bar} \def\H{\hat m_3 =\fr {1}{10},$ the normalized moments are
$$\widehat{\tilde} \def\B{\bar} \def\H{\hat m}_2= \fr{\tilde} \def\B{\bar} \def\H{\hat m_2}{\tilde} \def\B{\bar} \def\H{\hat m_1^2}=\fr 32, \widehat{\tilde} \def\B{\bar} \def\H{\hat m}_3= \fr{\tilde} \def\B{\bar} \def\H{\hat m_3 }{\tilde} \def\B{\bar} \def\H{\hat m_2 \tilde} \def\B{\bar} \def\H{\hat m_1}= \fr 95$$ and the  Johnson-Taaffe index is $4$.  The Bobbio, Horvath and Telek  method
    yields the minimal order of $3$ phases.

Finally consider  the JT approach for a "Beekman-Bowers"  approximation.
 Computing the partial JT index, we find
$$ \frac{4 \left(9 \theta ^2+30
   \theta +10\right)}{9 \theta
   ^2+40}.$$  By Theorem \ref{t:new} a), this is strictly less than the $4$ required for the Ramsay approach  when $\th \leq 1$. Taking into account
   the integer part, we find out that $\th \leq 1.5$ still yields a J index
   less or equal  to $4$.
   \end{ex}

\begin{figure}
\begin{center}
\includegraphics[scale=0.7]{uni.pdf}
\end{center}
\caption{Ruin probabilities for  uniform claims.
Red line: original calculated by numerical inverse Laplace;
green line: by 3 moments \cite{johnson1989matching} on the uniform itself;
blue line: by 3 moments \cite{bobbio2005matching} on the uniform itself;
yellow line: by 3 moments \cite{johnson1989matching} on the equilibrium of the uniform;
black line: by 3 moments \cite{bobbio2005matching} on the equilibrium of the uniform.\label{fig:5appr}}
\end{figure}




\begin{comment}
   \begin{ex}  It is possible  that  a second order  Coxian distribution matching the
first three moments of the aggregate loss $L$ exists as well, if

If $\n(L_i)=\fr{\widehat{\tilde} \def\B{\bar} \def\H{\hat m}_3} {\widehat{\tilde} \def\B{\bar} \def\H{\hat m}_2} < \fr 43,$ but $\th \leq \frac{\sqrt{6} \sqrt{2 -\fr{\widehat{\tilde} \def\B{\bar} \def\H{\hat m}_3} {\widehat{\tilde} \def\B{\bar} \def\H{\hat m}_2}} -2}{4 \widehat{\tilde} \def\B{\bar} \def\H{\hat m}_2-3 \widehat{\tilde} \def\B{\bar} \def\H{\hat m}_3},$ then
 there exists a  third order  Coxian distribution matching exactly the
first three moments of $L,$ but not those of $L_i$.


  \end{ex}
  \end{comment}
\begin{comment}
   \begin{rem}} \def\eeR{\end{rem}  To prevent oscillation (see Example \ref{ex:uni}), it might be  useful to include other information than the moments,  for example the "Markov parameters" (the power
series coefficients of the Laplace transform at infinity), giving rise to two point Pad\'e approximations -- see Section \ref{s:prop}.

In our problem, it makes sense to utilize for the "Beekman-Bowers type approach" at least the first  expansion coefficient of the
Laplace transform $\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q^*(s)$ at $\infty$; this is especially true in the perturbed case, when this coefficient distinguishes "smooth crossing" of "crossing by jump" -- see Example \ref{ex:per}.

Using only the first coefficient amounts, by the Pollaczek Khinchine  formula \ref{PK}, to ensuring the
admissibility  of the stationary excess density $f_e(x)$, while
using the second   coefficient as well amounts to ensuring the
admissibility  of  $f_e(x)$ under the additional constraint $f_e(0)=m_1^{-1}$ -- see \eqref{eq:twoc}.
\eeR


\begin{rem}} \def\eeR{\end{rem}
 Another classic approach which prevents oscillation is by selecting special nodes (for example   the Gauss-Legendre nodes, or the Gauss-Jacobi nodes, which may be used after a rescaling to $[0,1]$ -- see Dufresne \cite{dufresne2007fitting}). 
  \eeR


\end{comment}



\sec{Admissible two point Pad{\'e} approximation
\label} \def\fr{\frac} \def\im{\item{s:prop}}

We consider now briefly  including further information beyond moments,  for example
 by using two-point Pad\'e approximations (of special interest in the case of Levy processes with
 infinite variation paths -- see Example ). Ensuring admissibility in this case is however an open problem.

 We may
formally expand the ruin transform at infinity as well:
\begin{eqnarray}\label} \def\fr{\frac} \def\im{\item{rexp} \Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q^*(s)=\int_0^\infty e^{-s x} \sum_{k=0}^\infty
\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q^{(k)} (0) \frac{ x^k}{k!}=\sum_{k=0}^\infty \Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q^{(k)} (0)
s^{-k-1}\end{eqnarray} Here,
$$\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q (0)=\begin{cases}1 & \text{if } \s >0\\
\frac{\n_1}{c} \def\w{\omega} \def\q{\qquad} \def\D{\Delta}\def\Si{\Sigma}=\rho} \def\s{\sigma &\text{if } \s =0\end{cases}$$ is well known,
and
 the derivatives at $0$: \begin{eqnarray}\label} \def\fr{\frac} \def\im{\item{der}
 &&\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q' (0)=-\frac{\n_0}{c} \def\w{\omega} \def\q{\qquad} \def\D{\Delta}\def\Si{\Sigma}(1-\frac{\n_1}{c} \def\w{\omega} \def\q{\qquad} \def\D{\Delta}\def\Si{\Sigma})=- \frac{p \l}{c} \def\w{\omega} \def\q{\qquad} \def\D{\Delta}\def\Si{\Sigma^2}=- \frac{\rho} \def\s{\sigma(1-\rho} \def\s{\sigma) }{m_1},\\
&&\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q''(0) = -\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q' (0) ( f(0)-\frac{\n_0}{c} \def\w{\omega} \def\q{\qquad} \def\D{\Delta}\def\Si{\Sigma} ),
\nonumber} \def\bs{\bigskip\\&&\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q^{(3)}(0) = -\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q' (0) \left(f'(0) +2 \frac{\n_0}{c} \def\w{\omega} \def\q{\qquad} \def\D{\Delta}\def\Si{\Sigma}
f(0)  -(\frac{\n_0}{c} \def\w{\omega} \def\q{\qquad} \def\D{\Delta}\def\Si{\Sigma} )^2\right),\nonumber} \def\bs{\bigskip\\&& \Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q^{(4)}(0) = -\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q'
(0) \left(f''(0)+2 \frac{\n_0}{c} \def\w{\omega} \def\q{\qquad} \def\D{\Delta}\def\Si{\Sigma} f'(0)-\frac{\n_0}{c} \def\w{\omega} \def\q{\qquad} \def\D{\Delta}\def\Si{\Sigma} f(0)^2 +3
(\frac{\n_0}{c} \def\w{\omega} \def\q{\qquad} \def\D{\Delta}\def\Si{\Sigma}) ^2 f(0) -(\frac{\n_0}{c} \def\w{\omega} \def\q{\qquad} \def\D{\Delta}\def\Si{\Sigma})
   ^3\right) \nonumber} \def\bs{\bigskip\end{eqnarray} may be obtained recursively, by
differentiating the integro-differential equation
for
$\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q(x)$.
\begin{comment}
Alternatively, they may be obtained by computing
recursively the limits: \begin{eqnarray*} &&d_k=\lim_{s \to \infty} s f_k, f_k =
s(f_{k-1}- d_{k-1}), k=1, ..., \\&& f_0=\frac{p}{c} \def\w{\omega} \def\q{\qquad} \def\D{\Delta}\def\Si{\Sigma} -
\frac{p}{\m-\l \bar{F}^*(s)}, d_0=0,\\&&f_1=-\frac{\l p}{c} \def\w{\omega} \def\q{\qquad} \def\D{\Delta}\def\Si{\Sigma}
\frac{1-f^*(s)}{\m-\l \bar{F}^*(s)}, d_1=-\frac{\l p}{c} \def\w{\omega} \def\q{\qquad} \def\D{\Delta}\def\Si{\Sigma^2}\\&&
f_2=s \frac{\l p}{\m^2}\left} \def\ri{\right} \def\i{\infty(1- \frac{1-f^*(s)}{1-\frac{\l}{c} \def\w{\omega} \def\q{\qquad} \def\D{\Delta}\def\Si{\Sigma}
\bar{F}^*(s)}\ri)=- d_1 \left} \def\ri{\right} \def\i{\infty(\frac{ s
f^*(s)-\frac{\l}{c} \def\w{\omega} \def\q{\qquad} \def\D{\Delta}\def\Si{\Sigma}(1-f^*(s))}{1-\frac{\l}{c} \def\w{\omega} \def\q{\qquad} \def\D{\Delta}\def\Si{\Sigma}
\bar{F}^*(s)}\ri),\\&& d_2=- d_1 \left} \def\ri{\right} \def\i{\infty(  f(0)-\frac{\l}{c} \def\w{\omega} \def\q{\qquad} \def\D{\Delta}\def\Si{\Sigma}\ri),
 \end{eqnarray*}
if $f'(x) \in L_1$.

 This opens the way for developing two point Pad{\'e}
approximations, which will approximate $\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q^*(s)$ both at $0$ and
at $\infty$.

\begin{rem}} \def\eeR{\end{rem} Note that while the classic Pad\'e  approximation around $0$
may be obtained directly in many CAS systems (in Mathematica, for
example, we may use  the command PadeApproximant), for two point
Pad\'e, we need to solve  the Pad\'e equations, taking into
account the additional  constraints. \eeR
\end{comment}

\begin{thm \label} \def\fr{\frac} \def\im{\item{t:Ram2} {\bf Two-point Pad\'e-Ramsay approximation}
      Imposing both the correct limiting behavior $\lim_{s \to \infty}
      s \Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q^*(s)=\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q(0)=\rho} \def\s{\sigma=1/(1+\th)$  at $\infty$, and  the first
      derivative at $\infty$ $\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q' (0)=- \frac{\rho} \def\s{\sigma (1-\rho} \def\s{\sigma)}{m_1}$ leads to the two-point Pad\'e $(1,2)$  approximation:
      \begin{eqnarray}} \def\eeq{\end{eqnarray} \label} \def\fr{\frac} \def\im{\item{eq:twoc} && f_e^*(s)
    \approx \frac{b_0 + a_1  s }{b_0 + b_1 s +b_2  s^2}, \; \; \Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q^*(s)
    \approx \rho} \def\s{\sigma \frac{ b_2 s + b_1   -
    {a_1}  }{\tilde} \def\B{\bar} \def\H{\hat b_0 +\tilde} \def\B{\bar} \def\H{\hat b_1 s + b_2 s^2},\eeq
     where $a_1= \fr{b_2}{m_1}  $ and the coefficients $b_i, \tilde} \def\B{\bar} \def\H{\hat b_i,$  obtained fitting the first two equilibrium/aggregate loss moments, are:
    \begin{eqnarray*}&& b_2= (2 m_1 m_3-3  m_2^2)/6, \; b_1=
   (m_3-3 m_1 m_2)/3, \; b_0= m_2 -
    2  m_1^2, \\ && \tilde} \def\B{\bar} \def\H{\hat b_1 =b_1-\rho} \def\s{\sigma a_1, \q  \tilde} \def\B{\bar} \def\H{\hat b_0 =(1-\rho} \def\s{\sigma) b_0.
   \end{eqnarray*}



 \begin{comment}
 b) For the approximate inverse transform $f_e(x)$  to be hyperexponential (in particular, { nonnegative and completely monotone}), it is sufficient that $b_0 >0, b_2>0$.

 c)   The admissibility domain for $f_e(x)$ is:
   \begin{eqnarray*} &&\left(1<\widehat m_2<2\land  \widehat m_3>\frac{3 \left(3
   \widehat m_2-4\right)}{\widehat m_2}+3 \sqrt{2}
   \sqrt{-\frac{\widehat m_2^3-6 \widehat m_2^2+12
   \widehat m_2-8}{\widehat m_2^2}}\right)\lor \\&&\left(2<\widehat m_2<4\land \widehat m_2<\widehat m_3<\frac{3
   \widehat m_2}{4-\widehat m_2}\right)\lor \left(\widehat m_2\geq 4\land \widehat m_3>\widehat m_2\right)\end{eqnarray*}
   or
   \begin{eqnarray*} && \left(\sqrt{3}<\widehat m_2<2\land \widehat m_2<\widehat m_3<\frac{3 \left(3
   \widehat m_2-4\right)}{\widehat m_2}-3 \sqrt{2}
   \sqrt{-\frac{\widehat m_2^3-6 \widehat m_2^2+12
   \widehat m_2-8}{\widehat m_2^2}}\right)\\&&\lor \left(2<\widehat m_2<8-4
   \sqrt{2}\land \widehat m_3>3 \sqrt{2}
   \sqrt{\frac{\widehat m_2^5-6 \widehat m_2^4+12 \widehat m_2^3-8
   \widehat m_2^2}{\left(\widehat m_2^2-16
   \widehat m_2+32\right){}^2}}-\frac{3 \left(5 \widehat m_2^2-12
   \widehat m_2\right)}{\widehat m_2^2-16 \widehat m_2+32}\right),\end{eqnarray*}
   where $\widehat m_k=\frac{{m}_{k}}{{m}_{k-1} m_1}, k=2,3,...$ denote "normalized moments".
   \end{comment}
   \end{thm}




    {\bf Proof:}  The theorem follows from the identity
   \begin{eqnarray}} \def\eeq{\end{eqnarray} \label} \def\fr{\frac} \def\im{\item{eq:id} && s \Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q^*(s)= \rho} \def\s{\sigma \fr{1- f_e^*(s)}{1- \rho} \def\s{\sigma f_e^*(s)}=\rho} \def\s{\sigma \frac{b_2  s^2 + (b_1 -a_1)  s }{b_2  s^2 + (b_1 - \rho} \def\s{\sigma a_1)  s + b_0(1-\rho} \def\s{\sigma)} \eeq
   by taking into account that we are looking for an approximation of the form
      \begin{eqnarray*} \Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q^*(s) \approx \rho} \def\s{\sigma \frac{ b_2 s +  b_1   -  \fr{ b_2} {m_1} }{  b_2 s^2+ \tilde} \def\B{\bar} \def\H{\hat b_1 s+  \tilde} \def\B{\bar} \def\H{\hat b_0 }:=\rho} \def\s{\sigma \fr{\tilde} \def\B{\bar} \def\H{\hat N}{\tilde} \def\B{\bar} \def\H{\hat D},\end{eqnarray*}
      which will satisfy $\lim_{s \to \infty} s
   \frac{\tilde} \def\B{\bar} \def\H{\hat N} {\tilde} \def\B{\bar} \def\H{\hat D}=1, \lim_{s \to \infty} s(s
   \frac{\tilde} \def\B{\bar} \def\H{\hat N} {\tilde} \def\B{\bar} \def\H{\hat D}-1) =\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q'(0)/\rho} \def\s{\sigma= -\frac{1-\rho} \def\s{\sigma}{m_1}.
   $
\begin{comment}
   b) By Theorem \ref{t:fea2} b), nonnegativity of $f_e(x)$ holds iff the discriminant \begin{eqnarray*} &&\left(m_3-3 m_1 m_2\right){}^2+ 6 \left(2
   m_1^2- m_2\right) \left( 2 m_1 m_3-3
   m_2^2\right
  
   
 
  \end{eqnarray*}
    of the denominator of the Laplace transform is nonnegative, and if the root of the numerator is larger than the smallest root of the
   denominator. A tedious calculation (doable with Mathematica) yields the result.



\begin{Qu}  Find high-order admissible   two-point Pad\'e
approximations of ruin probabilities. \end{Qu}
\end{comment}


\section{Numerical results\label{sec:num}}

\subsection{Without perturbation}

In the following we illustrate the application of the approximations
provided in Theorem~\ref{t:Ram}~and~\ref{t:Ram2}.

\subsubsection{Mixed exponential claim distribution}

With mixed exponential claim distribution
\[
f(x)=
\frac{315 e^{-5 x}}{128} + \frac{7 e^{-4 x}}{8} + \frac{27 e^{-3 x}}{64} +
  \frac{3 e^{-2 x}}{16} + \frac{7 e^{-x}}{128}
\]
and $\lambda=1, c=2/5$ the ruin probability is
\[
\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q(x)=
\frac{245e^{-9x/2}}{32768}
+\frac{135e^{-7x/2}}{8192}
+\frac{567e^{-5x/2}}{16384}
+\frac{735e^{-3x/2}}{8192}
+\frac{19845e^{-x/2}}{32768}
\]

We approximated the ruin probabilities by the  Renyi and De Vylder's first order exponential
approximation given in \eqref{e:Re} and \eqref{e:DVL}, and by the approximations provided in
Theorem~\ref{t:Ram}~and~\ref{t:Ram2}.  In Figure~\ref{fig:mixexpruin} we
show the ruin probabilities with the different approximations and in
Figure~\ref{fig:mixexperror} the relative error with respect to the exact
solution.  As expected, the approximation in Theorem~\ref{t:Ram2}  is better near $x=0,$ and that may be exploited to obtain a better approximation by switching between the formulas in Theorem~\ref{t:Ram2}, Theorem~\ref{t:Ram} when they cross, starting with the first.

\begin{figure}
\begin{center}
\includegraphics[scale=0.6]{mixedexpruin.pdf}
\end{center}
\caption{Ruin probabilities with mixed exponential claims.\label{fig:mixexpruin}}
\end{figure}

\begin{figure}
\begin{center}
\includegraphics[scale=0.6]{mixedexp.pdf}
\end{center}
\caption{Absolute relative error of ruin probabilities with mixed exponential claims.\label{fig:mixexperror}}
\end{figure}

   \subsubsection{Gamma distribution with $\alpha=0.01, \beta=100$}

Let us consider now the Gamma claim distribution
\[
\frac{e^{-\frac{x}{\beta }} x^{\alpha -1} \beta ^{-\alpha }}{\Gamma (\alpha )}
\]
with parameters $\alpha=0.01, \beta=100$ and with claim arrival intensity
$\l=1$ and loading factor $\th=0.1$, which appeared frequently in the
literature (\cite{ramsay1992practical}, \cite{grandell2000simple}).  As for the previous example, we approximated the
ruin probabilities by the formulas given in Renyi,
\eqref{e:DVL} and Theorem~\ref{t:Ram}~and~\ref{t:Ram2}.  Ruin
probabilities themselves are given in Table~\ref{tab:gamma1ruin} while in
Table~\ref{tab:gamma1error} we provide absolute relative errors.


\begin{table}
\begin{center}
\footnotesize
\begin{tabular}{llllll}
\hline
 $x$ & exact~$\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q(x)$ & Renyi & DeVylder &
  Theorem~\ref{t:Ram} & Theorem~\ref{t:Ram2} \\
 \hline
0. & 0.909091 & 0.909091 & 0.882867 & 0.909091 & 0.909091 \\
 300. & 0.521143 & 0.529743 & 0.522539 & 0.521107 & 0.522526 \\
 600. & 0.308668 & 0.30869 & 0.309273 & 0.308713 & 0.309268 \\
 900. & 0.182866 & 0.179879 & 0.183048 & 0.182888 & 0.183047 \\
 1200. & 0.108338 & 0.104818 & 0.10834 & 0.108347 & 0.10834 \\
 1500. & 0.0641841 & 0.0610794 & 0.0641226 & 0.0641869 & 0.0641233 \\
 1800. & 0.0380254 & 0.035592 & 0.037952 & 0.0380257 & 0.0379527 \\
 2100. & 0.0225279 & 0.0207401 & 0.0224625 & 0.0225272 & 0.0224631 \\
 2400. & 0.0133465 & 0.0120856 & 0.0132948 & 0.0133456 & 0.0132953 \\
 2700. & 0.00790706 & 0.00704247 & 0.00786872 & 0.0079062 & 0.00786908 \\
 3000. & 0.00468448 & 0.00410377 & 0.00465722 & 0.0046838 & 0.00465748\\
\hline
\end{tabular}
\caption{Ruin probabilities with Gamma claims ($\alpha=0.01, \beta=100$)
  and its approximations\label{tab:gamma1ruin}}
\end{center}
\end{table}

\begin{table}
\begin{center}
\footnotesize
\begin{tabular}{lllll}
\hline
 $x$ & Renyi & DeVylder &
  Theorem~\ref{t:Ram} & Theorem~\ref{t:Ram2} \\
\hline
 0 & 0 & 0.0288462 & 0 & 0 \\
 300 & 0.0165011 & 0.00267814 & 0.0000688428 & 0.00265325 \\
 600 & 0.0000714085 & 0.0019599 & 0.000146799 & 0.00194387 \\
 900 & 0.0163373 & 0.000993297 & 0.000119325 & 0.000986105 \\
 1200 & 0.0324864 & 0.0000177571 & 0.0000819863 & 0.000019394 \\
 1500 & 0.0483709 & 0.000957418 & 0.0000440626 & 0.00094697 \\
 1800 & 0.0639946 & 0.00193166 & $6.12674\times 10^{-6}$ & 0.00191241 \\
 2100 & 0.0793618 & 0.00290493 & 0.000031798 & 0.00287691 \\
 2400 & 0.0944767 & 0.00387725 & 0.0000697035 & 0.00384047 \\
 2700 & 0.109343 & 0.00484863 & 0.000107635 & 0.00480311 \\
 3000 & 0.123966 & 0.00581907 & 0.000145554 & 0.00576482\\
\hline
\end{tabular}
\caption{Relative error of approximate ruin probabilities with Gamma claims
  ($\alpha=0.01,\beta=100$) \label{tab:gamma1error}}
\end{center}
\end{table}

\subsubsection{Gamma distribution with $\alpha=2.5, \beta=1$}

Next we consider Gamma distributed claims with $\alpha=2.5, \beta=1$ and
$\l=2/5, c=\frac{4}{5} \left(-1+4 \sqrt{2}\right)$.  This case is
interesting as it has been shown in \cite{avram2011moments} that direct,
moment based Pad{\'e} approximation of the claim distribution does not
result in valid distributions. The approximations presented in this paper
leads instead to valid ruin probabilities.  Ruin probabilities themselves
are given in Table~\ref{tab:gamma2ruin} while in
Table~\ref{tab:gamma2error} we provide absolute relative errors.


\begin{table}
\begin{center}
\footnotesize
\begin{tabular}{llllll}
\hline
 $x$ & exact~$\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q(x)$ & Renyi & DeVylder &
  Theorem~\ref{t:Ram} & Theorem~\ref{t:Ram2} \\
 \hline
 0. & 0.268422 & 0.268422 & 0.299749 & 0.268422 & 0.268422 \\
 0.5 & 0.22854 & 0.217791 & 0.237348 & 0.22894 & 0.228126 \\
 1. & 0.189678 & 0.176711 & 0.187938 & 0.189655 & 0.189069 \\
 1.5 & 0.154441 & 0.143379 & 0.148813 & 0.154172 & 0.154016 \\
 2. & 0.124037 & 0.116334 & 0.117834 & 0.123743 & 0.123926 \\
 2.5 & 0.0986589 & 0.0943911 & 0.0933036 & 0.0984496 & 0.0988216 \\
 3. & 0.0779451 & 0.0765868 & 0.07388 & 0.0778418 & 0.0782763 \\
 3.5 & 0.0612929 & 0.0621407 & 0.0584999 & 0.0612758 & 0.0616894 \\
 4. & 0.0480435 & 0.0504196 & 0.0463215 & 0.0480817 & 0.04843 \\
 4.5 & 0.0375759 & 0.0409093 & 0.0366785 & 0.0376414 & 0.0379079 \\
 5. & 0.0293456 & 0.0331929 & 0.0290429 & 0.0294185 & 0.0296037\\
\hline
\end{tabular}
\caption{Ruin probabilities with Gamma claims ($\alpha=2.5, \beta=1$)
  and its approximations\label{tab:gamma2ruin}}
\end{center}
\end{table}

\begin{table}
\begin{center}
\footnotesize
\begin{tabular}{lllll}
\hline
 $x$ & Renyi & DeVylder &
  Theorem~\ref{t:Ram} & Theorem~\ref{t:Ram2} \\
\hline
0. & 0. & 0.116709 & 0. & 0. \\
 0.5 & 0.0470337 & 0.0385392 & 0.00175079 & 0.0018135 \\
 1. & 0.0683674 & 0.00917815 & 0.000125885 & 0.00321361 \\
 1.5 & 0.0716262 & 0.0364389 & 0.00174365 & 0.00275474 \\
 2. & 0.062096 & 0.0500072 & 0.00236637 & 0.000891164 \\
 2.5 & 0.0432581 & 0.0542808 & 0.00212105 & 0.00164952 \\
 3. & 0.0174272 & 0.0521542 & 0.00132515 & 0.00424903 \\
 3.5 & 0.0138333 & 0.045568 & 0.000278168 & 0.00647046 \\
 4. & 0.0494556 & 0.0358424 & 0.000794076 & 0.00804433 \\
 4.5 & 0.0887094 & 0.023884 & 0.00174297 & 0.00883326 \\
 5. & 0.1311 & 0.0103171 & 0.0024838 & 0.00879521 \\
\hline
\end{tabular}
\caption{Relative error of approximate ruin probabilities with Gamma claims
  ($\alpha=2.5,\beta=1$) \label{tab:gamma2error}}
\end{center}
\end{table}

\subsection{With perturbation}

In this section we illustrate the application of the two approximations
given in Theorem~\ref{t:per} and Remark~\ref{r:per}, respectively.

\subsubsection{Mixed exponential claim distribution}

As without perturbation, in this case the exact ruin probabilities can be
calculated by symbolic inversion of the Laplace transform.  We applied the
approximations given in Theorem~\ref{t:per} and Remark~\ref{r:per}  for three values of
$\sigma$, namely, 0.1, 0.5 and 2.  Figures~\ref{fig:rp0.1}-\ref{fig:rp2}
show the exact and approximate ruin probabilities and we depicted the two
components of the ruin probability (by diffusion and by jump) as well.  The two
approximations result in distinguishable curves only in the case
$\sigma=0.1$.  The relative errors for the three values of $\sigma$ are
provided in Figure~\ref{fig:rperr0.1}-\ref{fig:rperr2}.  The error is
smaller for larger values of $\sigma$ and for most values of $x$ the better
approximation of $\Psi} \def\sRui{\overline{\Rui}} \def\ruik{\psi_q(x)$ is by the approach of Theorem~\ref{t:per}.


\begin{figure}
\begin{center}
\begin{minipage}{0.495\textwidth}
\includegraphics[width=\textwidth]
{pme_01_1.pdf}
\end{minipage}
\begin{minipage}{0.495\textwidth}
\includegraphics[width=\textwidth]
{pme_01_2.pdf}
\end{minipage}
\end{center}
\caption{Ruin probabilities with mixed exponential claim distribution with
  $\sigma=0.1$, with the approximation given in Theorem~\ref{t:per} (left)
  and the one given in Remark~\ref{r:per} (right)
\label{fig:rp0.1}}
\end{figure}

\begin{figure}
\begin{center}
\begin{minipage}{0.495\textwidth}
\includegraphics[width=\textwidth]{pme_05.pdf}
\end{minipage}
\end{center}
\caption{Ruin probabilities with mixed exponential claim distribution with
  $\sigma=0.5$, with the  approximations
  of Theorem~\ref{t:per} and Remark~\ref{r:per}  (the two approximations look the same)
\label{fig:rp0.5}}
\end{figure}

\begin{figure}
\begin{center}
\begin{minipage}{0.495\textwidth}
\includegraphics[width=\textwidth]{pme_2.pdf}
\end{minipage}
\end{center}
\caption{Ruin probabilities with mixed exponential claim distribution with
  $\sigma=2$, with the  approximations
  of Theorem~\ref{t:per} and Remark~\ref{r:per}  (the two approximations look the same)
\label{fig:rp2}}
\end{figure}

\begin{figure}
\begin{center}
\includegraphics[width=0.495\textwidth]
{pme_01_error.pdf}
\end{center}
\caption{Absolute relative error of the approximations given in
  Theorem~\ref{t:per} and Remark~\ref{r:per}  with mixed exponential claims and $\sigma=0.1$.
\label{fig:rperr0.1}}
\end{figure}

\begin{figure}
\begin{center}
\includegraphics[width=0.495\textwidth]
{pme_05_error.pdf}
\end{center}
\caption{Absolute relative error of the approximations given in
  Theorem~\ref{t:per} and Remark~\ref{r:per}  with mixed exponential claims and $\sigma=0.5$.
\label{fig:rperr0.5}}
\end{figure}

\begin{figure}
\begin{center}
\includegraphics[width=0.495\textwidth]
{pme_2_error.pdf}
\end{center}
\caption{Absolute relative error of the approximations given in
  Theorem~\ref{t:per} and Remark~\ref{r:per}  with mixed exponential claims and $\sigma=2$.
\label{fig:rperr2}}
\end{figure}




{\bf Acknowledgement:}  We thank Jiandong Ren for useful comments

\bibliographystyle{alpha}
