Ruin probabilities as recurrence sequences in a discrete-time risk process

The Gerber-Dickson risk process ([2017dickson]) is given by

$$U(t)=u+t-\sum_{j=1}^{t}Y_{j},\quad\mbox{for }\ t=0,1,\ldots$$

(1)

This is a simple stochastic process that represents the evolution in time of the capital of an insurance company. Here $U(0)=u\geq 0$ is an integer that stands for the initial capital, the insurance company receives one unit of currency as premium in each unit time period and $Y_{1},Y_{2},\ldots$ are independent, identically distributed discrete random variables taking values on $\{0,1,\ldots\}$. These are the random claims amounts payable at the end of each period and we denote any of them by $Y$. The probability function of the claims is $f(y)={\mathds{P}}(Y=y)$ and the distribution function is $F(y)=P(Y\leq y)$, for $y=0,1,\ldots$ We will assume that the mean value of the claims, denoted by $\mu_{Y}$, satisfies the net profit condition, that is, $\mu_{Y}<1$. When this restriction is not set, in the long run the capital reaches values below zero with probability $1$, which is not desirable. Observe that the net profit condition implies that $f(0)>0$.

The time of ruin is defined as the stopping time $\tau=\min\,\{t\geq 1:U(t)\leq 0\}$, where the minimum of the empty set is defined as infinity. Given an arbitrary distribution for the claims, a central problem in the theory of ruin is to find the probability of ultimate ruin (infinite time horizon), which is defined by

$$\psi(u)={\mathds{P}}(\tau<\infty\mid U(0)=u).\\$$

The risk process (1) is rather elementary and several generalizations have been proposed. For example, let $\{N(t):t=0,1,\ldots\}$ be a counting process such that $N(t)\sim\mbox{Binomial}(t,p)$. The compound binomial risk process is defined as

$$U(t)=u+t-\sum_{j=1}^{N(t)}X_{j},\quad\mbox{for }\ t=0,1,\ldots$$

(2)

where now claims $X_{j}$ take only positive values $1,2,\ldots$ and they are independent and identically distributed as before. The risk process (2) was introduced and studied by H. U. Gerber in ([1988gerber]) and is not really a generalization of (1) since it can be shown (see [Santana-Rincon-2020]) that models (1) and (2) are equivalent in the sense that $U(t)$ has the same distribution in both models when the claims are related by $Y_{j}=R_{j}\,X_{j}$, where $R_{1},R_{2},\ldots$ is a sequence of i.i.d. $\mbox{Bernoulli}(p)$ random variables. As an example of generalization of (1), see [2009Trufin-Loisel] where authors consider a more elaborate discrete-time risk process where experience rating is taken into account. In particular, they derive the logarithmic asymptotic behavior of the ultimate ruin probability for large values of $u$. The probability of ruin in finite horizon where the reserve is invested in a risky asset is studied in [2003Tang-Tsitsiashvili]. For some other interesting extensions of the process (1) see also [2018Liu-Wang-Guo], [2011Peng-Huang-Wang], [cai_2002], [cossette_marceau_maume-deschamps_2010], [diasparra_romera_2009], [Jasiulewicz-Kordecki-2015], [sun_yang_2003], [wu_chen_guo_jin_2015], [YANG-etal200921]. A comprehensive survey of results on several discrete-time risk processes can be found in ([2009Li-Lu-Garrido]).

The basic question of finding the ruin probability $\psi(u)$ for our model (1) and its generalizations remains a hard problem and in this work we give a hint at why this is so. The technique that we here present to solve the problem makes use of the theory of linear recurrence sequences. We have applied this technique to both discrete and continuous time models in ([Rincon-Santana-2021, Santana-Rincon-4, Rincon-Santana-2022]). In those works, however, the recurrence sequences arise from a specific distribution assumed for the claims. The formulas found are thus specialized to that claims distribution. On the contrary, in the present work we will see that the structure of the discrete-time risk process allows to directly identify a recurrence sequence for the ruin probabilities themselves. The claims distribution is only asked to have finite support but else arbitrary. Thus, we will see that solving the recurrence sequence will yield a new formula for the probability of ruin.

Before presenting our results, the reader should be aware that there exist in the actuarial literature several formulas for the ruin probability for discrete-time risk processes. Particularly, for the compound binomial risk process (2), we have the H. U. Gerber’s formula ([1988gerber]). When the definition of time of ruin is modified so that ruin occurs only when the risk process is strictly negative, we have the E. Shiu’s formula ([1989shiu]). We also have a formula by S. Li and J. Garrido ([2002Li]). All those formulas were obtained by different methods and are rather elaborate as they are expressed as infinite series involving the partial sums $S_{k}=X_{1}+\cdots+X_{k}$. These expressions can also be consulted in the excellent survey ([2009Li-Lu-Garrido]).

We here present the procedure that leads to a new expression for the ultimate ruin probability for the risk process (1). Conditioning on the outcome of the first claim (method also known as first step analysis), it is easy to show (see [2017dickson]) that the ruin probability $\psi(u)$ for model (1) satisfies the recursive relation

$$\psi(u)=\sum\limits_{k=0}^{u}\psi(u+1-k)f(k)+\overline{F}(u),\quad\mbox{for }u\geq 0,$$

(3)

where $\overline{F}(u)=1-F(u)$. There are several equivalent forms the relation (3) can be written, another version is the equation

$$\psi(u)=\sum\limits_{k=0}^{u-1}\psi(u-k)\overline{F}(k)+\sum\limits_{k=u}^{\infty}\overline{F}(k),\quad\mbox{for }u\geq 0.$$

(4)

Evaluating (4) at $u=0$ (the empty sum is defined as null) yields the well known result $\psi(0)=\mu_{Y}=\sum_{k=0}^{\infty}\overline{F}(k)$. Thus, in the ensuing calculations we will assume that the value of $\psi(0)$ is known to be the average value of the claims.

Let $m\geq 2$ be an integer and suppose that the claims take only the first $m+1$ values $0,1,\ldots,m$ with $f(m)$ strictly positive. Then we have $\overline{F}(u)=0$ for $u\geq m$ and the relation (3) takes the simpler form

$$\psi(u)=\sum\limits_{k=0}^{m}\psi(u+1-k)f(k),\quad\mbox{for }u\geq m.$$

(5)

Observe that the upper limit of the sum is now the parameter $m$ not $u-1$. Solving for $\psi(u+1)$ gives

$$\psi(u+1)=\frac{1}{f(0)}\,[\ \psi(u)[1-f(1)]-\psi(u-1)f(2)-\cdots-\psi(u-(m-1))f(m)\ ],$$

(6)

for $u\geq m$. This is a linear recurrence sequence of order $m$ for $\psi(u)$ with initial data $\psi(1),\ldots,\psi(m)$. These initial values of the sequence are given by the first terms of (3), namely,

	$\displaystyle\psi(0)$	$\displaystyle=$	$\displaystyle\psi(1)f(0)+\overline{F}(0)$		(7)
	$\displaystyle\psi(1)$	$\displaystyle=$	$\displaystyle\psi(2)f(0)+\psi(1)f(1)+\overline{F}(1)$
	$\displaystyle\psi(2)$	$\displaystyle=$	$\displaystyle\psi(3)f(0)+\psi(2)f(1)+\psi(1)f(2)+\overline{F}(2)$
		$\displaystyle\vdots$
	$\displaystyle\psi(m-1)$	$\displaystyle=$	$\displaystyle\psi(m)f(0)+\psi(m-1)f(1)+\cdots+\psi(1)f(m-1)+\overline{F}(m-1).$

This is a set of $m$ linear equations for the $m$ unknown $\psi(1),\ldots,\psi(m)$, where, as said earlier, $\psi(0)$ is taken as known and equal to $\mu_{Y}$. In matrix form, the system (7) can be written as follows

$$\left(\begin{array}[]{ccccc}f(0)&0&0&\cdots&0\\ f(1)-1&f(0)&0&\cdots&0\\ f(2)&f(1)-1&f(0)&\cdots&0\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ f(m-1)&f(m-2)&f(m-3)&\cdots&f(0)\end{array}\right)\left(\begin{array}[]{c}\psi(1)\\ \psi(2)\\ \psi(3)\\ \vdots\\ \psi(m)\end{array}\right)=\left(\begin{array}[]{c}\psi(0)-\overline{F}(0)\\ -\overline{F}(1)\\ -\overline{F}(2)\\ \vdots\\ -\overline{F}(m-1)\end{array}\right).$$

(8)

The determinant of the above matrix is $(f(0))^{m}$, which is strictly positive by the net profit condition. Hence, a unique collection of values for $\psi(1),\ldots,\psi(m)$ can be determined from (8). Given these initial values, the recursive relation (6) can now be used to find $\psi(u)$ for all $u\geq 1$. To this end, define the $m$ coefficients

$$\alpha_{0}=\frac{1-f(1)}{f(0)},\ \alpha_{1}=-\frac{f(2)}{f(0)},\ \alpha_{2}=-\frac{f(3)}{f(0)},\ \ldots,\ \alpha_{m-1}=-\frac{f(m)}{f(0)}.$$

(9)

Observe that $\alpha_{0}>1$ and $\alpha_{1}\leq 0,\ldots,\alpha_{m-2}\leq 0$ but $\alpha_{m-1}<0$ since $f(m)>0$. Hence, there is exactly one change of sign in the sequence $\alpha_{0},\alpha_{1},\ldots,\alpha_{m-1}$. Observe also that $\alpha_{0}+\alpha_{1}+\cdots+\alpha_{m-1}=1$. Then, the recursive relation (6) can be written as

$$\psi(u+1)=\sum_{k=0}^{m-1}\alpha_{k}\,\psi(u-k),\quad\mbox{for }\ u\geq m,$$

(10)

and writing $u+m$ instead of $u$ yields

$$\psi(u+m+1)=\sum_{k=0}^{m-1}\alpha_{k}\,\psi(u+m-k),\quad\mbox{for }\ u\geq 0.$$

(11)

The next step is to solve the recurrence equation (11). We will do that using the method of characteristic polynomials (see [1971Brousseau, 2013Sedgewick-Flajolet]). The characteristic polynomial associated with (11) is

$$p(y)=y^{m}-\sum_{k=0}^{m-1}\alpha_{k}\,y^{m-1-k}.$$

(12)

Being $p(y)$ an $m$ degree polynomial, it has at most $m$ zeroes, which can be real or complex and can have any multiplicity. Let us suppose that $z_{1},\ldots,z_{\ell}$ are the roots of $p(y)=0$, where $1\leq\ell\leq m$, with multiplicities $n_{1},\ldots,n_{\ell}$, respectively. Observe that $n_{1}+\cdots+n_{\ell}=m$. It is known ([2013Sedgewick-Flajolet]) and this is the important fact, that the general solution to the recurrence sequence (11) is given by

$$\psi(u)=\sum_{k=1}^{\ell}\sum_{j=1}^{n_{k}}b_{k,j}\,u^{j-1}\,z_{k}^{u},\quad\mbox{for }\ u\geq 1,$$

(13)

where $b_{k,j}$ are $m$ constants which can be chosen so that the expression (13) matches the initial data $\psi(1),\ldots,\psi(m)$. From (13), observe that each root $z_{k}$ with multiplicity $n_{k}$ will have $n_{k}$ associated constants $b_{k,1},\ldots,b_{k,n_{k}}$. When the multiplicity of $z_{k}$ is $1$ there is only one constant and we write $b_{k}$ instead of $b_{k,1}$. For example, suppose that $m=6$ and that we have $2$ simple positive roots $z_{1}$ and $z_{2}$, and $2$ roots $z_{3}$ and $z_{4}$ each with multiplicity $2$. The $6$ constants $b_{k,j}$ of (13) that satisfy the initial conditions are given by the solution to the system

$$\left(\begin{array}[]{c|c|cc|cc}z_{1}&z_{2}&z_{3}&z_{3}&z_{4}&z_{4}\\ z_{1}^{2}&z_{2}^{2}&z_{3}^{2}&2z_{3}^{2}&z_{4}^{2}&2z_{4}^{2}\\ z_{1}^{3}&z_{2}^{3}&z_{3}^{3}&3z_{3}^{3}&z_{4}^{3}&3z_{4}^{3}\\ z_{1}^{4}&z_{2}^{4}&z_{3}^{4}&4z_{3}^{4}&z_{4}^{4}&4z_{4}^{4}\\ z_{1}^{5}&z_{2}^{5}&z_{3}^{5}&5z_{3}^{5}&z_{4}^{5}&5z_{4}^{5}\\ z_{1}^{6}&z_{2}^{6}&z_{3}^{6}&6z_{3}^{6}&z_{4}^{6}&6z_{4}^{6}\\ \end{array}\right)\left(\begin{array}[]{c}b_{1}\\ \mbox{-----}\\ b_{2}\\ \mbox{-----}\\ b_{3,1}\\ b_{3,2}\\ \mbox{-----}\\ b_{4,1}\\ b_{4,2}\end{array}\right)=\left(\begin{array}[]{c}\psi(1)\\ \psi(2)\\ \psi(3)\\ \psi(4)\\ \psi(5)\\ \psi(6)\\ \end{array}\right).$$

(14)

The vertical and horizontal lines help visually separate the terms associated with the different roots. In the general case, the $m$ constants $b_{k,j}$ are the solution to the system of linear equations

$$\mathbf{Z}\left(\begin{array}[]{l}b_{1}\\ b_{2}\\ b_{3,-}\\ \vdots\\ b_{\ell,-}\end{array}\right)=\left(\begin{array}[]{c}\psi(1)\\ \psi(2)\\ \psi(3)\\ \vdots\\ \psi(m)\end{array}\right),$$

(15)

with $b_{k,-}$ being the column vector $(b_{k,1},\ldots,b_{k,n_{k}})^{t}$ for $k=3,\ldots,\ell$, and $\mathbf{Z}$ is the $m\times m$ matrix of the form $(\mathbf{\underline{Z}_{1}},\mathbf{\underline{Z}_{2}},\ldots,\mathbf{\underline{Z}_{\ell}})$, where the $k$-entry is given by the following $m\times n_{k}$ submatrix

$$\quad\mathbf{\underline{Z}}_{k}=\left(\begin{array}[]{cccc}z_{k}&z_{k}&\cdots&z_{k}\\ z_{k}^{2}&2z_{k}^{2}&\cdots&2^{n_{k}-1}z_{k}^{2}\\ z_{k}^{3}&3z_{k}^{3}&\cdots&3^{n_{k}-1}z_{k}^{3}\\ \vdots&\vdots&\vdots&\vdots\\ z_{k}^{m}&mz_{k}^{m}&\cdots&m^{n_{k}-1}z_{k}^{m}\\ \end{array}\right)_{m\times n_{k}}.$$

Observe that each root $z_{k}$ has the associated submatrix $\mathbf{\underline{Z}}_{k}$. It is clear that when a root $z_{k}$ is simple ($n_{k}=1$), the submatrix $\mathbf{\underline{Z}}_{k}$ is a column vector. As a summary, we now write the statement of our main result.

In a nutshell, the general procedure to obtain the ruin probability values in (13) is as follows: [a] Given a probability function $f(x)$ for the claims, calculate the coefficients $\alpha_{0},\alpha_{1},\ldots,\alpha_{m-1}$ using (9) and construct the characteristic polynomial (12). [b] Find the roots $z_{1},\ldots,z_{m}$ of (12). [c] Compute the initial values $\psi(1),\ldots,\psi(m)$ using (8). [d] Finally, find the constant $b_{k,j}$ using (15). These steps yield all the terms appearing in our general formula (13).

As a special case, suppose that all the roots of the polynomial (12) are simple, then there are $m$ different roots $z_{1},\ldots,z_{m}$ each with multiplicity $1$ and the ruin probability formula (13) reduces to

$$\psi(u)=\sum_{k=1}^{m}b_{k}\,z_{k}^{u},\quad\mbox{for }\ u\geq 1,$$

(16)

where $b_{1},\ldots,b_{m}$ are the solution to the linear system

$$\left(\begin{array}[]{cccc}z_{1}&z_{2}&\cdots&z_{m}\\ z_{1}^{2}&z_{2}^{2}&\cdots&z_{m}^{2}\\ \vdots&\vdots&\ddots&\vdots\\ z_{1}^{m}&z_{2}^{m}&\cdots&z_{m}^{m}\end{array}\right)\left(\begin{array}[]{c}b_{1}\\ b_{2}\\ \vdots\\ b_{m}\end{array}\right)=\left(\begin{array}[]{c}\psi(1)\\ \psi(2)\\ \vdots\\ \psi(m)\end{array}\right).$$

(17)

In this way we have translated a general problem from the theory of ruin into two classical mathematical problems: finding the roots of a polynomial and solving a system of linear equations. Of course, the main difficulty in applying this method lies in finding the roots of the polynomial.

In this section we make some comments on the main formula (13) and its components. Despite the fact that some of the roots $z_{k}$ and the coefficients $b_{k,j}$ in (13) can be complex, the whole expression for $\psi(u)$ always gives a real number in $(0,1)$, as expected. We will explain why this is so below.

The numerical examples presented in the next section show that the second positive root $z_{2}$ and its associated coefficient $b_{2}$ play a leading role in the ruin probability formula (13). This justifies looking for an approximate solution to the linear system (15) of the form $(0,b_{2},0,\ldots,0)$, which gives the following approximation for the ruin probability.

In the following section we will numerically show how good the approximations $\hat{\psi}_{1}(u)$ and $\hat{\psi}_{2}(u)$ are in some particular cases.

We now present some numerical examples of particular distributions for the claims which yield exact ruin probabilities $\psi(u)$ using formula (13). The roots of the characteristic polynomial (12) and the solution to the linear system (15) are found using the R pracma package ([R-pracma]). In each case observe that the approximation $\hat{\psi}_{1}(u)=b_{2}\,z_{2}^{u}$, for $u\geq 1$, seems to be accurate. This means that the term $b_{2}z_{2}^{u}$ is the one making the major contribution to the exact ruin probability formula (13) for $\psi(u)$.

It is apparent that a correct calculation of the roots of (12) is absolutely essential for providing exact ruin probabilities trough formula (13). When using a computer, numerical errors can lead to inaccurate ruin probabilities. For example, in our numerical experimentation we found that for claims with small mean $\mu_{Y}$, the function polyroots of the R pracma package, sometimes produced some lack of precision in the calculations of the roots of the characteristic polynomial. Hence, for any numerical application of the ruin probability formula (13), it is advisable to make sure the roots are properly calculated. In the examples shown below and to save some space, most of the numbers are rounded up to three or four decimal digits.