2. A dual risk model with causal dependence structure and proportional gains

In this section, we generalize the model in [Reference Boxma, Frostig and Palmowski18] by incorporating a causal dependence structure, namely that the distribution of the gain size depends on the gain interarrival time. This type of dependence it also affects the proportional parameter, which in turn affects the size of the capital jump. In [Reference Albrecher and Boxma6], the authors considered the classical ruin model with causal dependencies (motivated by the work in [Reference Boxma and Perry23]), where the distribution of the interclaim time depends on the actual size of the previous claim based on a (random) threshold type policy. In the following, we consider the dual model with the additional feature of the proportional gains.

Assume that gains arrive according to a renewal process $\{N(t);t\geq 0\}$ with independent and identically distributed (i.i.d.) interarrival times $B_{i}:=S_{i}-S_{i-1}$ ( $S_{0}=0$), that is, the interarrival time between jumps i − 1 and i ( $\geq 1$). We assume that B_i are i.i.d. random variables having cumulative distribution function (c. distribution function) $B(.)$, density $b(.)$, and Laplace–Stieltjes transform $\phi(.)$. If B_i is larger than a threshold T_i, then, the capital jumps up to the level $(1+a_{0})u+C_{i}^{(0)}$, $a_{0} \gt 0$, where $C_{i}^{(0)}$ are i.i.d. random variables that follow a hyperexponential distribution with c. distribution function $C_{0}(x)=\sum_{k=1}^{K}q_{k}(1-e^{-\mu_{j}x})$. Otherwise, the capital jumps up to the level $(1+a_{1})u+C_{i}^{(1)}$, $a_{1} \gt 0$, where $C_{i}^{(1)}$ are i.i.d. random variables that follow a hyperexponential distribution with c. distribution function $C_{1}(x)=\sum_{l=1}^{L}h_{l}(1-e^{-\nu_{l}x})$. The thresholds T_i are assumed to be i.i.d. random variables with c. distribution function $T(.)$. Thus, for $U(0)=x \gt 0$, the surplus process U(t), $t\geq 0$, satisfies:

(2.1)\begin{equation} U(t)=x-t+\sum_{i=1}^{N(t)}[(C_{i}^{(0)}+a_{0}U(S_{i}^{-}))1_{\{B_{i} \gt T_{i}\}}+(C_{i}^{(1)}+a_{1}U(S_{i}^{-}))1_{\{B_{i}\leq T_{i}\}}]. \end{equation}

Contrary to the case in [Reference Boxma, Frostig and Palmowski18] where the authors assumed exponentially distributed gain sizes (they just mentioned in [Reference Boxma, Frostig and Palmowski18, Rem. 2.4] that the model can be generalized to the case of hyperexponentially distributed gain sizes), we consider in Remark 2.1 the case where the gain sizes follow a mixed Erlang distribution, a class of the phase-type distributions that can be used to approximate any given continuous distribution in $[0,\infty)$.

2.1. The ruin probability

We focus on deriving the Laplace transform of the ruin probability R(x) when starting in x, by distinguishing the two cases in which no jump up occurs before x (hence ruin occurs at time x) and in which a jump up occurs at some time $t\in(0,x)$. Then,

(2.2)\begin{equation} \begin{array}{rl} R(x)=&1-B(x)+\int\limits_{t=0}^{x}\int\limits_{y=0}^{\infty}\left(P(T_{i} \lt t)R((1+a_{1})(x-t)+y)\sum_{l=1}^{L}h_{l}\nu_{l}e^{-\nu_{l}y}\right. \\ &\left.+ P(T_{i}\geq t)R((1+a_{0})(x-t)+y)\sum_{k=1}^{K}q_{k}\mu_{k}e^{-\mu_{k}y}\right)dydB(t). \end{array} \end{equation}

Denote the Laplace transform,

\begin{equation*} \rho(s):=\int\limits_{x=0}^{\infty}e^{-sx}R(x)dx. \end{equation*}

Then, (2.2) becomes,

(2.3)\begin{equation} \rho(s)=\frac{1-\phi(s)}{s}+I_{0}(s)+I_{1}(s), \end{equation}

where

\begin{equation*} \begin{array}{rl} I_{0}(s):=&\int\limits_{x=0}^{\infty}e^{-sx}\int\limits_{t=0}^{x}\int\limits_{z=(1+a_{0})(x-t)}^{\infty}P(T_{i}\geq t) R(z)\sum_{k=1}^{K}q_{k}\mu_{k}e^{-\mu_{k}(z-(1+a_{0})(x-t))}dzdB(t)dx \\ =& \int\limits_{t=0}^{\infty}e^{-st}P(T_{i}\geq t) \sum_{k=1}^{K}q_{k}\int\limits_{x=t}^{\infty}e^{-(x-t)[s-\mu_{k}(1+a_{0})]}\int\limits_{z=(1+a_{0})(x-t)}^{\infty}R(z)\mu_{k}e^{-\mu_{k}z}dzdB(t)dx \\ =&\chi_{0}(s)\sum_{k=1}^{K}q_{k}\mu_{k}\int\limits_{w=0}^{\infty}e^{-w[s-\mu_{k}(1+a_{0})]}\int\limits_{z=(1+a_{0})w}^{\infty}R(z)e^{-\mu_{k}z}dzdw \\ =&\chi_{0}(s)\sum_{k=1}^{K}q_{k}\mu_{k}\int\limits_{z=0}^{\infty}R(z)e^{-\mu_{k}z}\left(\frac{e^{-\frac{z}{1+a_{0}}(s-\mu_{k}(1+a_{0}))}-1}{\mu_{k}(1+a_{0})-s}\right)dz\\ =&\chi_{0}(s)\frac{\sum_{k=1}^{K}q_{k}\mu_{k}}{\mu_{k}(1+a_{0})-s}[\rho(\frac{s}{1+a_{0}})-\rho(\mu_{k})], \end{array} \end{equation*}

where $\chi_{0}(s):=E(e^{-sB}1(T\geq B))=\int\limits_{x=0}^{\infty}e^{-sx}(1-T(x))dB(x)$.

Similarly,

\begin{equation*} I_{1}(s):=\chi_{1}(s)\frac{\sum_{l=1}^{L}h_{l}\nu_{l}}{\nu_{l}(1+a_{1})-s}[\rho(\frac{s}{1+a_{1}})-\rho(\nu_{l})], \end{equation*}

where now $\chi_{1}(s):=E(e^{-sB}1(T \lt B))=\int\limits_{x=0}^{\infty}e^{-sx}T(x)dB(x)$. Note that $\chi_{0}(s)+\chi_{1}(s)=\phi(s)$.

By setting $a_{i}(s):= \frac{s}{1+a_{i}}$, $ \widehat{\chi}_{i}(s):=\frac{\chi_{i}(s)}{1+a_{i}}$. $\psi_{i}(s):=\frac{\mu_{i}}{\mu_{i}-s}$, $\omega_{i}(s):=\frac{\nu_{i}}{\nu_{i}-s}$, $i=0,1,$ Eq. (2.3) becomes

(2.4)\begin{equation} \rho(s)=\rho(a_{0}(s))H_{0}(s)+\rho(a_{1}(s))H_{1}(s)+K(s), \end{equation}

where $H_{0}(s):=\widehat{\chi}_{0}(s)\sum_{k=1}^{K}q_{k}\psi_{k}(a_{0}(s))$, $H_{1}(s):=\widehat{\chi}_{1}(s)\sum_{l=1}^{L}h_{i}\omega_{l}(a_{1}(s))$,

\begin{equation*} K(s):=\frac{1-\phi(s)}{s}-\widehat{\chi}_{0}(s)\sum_{k=1}^{K}q_{k}\psi_{k}(a_{0}(s))\rho(\mu_{k})-\widehat{\chi}_{1}(s)\sum_{l=1}^{L}h_{l}\omega_{l}(a_{1}(s))\rho(\nu_{l}). \end{equation*}

Note that $a_{i}(a_{j}(s))=\frac{s}{(1+a_{0})(1+a_{1})}=a_{j}(a_{i}(s))$, thus the mappings $a_{i}(s)$ commute. Moreover, $|a_{i}(s)-a_{j}(u)|\leq \psi|s-u|$, where $\psi:=max\{1/(1+a_{0}),1/(1+a_{1})\}$, so that the mappings $a_{i}(s)$ are contraction mappings.

After n − 1 iterations of (2.4), we have

(2.5)\begin{equation} \rho(s)=\sum_{m=0}^{n}L_{m,n-m}(s)\rho(a_{m,n-m}(s))+\sum_{j=0}^{n-1}\sum_{m=0}^{j}L_{m,j-m}(s)K(a_{m,j-m}(s)), \end{equation}

where $a_{m,j-m}(s):=a_{0}^{m}(a_{1}^{j-m}(s))$, with $a_{0,0}(s)=s$ and $a_{i}^{m}(s)$ is the mth iterate of $a_{i}(s)$, $i=0,1$, $m=0,1,2,\ldots$. Moreover, the functions $L_{m,j-m}(s)$ are computed recursively, with $L_{0,0}(s)=1$, $L_{1,0}(s):=H_{0}(s)$, $L_{0,1}(s):=H_{1}(s)$, and

(2.6)\begin{eqnarray} \begin{array}{rl} L_{m+1,j-m}(s)= & L_{m,j-m}(s) L_{1,0}(a_{m,j-m}(s))+ L_{m+1,j-m-1}(s) L_{0,1}(a_{m+1,j-m-1}(s)),\,j-m\geq m+1, \\ L_{m,j-m+1}(s)= & L_{m,j-m}(s) L_{0,1}(a_{m,j-m}(s))+ L_{m-1,j-m+1}(s) L_{1,0}(a_{m-1,j-m+1}(s)),\,j-m\leq m-1. \end{array}\nonumber\\ \end{eqnarray}

By observing that $\psi:=max\{1/(1+a_{0}),1/(1+a_{1})\} \lt 1$, it is seen following the lines in [Reference Adan, Boxma and Resing3, Sect. 2] that $\rho(a_{m,j-m}(s))=\rho(\frac{s}{(1+a_{0})^{m}(1+a_{1})^{j-m}})$ converges geometrically fast to $\rho(0)=1$. Moreover,

(2.7)\begin{equation} |L_{m,j-m}(s)|\leq \left(\!{j \atop m,j-m}\!\right), \end{equation}

and hence, the sum in the first term of the right hand side in (2.5) is bounded by one, since $L_{m,j-m}(s)$ contains $(\!{j \atop m}\!)$ products of $H_{0}(.)$, $H_{1}(.)$, while $H_{m}(.)$, $m=0,1$, is related to the Laplace–Stieltjes transform of a nonnegative random variable. Thus, the first term in the right hand side in (2.5) converges as $n\to\infty$. On the other hand, $K(a_{m,j-m}(s))$ for large values of j converges to a constant that is smaller than one. Indeed, for large values of j, $a_{m,j-m}(s)\to 0$, and K(0) converges to

\begin{equation*} \bar{b}-\frac{\chi_{0}(0)}{1+a_{0}}\sum_{k=1}^{K}q_{k}\rho(\mu_{k})-\frac{\chi_{1}(0)}{1+a_{1}}\sum_{l=1}^{L}h_{i}\rho(\nu_{l})= \bar{b}-\frac{P(B\leq T)}{1+a_{0}}\sum_{k=1}^{K}q_{k}\rho(\mu_{k})-\frac{P(B \gt T)}{1+a_{1}}\sum_{l=1}^{L}h_{l}\rho(\nu_{l}). \end{equation*}

Note that setting s = 0 in (2.4) results in

\begin{equation*} \begin{array}{rl} \bar{b}-\frac{P(B\leq T)}{1+a_{0}}\sum_{k=1}^{K}q_{k}\rho(\mu_{k})-\frac{P(B \gt T)}{1+a_{1}}\sum_{l=1}^{L}h_{l}\rho(\nu_{l})=&1-\frac{P(B\leq T)}{1+a_{0}}-\frac{P(B \gt T)}{1+a_{1}} \lt 1. \end{array} \end{equation*}

Thus, $|K(a_{m,j-m}(s)-K(0))|\leq\psi^{j}|s|D$, where D is the maximum value of $|K^{\prime}(s)|$. Thus the jth term in the double sum in (2.5) is bounded (having also in mind (2.7)), and thus, the double sum converges for $n\to\infty$. Therefore,

(2.8)\begin{equation} \rho(s)=\lim_{n\to\infty}\sum_{m=0}^{n}L_{m,n-m}(s)+\sum_{j=0}^{\infty}\sum_{m=0}^{j}L_{m,j-m}(s)K(a_{m,j-m}(s)). \end{equation}

We still need to obtain the values of $\rho(\mu_{k})$, $k=1,\ldots,K$, and $\rho(\nu_{l})$, $l=1,\ldots,L$. Setting $s=\mu_{k}$, $k=1,\ldots,K$, and $s=\nu_{l}$, $l=1,\ldots,L$, in (2.8), we can obtain a set of equations to derive these unknowns.

Remark 2.1. We now consider the case where $C_{i}^{(0)}$, $C_{i}^{(1)}$ are i.i.d. random variables that follow a mixed Erlang distribution. More precisely, we assume that when $B_{i}\geq T_{i}$, then the gain sizes $C_{i}^{(0)}$ have c. distribution function

\begin{equation*} C_{0}(x)=\sum_{n_{0}=1}^{N_{0}}k_{n_{0}}(1-e^{-\mu_{0}x}\sum_{l=0}^{n_{0}-1}\frac{(\mu_{0}x)^{l}}{l!}),\,x\geq 0, \end{equation*}

where $\sum_{n_{0}=1}^{N_{0}}k_{n_{0}}=1$. Similarly, when $B_{i} \lt T_{i}$, then the gain sizes $C_{i}^{(1)}$ have c. distribution function

\begin{equation*} C_{1}(x)=\sum_{n_{1}=1}^{N_{0}}k_{n_{1}}(1-e^{-\mu_{1}x}\sum_{l=0}^{n_{1}-1}\frac{(\mu_{1}x)^{l}}{l!}),\,x\geq 0, \end{equation*}

with $\sum_{n_{1}=1}^{N_{1}}k_{n_{1}}=1$. In other words, we assume that the gain sizes follow with probability $k_{n_{m}}$, $n_{m} = 1,\ldots,N_{m}$, $m=0,1,$ an Erlang distribution of scale parameter µ_m and n_m stages.

The class of the phase-type distributions of the above form is dense in the space of distribution functions defined on $[0,\infty)$, and thus, for any distribution function B, there is a sequence B_n of phase-type distributions of this class that converges weakly to B as n goes to infinity; see [Reference Schassberger47]. This class, since it may be used to approximate any given continuous distribution on $[0,\infty)$ arbitrarily close.

The whole analysis can be repeated, although there will be some difficulties. More precisely, by applying the same steps, (2.3) is still valid, but now,

\begin{equation*} \begin{array}{rl} I_{0}(s)= &\chi_{0}(s)\sum_{n_{0}=1}^{N_{0}}k_{n_{0}}\mu_{n_{0}}^{n_{0}} \int\limits_{w=0}^{\infty}e^{-w[s-\mu_{n_{0}}(1+a_{0})]}\int\limits_{z=(1+a_{0})w}^{\infty}R(z)e^{-\mu_{0}z}(z-(1+a_{0})w)^{n_{0}-1}dzdw \\ = & \chi_{0}(s)\sum_{n_{0}=1}^{N_{0}}k_{n_{0}}\mu_{n_{0}}^{n_{0}}\sum_{l=0}^{n_{0}-1}\frac{(-1)^{n_{0}-1-l}(1+a_{0})^{n_{0}-1-l}}{l!(n_{0}-1-l)!}\int\limits_{z=0}^{\infty}z^{l}e^{-\mu_{n_{0}}z}R(z)\int\limits_{w=0}^{\frac{z}{1+a_{0}}}w^{n_{0}-1-l}e^{-w[s-\mu_{n_{0}}(1+a_{0})]}\\ & \times\ dwdz\\ =& \chi_{0}(s)\sum_{n_{0}=1}^{N_{0}}k_{n_{0}}\mu_{n_{0}}^{n_{0}}\sum_{l=0}^{n_{0}-1}\frac{(-1)^{n_{0}-1-l}(1+a_{0})^{n_{0}-1-l}}{l!(s-\mu_{n_{0}}(1+a_{0}))^{n-l}}\int\limits_{z=0}^{\infty}z^{l}e^{-\mu_{0}z}R(z)\\ &\times\left[1-e^{-\frac{z}{1+a_{0}}(s-\mu_{n_{0}}(1+a_{0}))}\sum_{i=0}^{n_{0}-1-l}\frac{(\frac{z}{1+a_{0}}(s-\mu_{n_{0}}(1+a_{0})))^{i}}{i!}\right]dz\\ =&\chi_{0}(s)\sum_{n_{0}=1}^{N_{0}}k_{n_{0}}(-1)^{n_{0}-1}\sum_{l=0}^{n_{0}-1}\left(\frac{\mu_{n_{0}}}{s-(1+a_{0})\mu_{n_{0}}}\right)^{n_{0}-l}\frac{(1+a_{0})^{n_{0}-1-l}}{l!}\rho^{(l)}(\mu_{n_{0}})\\ &-\rho(\frac{s}{1+a_{0}})\chi_{0}(s)\sum_{n_{0}=1}^{N_{0}}k_{n_{0}}(1+a_{0})^{n_{0}-1}(-1)^{n_{0}}\left(\frac{\mu_{n_{0}}}{s-(1+a_{0})\mu_{n_{0}})}\right)^{n_{0}}, \end{array} \end{equation*}

where $\rho^{(l)}(\mu_{n_{0}})$ denotes the lth derivative of $\rho(s)$ at point $s=\mu_{n_{0}}$, $n_{0}=1,\ldots,N_{0}$. Similarly,

\begin{equation*} \begin{array}{rl} I_{1}(s)= &\chi_{1}(s)\sum_{n_{1}=1}^{N_{1}}k_{n_{1}}(-1)^{n_{1}-1}\sum_{l=0}^{n_{1}-1}\left(\frac{\mu_{n_{1}}}{s-(1+a_{0})\mu_{n_{1}}}\right)^{n_{1}-l}\frac{(1+a_{1})^{n_{1}-1-l}}{l!}\rho^{(l)}(\mu_{n_{1}})\\ &-\rho(\frac{s}{1+a_{1}})\chi_{1}(s)\sum_{n_{1}=1}^{N_{1}}k_{n_{1}}(1+a_{1})^{n_{1}-1}(-1)^{n_{1}}\left(\frac{\mu_{n_{1}}}{s-(1+a_{1})\mu_{n_{1}})}\right)^{n_{1}}, \end{array} \end{equation*}

where $\rho^{(l)}(\mu_{n_{1}})$ denotes the lth derivative of $\rho(s)$ at point $s=\mu_{n_{1}}$, $n_{1}=1,\ldots,N_{1}$. Therefore, we come up with the following functional equation:

(2.9)\begin{equation} \rho(s)=\rho(a_{0}(s))\widehat{\chi}_{0}(s)\sum_{n_{0}=1}^{N_{0}}k_{n_{0}}\left(\frac{\mu_{n_{0}}}{\mu_{n_{0}}+a_{0}(s)}\right)^{n_{0}}+\rho(a_{1}(s))\widehat{\chi}_{1}(s)\sum_{n_{1}=1}^{N_{1}}k_{n_{1}}\left(\frac{\mu_{n_{1}}}{\mu_{n_{1}}+a_{1}(s)}\right)^{n_{1}}+K(s), \end{equation}

where now,

\begin{equation*} \begin{array}{rl} K(s):=&\frac{1-\phi(s)}{s}-\widehat{\chi}_{0}(s)\sum_{n_{0}=1}^{N_{0}}k_{n_{0}}\sum_{l=0}^{n_{0}-1}\frac{(-1)^{l}}{l!}\left(\frac{\mu_{n_{0}}}{\mu_{n_{0}}+a_{0}(s)}\right)^{n_{0}-l}\rho^{(l)}(\mu_{n_{0}})\\ &-\widehat{\chi}_{1}(s)\sum_{n_{1}=1}^{N_{1}}k_{n_{1}}\sum_{l=0}^{n_{1}-1}\frac{(-1)^{l}}{l!}\left(\frac{\mu_{n_{1}}}{\mu_{n_{1}}+a_{1}(s)}\right)^{n_{1}-l}\rho^{(l)}(\mu_{n_{1}}).\end{array} \end{equation*}

The form of Eq. (2.9) is the same as the one in (2.4), so we can apply the same iterative procedure to obtain $\rho(s)$. Note, however, that we need to obtain the terms $\rho^{(d_{k})}(\mu_{n_{k}})$, $k=0,1$, $d_{k}=0,1,\ldots,N_{k}-1$, $n_{k}=1,2,\ldots,N_{k}$, that is, the values of the d_kth derivative of $\rho(s)$ at points $s=\mu_{n_{k}}$, $k=0,1$. This task can be accomplished by differentiating the equivalent of (2.8) $0,1,\ldots,N_{k}-1$, $k=0,1,$ times with respect to s, each time followed by a substitution of $s=\mu_{n_{k}}$, respectively. This procedure results in $N_{0}+N_{1}$ equations for these unknowns.

2.2. The time to ruin

We now turn our attention to $\tau_{x}=\inf\{t\geq 0:U(t)=0|U(0)=x\}$, that is, the time to ruin starting at level x. Let $T(x,\alpha):=E(e^{-\alpha\tau_{x}}1(\tau_{x} \lt \infty)|U(0)=x)$ the Laplace transform of the time to ruin when the initial surplus equals x, where α > 0, and $1(.)$ is the indicator function. Then,

\begin{equation*} \begin{array}{rl} T(x,\alpha)=&e^{-ax}(1-B(x))+\int\limits_{t=0}^{x}e^{-\alpha t}\int\limits_{y=0}^{\infty}\left(P(T \lt t) E(e^{\alpha \tau_{(1+a_{0})(x-t)+y}})\sum_{l=1}^{L}h_{l}\nu_{l}e^{-\nu_{l}y}\right.\\ &\left.+P(T\geq t) E(e^{\alpha \tau_{(1+a_{1})(x-t)+y}})\sum_{k=1}^{K}q_{k}\mu_{k}e^{-\mu_{k}y}\right)dydB(t). \end{array} \end{equation*}

Remark 2.2. Note that by taking α = 0 in the expression above, we obtain the expression for the ruin probability R(x) given in (2.2).

Letting,

(2.10)\begin{equation} \tau(s,\alpha):=\int\limits_{x=0}^{\infty}e^{-sx}T(x,\alpha)dx, \end{equation}

we obtain after some algebra:

(2.11)\begin{equation} \tau(s,\alpha)=\tau(a_{0}(s),\alpha)H_{0}(s,\alpha)+\tau(a_{1}(s),\alpha)H_{1}(s,\alpha)+\tilde{K}(s,\alpha), \end{equation}

where,

\begin{equation*} \begin{array}{rl} H_{0}(s,\alpha)=&\widehat{\chi}_{0}(s+\alpha)\sum_{k=1}^{K}q_{k}\psi_{k}(a_{0}(s)), \\ H_{1}(s,\alpha)=&\widehat{\chi}_{1}(s+\alpha)\sum_{l=1}^{L}h_{l}\omega_{l}(a_{1}(s)),\\ \tilde{K}(s,\alpha)=&\frac{1-\phi(s+\alpha)}{s+\alpha}-\widehat{\chi}_{0}(s+\alpha)\sum_{k=1}^{K}q_{k}\psi_{k}(a_{0}(s))\tau(\mu_{k},\alpha)-\widehat{\chi}_{1}(s+\alpha)\sum_{l=1}^{L}h_{l}\omega_{l}(a_{1}(s))\tau(\nu_{l},\alpha). \end{array} \end{equation*}

After n − 1 iterations, (2.11) yields

(2.12)\begin{equation} \tau(s,\alpha)=\sum_{l=0}^{n}\tilde{L}_{l,n-l}(s,a)\tau(a_{l,n-l}(s),\alpha)+\sum_{j=0}^{n-1}\sum_{l=0}^{j}\tilde{L}_{l,j-l}(s,\alpha)\tilde{K}(a_{l,j-l}(s),\alpha), \end{equation}

where now the functions $\tilde{L}_{l,j-l}(s,\alpha)$ are computed recursively, with $\tilde{L}_{0,0}(s,\alpha)=1$, $\tilde{L}_{1,0}(s,\alpha):= H_{0}(s,\alpha)$, $\tilde{L}_{0,1}(s,\alpha):= H_{1}(s,\alpha)$, and

\begin{equation*} \begin{array}{rl} \tilde{L}_{l+1,j-l}(s,\alpha)= & \tilde{L}_{l,j-l}(s,\alpha) \tilde{L}_{1,0}(a_{l,j-l}(s),\alpha)+ \tilde{L}_{l+1,j-l-1}(s,\alpha) \tilde{L}_{0,1}(a_{l+1,j-l-1}(s),\alpha),\,j-l\geq l+1, \\ \tilde{L}_{l,j-l+1}(s,\alpha)= & \tilde{L}_{l,j-l}(s,\alpha) \tilde{L}_{0,1}(a_{l,j-l}(s),\alpha)+ \tilde{L}_{l-1,j-l+1}(s,\alpha) \tilde{L}_{1,0}(a_{l-1,j-l+1}(s),\alpha),\,j-l\leq l-1. \end{array} \end{equation*}

For large j, $\tilde{K}(a_{l,j-l}(s),\alpha)$ approaches a function of α. Following similar arguments as those in subsection 2.1, we ensure the convergence of the infinite sums of products, as $n\to\infty$, in (2.12), and obtain

(2.13)\begin{equation} \tau(s,\alpha)=\lim_{n\to\infty}\sum_{l=0}^{n}\tilde{L}_{l,n-l}(s,\alpha)+\sum_{j=0}^{\infty}\sum_{l=0}^{j}\tilde{L}_{l,j-l}(s,\alpha)\tilde{K}(a_{l,j-l}(s),\alpha). \end{equation}

The values of $\tau(\mu_{j})$, $j=1,\ldots,K$, $\tau(\nu_{i})$, $i=1,\ldots,L$, can be obtained by solving a system of K + L equations, which is derived by substituting $s=\mu_{j}$, $j=1,\ldots,K$, and $s=\nu_{i}$, $i=1,\ldots,L$, in (2.13).

Remark 2.3. Similar arguments can be used to cope with the case where $C_{i}^{k}$, $k=0,1,$ are i.i.d. random variables that follow a mixed Erlang distribution, although it would be trickier due to the presence of additional unknown terms that corresponds to the derivatives of $\tau(s,a)$ at specific points (see Remark 2.1), and further details are omitted.

3. The dual risk model with additive/proportional gains and dependencies based on FGM copula

We now generalize the model in [Reference Boxma, Frostig and Palmowski18] by assuming that the gain interarrival times B_i, and the gain sizes C_i are dependent based on the FGM copula. As indicated in [Reference Nelsen39], the FGM copula is a perturbation of the product copula and a first order approximation to the Ali Mikhail Haq, Frank, and Placket copulas. A dependent structure based on copulas has been considered so far in the classical risk reserve process (e.g., [Reference Cossette, Marceau and Marri27]), but to our best knowledge there is a lack of results regarding the use of copulas even in the standard (i.e., where a = 0) dual risk model. In this work we fill this gap, by considering copulas to describe the dependence structure in this general dual risk model with proportional gains.

In [Reference Cossette, Marceau and Marri27], the authors considered an extension of the classical compound Poisson risk model where the claim amounts and the claim inter-arrival times are dependent through an FGM copula. In their work, they assumed that the interclaim arrival times are exponentially distributed. Later, in [Reference Chadjiconstantinidis and Vrontos24], the authors generalized the work in [Reference Cossette, Marceau and Marri27], by considering an extension to the renewal or Sparre Andersen risk process, where the claim interarrivals are Erlang distributed. Later, in [Reference Cossette, Marceau and Marri26], the authors considered an extension to the classical compound Poisson risk model in which, the dependence structure between the claim amounts and the interclaim time is embedded via a generalized FGM (GFGM) copula. We also mentioned the recent work in [Reference Boxma and Mandjes21], where some risk models with a semi-linear dependent structure (no copula dependent) were discussed.

As stated in [Reference Nelsen39], “copulas are functions which ‘join’ or ‘couple’ multivariate distribution functions to their one-dimensional marginal distribution functions.” In other words, a copula itself is a multivariate distribution function whose inputs are the respective marginal cumulative probability distribution functions for the random variables of interest. A bivariate copula C is a joint distribution function on $[0,1]\times[0,1]$ with standard uniform marginal distributions. Sklar’s theorem (see, e.g., [Reference Nelsen39]) states that any bivariate distribution function F with marginals F ₁ and F ₂ can be written as $F(x,y)=C(F_{1}(x),F_{2}(y))$, for some copula C. This copula is unique if F is continuous. Otherwise, it is uniquely defined on the range of the marginals. We refer the reader to [Reference Joe35, Reference Nelsen39] for further details on copulas. Modeling the dependence structure between random variables (r.v.) using copulas has become popular in actuarial science and financial risk management. The reader may refer to, for example, [Reference Albrecher and Teugels8, Reference Denuit, Dhaene, Goovaerts and Kaas28, Reference McNeil, Frey and Embrechts38] for applications of copulas in actuarial science and financial risk management.

In this section, we generalize the model in [Reference Boxma, Frostig and Palmowski18], by considering a dependence structure among interarrival gains, say B_i, and gain sizes C_i based on the FGM copula. We assume that B_i are i.i.d. random variables having c. distribution function $B(.)$, density $b(.)$, and Laplace–Stieltjes transform $\phi(.)$, and assuming C_i to follow Erlang distribution. In particular, we assume that $\{(B_{i},C_{i}),i\geq 1\}$ form a sequence of i.i.d. random vectors distributed as the canonical r.v. (B, C), in which the components are dependent. The joint density of (B, C) is denoted by $f_{B,C}(x,t)$, $x,t\geq 0$, the c. distribution function by $F_{B,C}(x,t)$. The FGM copula is used to define the joint distribution of (B, C). The FGM copula is defined by

\begin{equation*} C_{\theta}^{FGM}(u,v):=uv+\theta uv(1-u)(1-v), \end{equation*}

for $(u,v)\in[0,1]^{2}$, and $\theta\in(-1,1)$. The FGM copula allows positive and negative dependence, and it further includes the independence copula for θ = 0. Its tractability, and simplicity, makes FGM copula quite useful in applications to describe dependence structures.

The density associated with the above expression is $c_{\theta}^{FGM}(u,v):=\frac{\partial^{2}}{\partial u\partial v}C_{\theta}^{FGM}(u,v)=1+\theta(1-2u)(1-2v)$, and then, the bivariate density of (B, C) is given by

\begin{equation*} \begin{array}{rl} f_{B,C}(x,t)=&c_{\theta}^{FGM}(F_{B}(x),F_{C}(t))f_{B}(x)f_{C}(t)= f_{B}(x)f_{C}(t)+\theta h(x)(2\bar{F}_{C}(t)-1), \end{array} \end{equation*}

with $h(x):=f_{B}(x)(1-2F_{B}(x))$ with Laplace transform $h^{*}(s)=\int_{0}^{\infty}e^{-sx}h(x)dx$, and $\bar{F}_{C}(t):=1-F_{C}(t)$. In our case,

(3.1)\begin{equation} \begin{array}{rl} f_{B,C}(x,t)=&f_{B}(x)\frac{\mu^{n}}{(n-1)!}t^{n-1}e^{-\mu t}\\ &+\theta h(x)\left[2\frac{\mu^{n}}{(n-1)!}t^{n-1}e^{-2\mu t}(\sum_{i=0}^{n-1}\frac{(\mu t)^{i}}{i!})-\frac{\mu^{n}}{(n-1)!}t^{n-1}e^{-\mu t}\right]. \end{array} \end{equation}

Therefore,

(3.2)\begin{equation} \begin{array}{rl} R(x)=&1-B(x)+\int\limits_{t=0}^{x}\int\limits_{y=0}^{\infty}R((1+a)(x-t)+y)f_{B,C}(t,y)dydt. \end{array} \end{equation}

Taking Laplace transform, we come up with the following expression

\begin{equation*} \rho(s)=\frac{1-\phi(s)}{s}+\phi(s)I_{0}(s)+I_{1}(s)-I_{2}(s), \end{equation*}

where

\begin{equation*} \begin{array}{rl} I_{0}(s)=&\phi(s)\int\limits_{x=t}^{\infty}e^{-s(x-t)}\int\limits_{z=(1+a)(x-t)}^{\infty}R(z)\mu^{n}\frac{(z-(1+a)(x-t))^{n-1}}{(n-1)!}e^{-\mu(z-(1+a)(x-t)}dzdx \\ = &\sum_{l=0}^{n-1}\frac{(-1)^{n-1}(1+a)^{n-1-l}\mu^{n-l}}{l!(s-\mu(1+a))^{n-l}}\rho^{(l)}(\mu)+\frac{1}{1+a}\left(\frac{\mu(1+a)}{\mu(1+a)-s}\right)^{n}\rho(\frac{s}{1+a}),\\ I_{1}(s)=&\theta h^{*}(s)\int\limits_{x=t}^{\infty}e^{-s(x-t)}\int\limits_{z=(1+a)(x-t)}^{\infty}\frac{2R(z)\mu^{n}(z-(1+a)(x-t))^{n-1}}{(n-1)!}e^{-2\mu(z-(1+a)(x-t))}\\ &\times\sum_{i=0}^{n-1}\frac{(\mu(z-(1+a)(x-t)))^{i}}{i!}dzdx\\ =&2\theta h^{*}(s)\sum_{i=0}^{n-1}\left(\!{n-1+i \atop i}\!\right)\left[\sum_{l=0}^{n-1+i}\frac{(-1)^{n+i-1}}{l!2^{l}(1+a)}\left(\frac{\mu(1+a)}{s-2\mu(1+a)}\right)^{n+i}\rho^{(l)}(2\mu)\right.\\ &\left.+\frac{1}{1+a}\left(\frac{\mu(1+a)}{2\mu(1+a)-s}\right)^{n+i}\rho(\frac{s}{1+a})\right],\\ I_{2}(s)=&\theta h^{*}(s)I_{0}(s). \end{array} \end{equation*}

Combining the above, we come up with the following functional equation:

(3.3)\begin{equation} \rho(s)=J(s;\theta)\rho(\frac{s}{1+a})+H(s;\theta), \end{equation}

where

(3.4)\begin{equation} \begin{array}{rl} J(s;\theta)=&\frac{\phi(s)-\theta h^{*}(s)}{1+a}\left(\frac{\mu(1+a)}{\mu(1+a)-s}\right)^{n}+ \frac{2\theta h^{*}(s)}{1+a}\sum_{i=0}^{n-1}\left(\!{n-1+i \atop i}\!\right)\left(\frac{\mu(1+a)}{2\mu(1+a)-s}\right)^{n+i},\\ H(s;\theta)= & \frac{1-\phi(s)}{s}+\widehat{H}(s;\theta)\\=& \frac{1-\phi(s)}{s}+\frac{(-1)^{n-1}(\phi(s)-\theta h^{*}(s))}{1+a}\sum_{l=0}^{n-1}\frac{\rho^{(l)}(\mu)}{l!}\left(\frac{\mu(1+a)}{s-\mu(1+a)}\right)^{n-l}\\ &+2\theta h^{*}(s)\sum_{i=0}^{n-1}\left(\!{n-1+i\atop i}\!\right)\sum_{l=0}^{n-1+i}\frac{(-1)^{n+i-1}}{l!2^{l}(1+a)}\left(\frac{\mu(1+a)}{s-2\mu(1+a)}\right)^{n+i}\rho^{(l)}(2\mu). \end{array} \end{equation}

Note that (3.3) has the same form as in [Reference Boxma, Frostig and Palmowski18, Eq. (2.5)], so we can apply similar arguments to solve it. Iterating (3.3) N − 1 times results in

\begin{equation*} \rho(s)=\sum_{k=0}^{N-1}\prod_{j=0}^{k-1}J(\frac{s}{(1+a)^{j}} ;\theta)H(\frac{s}{(1+a)^{k}};\theta)+\rho(\frac{s}{(1+a)^{N}})\prod_{j=0}^{N-1}J(\frac{s}{(1+a)^{j}};\theta), \end{equation*}

with an empty product being equal to 1. Observe that for large k, $J(\frac{s}{(1+a)^{k}} ;\theta)$ approaches

\begin{equation*} \begin{array}{l} \frac{\phi(0)-\theta h^{*}(0)}{1+a}+\frac{2\theta h^{*}(0)}{1+a}\sum_{i=0}^{n-1}\left(\!{n-1+i\atop i}\!\right)(\frac{1}{2})^{n+i}=\frac{\phi(0)-\theta h^{*}(0)}{1+a}+\frac{\theta h^{*}(0)}{1+a}(\frac{1}{2}+\frac{1}{2})^{n-1}=\frac{\phi(0)}{1+a} \lt 1, \end{array} \end{equation*}

while $H(\frac{s}{(1+a)^{k}})$ approaches a constant (as $k\to\infty$). Thus,

(3.5)\begin{equation} \rho(s)=\sum_{k=0}^{\infty}\prod_{j=0}^{k-1}J(\frac{s}{(1+a)^{j}} ;\theta)H(\frac{s}{(1+a)^{k}};\theta). \end{equation}

The terms $\rho^{(l)}(\mu)$, $l=0,1,\ldots,n-1$, $\rho^{(l)}(2\mu)$, $l=0,1,\ldots,2n-2$, are obtained as follows: First, we simply differentiating l times (3.5) with respect to s, and setting $s=\mu$, and $s=2\mu$, respectively. Then, a system of $2n+1$ equations is derived with unknown the terms $\rho^{(l)}(\mu)$, $l=0,1,\ldots,n-1$, $\rho^{(l)}(2\mu)$, $l=0,1,\ldots,2n-2$.

Then, we summarize the above in the following main result:

3.1. The time to ruin

Working on the same spirit as in subsection 2.2, we focus on $\tau(s,\alpha)$ defined in (2.10), and we obtain

\begin{equation*} \tau(s,\alpha)=\frac{1-\phi(s+\alpha)}{s+\alpha}+\phi(s+\alpha)P_{0}(s,\alpha,\theta)+P_{1}(s,\alpha;\theta)-P_{2}(s,\alpha;\theta), \end{equation*}

where now,

\begin{equation*} \begin{array}{rl} P_{0}(s,\alpha,\theta)=&\frac{1}{1+a}\left\{\sum_{i=0}^{n-1}\frac{(-1)^{n-1}}{i!}(\frac{\mu(1+a)}{s-\mu(1+a)})^{n-i}\tau^{(i)}(\mu,\alpha)+(\frac{\mu(1+a)}{s-\mu(1+a)})^{n}\tau(\frac{s}{1+a},\alpha)\right\}, \\ P_{1}(s,\alpha,\theta)=& 2\theta h^{*}(s+\alpha)\sum_{i=0}^{n-1}\left(\!{n-1+i\atop i}\!\right)\left[\sum_{l=0}^{n-1+i}\frac{(-1)^{n+i-1}}{l!2^{l}(1+a)}\left(\frac{\mu(1+a)}{s-2\mu(1+a)}\right)^{n+i}\tau^{(l)}(2\mu,\alpha)\right.\\ &\left.+\frac{1}{1+a}\left(\frac{\mu(1+a)}{2\mu(1+a)-s}\right)^{n+i}\tau(\frac{s}{1+a},\alpha)\right],\\ P_{2}(s,\alpha;\theta)=&\theta h^{*}(s+\alpha)P_{0}(s,\alpha;\theta). \end{array} \end{equation*}

Substituting back, we finally obtain

(3.6)\begin{equation} \tau(s,\alpha)=\tau(\frac{s}{1+a},\alpha)\tilde{J}(s,\alpha;\theta)+\tilde{H}(s,\alpha;\theta), \end{equation}

where

\begin{equation*} \begin{array}{rl} \tilde{J}(s,\alpha;\theta)=& \frac{\phi(s+\alpha)-\theta h^{*}(s+\alpha)}{1+a}\left(\frac{\mu(1+a)}{\mu(1+a)-s}\right)^{n}+ \frac{2\theta h^{*}(s+\alpha)}{1+a}\sum_{i=0}^{n-1}\left(\!{n-1+i\atop i}\!\right)\left(\frac{\mu(1+a)}{2\mu(1+a)-s}\right)^{n+i}, \\ \tilde{H}(s,\alpha;\theta)= & \frac{1-\phi(s+\alpha)}{s+\alpha}+\frac{(-1)^{n-1}(\phi(s+\alpha)-\theta h^{*}(s+\alpha))}{1+a}\sum_{l=0}^{n-1}\frac{\tau^{(l)}(\mu,\alpha)}{l!}\left(\frac{\mu(1+a)}{s-\mu(1+a)}\right)^{n-l}\\ &+2\theta h^{*}(s+\alpha)\sum_{i=0}^{n-1}\left(\!{n-1+i \atop i}\!\right)\sum_{l=0}^{n-1+i}\frac{(-1)^{n+i-1}}{l!2^{l}(1+a)}\left(\frac{\mu(1+a)}{s-2\mu(1+a)}\right)^{n+i}\tau^{(l)}(2\mu,\alpha). \end{array} \end{equation*}

Remark 3.4. Consider the case where when the ith jump upward occur (while the surplus just before is u), the jump size is $au+C_{i}$ with probability p, and D_i with probability $q:=1-p$. However, there is a dependence based on FGM copula. In particular, we assume that C_i (resp. D_i) is Erlang distributed with parameters n, µ (resp. m, ν). The joint bivariate pdf of (B, C) is given by (3.1) (where $\theta\equiv\theta_{1}$), while the one of (B, D) by

(3.7)\begin{equation} \begin{array}{rl} f_{B,D}(x,t)=&f_{B}(x)\frac{\nu^{n}}{(m-1)!}t^{m-1}e^{-\nu t}\\ &+\theta_{2} h(x)\left[2\frac{\nu^{n}}{(m-1)!}t^{m-1}e^{-2\nu t}(\sum_{i=0}^{m-1}\frac{(\nu t)^{i}}{i!})-\frac{\nu^{n}}{(m-1)!}t^{m-1}e^{-\nu t}\right], \end{array} \end{equation}

with $\theta_{2}\in(-1,1)$. Then,

(3.8)\begin{equation} \begin{array}{rl} R(x)=&1-B(x)+p\int\limits_{t=0}^{x}\int\limits_{y=0}^{\infty}R((1+a)(x-t)+y)f_{B,C}(t,y)dydt,\\ &+q\int\limits_{t=0}^{x}\int\limits_{y=0}^{\infty}R(x-t+y)f_{B,D}(t,y)dydt. \end{array} \end{equation}

By taking p = 1, we get the previous model, while for a = 0, we get the classical dual risk model with dependence based on the FGM copula. Using similar calculations as those that leads to (3.3), we obtain after heavy but straightforward calculations

(3.9)\begin{equation} \rho(s)[1-q\tilde{J}(s;\theta_{2})]=\frac{1-\phi(s)}{s}+pJ(s;\theta_{1})\rho(\frac{s}{1+a})+p \widehat{H}(s;\theta_{1})+q\tilde{H}(s;\theta_{2}), \end{equation}

where $J(s;\theta_{1})$, $\widehat{H}(s,\theta_{1})$ as given in (3.4) for $\theta=\theta_{1}$, and

(3.10)\begin{equation} \begin{array}{rl} \tilde{J}(s;\theta_{2})= &(\phi(s)-\theta_{2}h^{*}(s))(\frac{\nu}{\nu-s})^{m}+2\theta_{2} h^{*}(s)\sum_{i=0}^{m-1}\left(\!{m-1+i\atop i}\!\right) (\frac{\nu}{2\nu-s})^{m+i}, \\ \tilde{H}(s;\theta_{2})= & (\phi(s)-\theta_{2}h^{*}(s))(-1)^{m-1}\sum_{i=0}^{m-1}(\frac{\nu}{s-\nu})^{m-i}\frac{\rho^{(i)}(\nu)}{i!}\\ &+2\theta_{2}h^{*}(s)\sum_{i=0}^{m-1}\left(\!{m-1+i\atop i}\!\right) \sum_{l=0}^{m-1+i}(\frac{\nu}{s-2\nu})^{m+i-l}\frac{(-1)^{m+i-l}}{l!2^{l}}\rho^{(l)}(2\nu). \end{array} \end{equation}

After simple algebraic arguments, (3.9) is rewritten as

(3.11)\begin{equation} \rho(s)=J_{1}(s;\theta_{1},\theta_{2})\rho(\frac{s}{1+a})+H_{1}(s;\theta_{1},\theta_{2}), \end{equation}

where

\begin{equation*} \begin{array}{rl} J_{1}(s;\theta_{1},\theta_{2})=&\frac{p J(s;\theta_{1})}{1-q\tilde{J}(s;\theta_{2})}, \\ H_{1}(s;\theta_{1},\theta_{2})= &\frac{\frac{1-\phi(s)}{s}+p\widehat{H}(s;\theta_{1})+q\tilde{H}(s;\theta_{2})}{1-q\tilde{J}(s;\theta_{2})}. \end{array} \end{equation*}

To proceed, we need to investigate the number of zeroes of $1-q\tilde{J}(s;\theta_{2})$ in the right-half complex plane.

Proposition 3.5. For $\theta_{2}\neq 0$, the equation $q\tilde{J}(s;\theta_{2})=1$ has exactly $3m-1$ roots, say $s_{1},s_{2},\ldots,s_{3m-1}$ in the right-half complex plane, that is, $Re(s_{j}) \gt 0$, $j=1,2,\ldots,3m-1$.

Proof. The proof is based on Rouché’s theorem [Reference Titchmarsh48]. Since $q\tilde{J}(s;\theta_{2})=1$ can be rewritten as

(3.12)\begin{eqnarray} \begin{array}{r} q\left[\nu^{m}\phi(s)(2\nu-s)^{2m-1}+\theta\nu^{m}h^{*}(s)\left(2\sum_{i=0}^{m-1}\left(\!{m+i-1 \atop i}\!\right)\nu^{i}(\nu-s)^{m}(2\nu-s)^{m-i-1}-(2\nu-s)^{2m-1}\right)\right]\\ =(\nu-s)^{m}(2\nu-s)^{2m-1}, \end{array}\nonumber\\ \end{eqnarray}

it suffices to show that (3.12) has exactly $3m-1$ roots with positive real parts. Let r > 0 be a sufficiently large number, and denote by C_r the contour containing the imaginary axis running from ir to −ir and a semicircle with radius r running clockwise from ir to −ir, that is, $C_{r}=\{s\in\mathbb{C}: |s|=r, Re(s)\geq 0, r \gt 0\,\text{fixed}\}$. Denote by C the limiting contour, by letting $r\to\infty$. We distinguish two cases according to as $Re(s) \gt 0$ or $Re(s)=0$.

For s, such that $Re(s) \gt 0$, we have that $|2\nu-s|\to\infty$, and $|\nu-s|\to\infty$ as $r\to\infty$. Thus,

\begin{equation*} \begin{array}{rl} |q\tilde{J}(s;\theta_{2})|\leq &q\left[|\phi(s)|\frac{\nu^{m}}{|\nu-s|^{m}}+|\theta_{2}||h^{*}(s)|(2\sum_{i=0}^{m-1}\left(\!{m+i-1\atop i}\!\right))\frac{\nu^{m+i}}{|2\nu-s|^{m+i}}+\frac{\nu^{m}}{|\nu-s|^{m}})\right]\to 0, \end{array} \end{equation*}

on C, that is, as $r\to \infty$, and thus, $|q\tilde{J}(s;\theta_{2})| \lt 1$, on C.

Denote by $d(s):=2\sum_{i=0}^{m-1}\left(\!{m+i-1\atop i}\!\right))(\frac{\nu}{2\nu-s})^{m+i}-(\frac{\nu}{\nu-s})^{m}$. Note that,

\begin{equation*} d(0)=2\sum_{i=0}^{m-1}\left(\!{m+i-1\atop i}\!\right))(\frac{1}{2})^{m+i}-1=2(\frac{1}{2})^{m}\sum_{i=0}^{m-1}\left(\!{m+i-1\atop i}\!\right))(\frac{1}{2})^{i}-1=2(\frac{1}{2})^{m}2^{m-1}-1=0. \end{equation*}

For s such that $Re(s)=0$,

\begin{equation*} \begin{array}{rl} |q\tilde{J}(s;\theta_{2})|\leq &q\left[|\phi(s)|\frac{\nu^{m}}{|\nu-s|^{m}}+|\theta_{2}||h^{*}(s)|(2\sum_{i=0}^{m-1}\left(\!{m+i-1\atop i}\!\right))\frac{\nu^{m+i}}{|2\nu-s|^{m+i}}+\frac{\nu^{m}}{|\nu-s|^{m}})\right]\\ \leq& q\left[\frac{\nu^{m}}{|\nu-s|^{m}}+|\theta_{2}||d(s)|\right]\leq q(1+|\theta_{2}||d(0)|)=q \lt 1. \end{array} \end{equation*}

Thus, in each case, we proved that $|q\tilde{J}(s;\theta_{2})| \lt 1$, or equivalently,

\begin{equation*} \begin{array}{r} |q\left[\nu^{m}\phi(s)(2\nu-s)^{2m-1}+\theta\nu^{m}h^{*}(s)\left(2\sum_{i=0}^{m-1}\left(\!{m+i-1\atop i}\!\right)\nu^{i}(\nu-s)^{m}(2\nu-s)^{m-i-1}-(2\nu-s)^{2m-1}\right)\right]|\\ \lt |(\nu-s)^{m}(2\nu-s)^{2m-1}|, \end{array} \end{equation*}

thus by Rouché’s theorem [Reference Titchmarsh48], it follows that (3.12) has the same number of roots with $(\nu-s)^{m}(2\nu-s)^{2m-1}=0$ inside C_r. Since the latter equation has exactly $3m-1$ positive roots inside C_r (with $r\to\infty$), we deduce that (3.12), or equivalently $1-q\tilde{J}(s;\theta_{2})=0$ has exactly $3m-1$ roots, say $s_{1},\ldots,s_{3m-1}$ with positive real parts.

Iterating (3.11) N − 1 times yields

\begin{equation*} \rho(s)=\sum_{k=0}^{N-1}\prod_{j=0}^{k-1}J_{1}(\frac{s}{(1+a)^{j}} ;\theta_{1},\theta_{2})H_{1}(\frac{s}{(1+a)^{k}};\theta_{1},\theta_{2})+\rho(\frac{s}{(1+a)^{N}})\prod_{j=0}^{N-1}J(\frac{s}{(1+a)^{j}};\theta_{1},\theta_{2}). \end{equation*}

Using similar arguments as above, we conclude that

(3.13)\begin{equation} \rho(s)=\sum_{k=0}^{\infty}\prod_{j=0}^{k-1}J_{1}(\frac{s}{(1+a)^{j}} ;\theta_{1},\theta_{2})H_{1}(\frac{s}{(1+a)^{k}};\theta_{1},\theta_{2}). \end{equation}

We still need to obtain $\rho^{(l)}(\mu)$, $l=0,1,\ldots,n-1$, $\rho^{(l)}(2\mu)$ $l=0,1,\ldots,2n-2$, and $\rho^{(d)}(\nu)$, $d=0,1,\ldots,m-1$, $\rho^{(d)}(2\nu)$, $d=0,1,\ldots,2m-2$. Differentiating l times ( $l=0,1,\ldots,n-1$) (3.13) and substituting $s=\mu$, we obtain n equations. Similarly, by differentiating l times ( $l=0,1,\ldots,2n-2$) (3.13) and substituting $s=2\mu$, we have other $2n-1$. We still need other $3m-1$ equations. To this end, we invoke Rouché’s theorem [Reference Titchmarsh48] and proposition 3.5, where we showed that $1-p\tilde{J}(s;\theta_{2})=0$ has exactly $3m-1$ roots with positive real parts. Setting $s=s_{j}$, $j=1,\ldots,3m-1$ in (3.9) yields (since $\rho(s)$ is analytic in the right half s-plane, $\rho(s_{j})$, $j=1,\ldots,3m-1$ is finite)

(3.14)\begin{equation} \frac{1-\phi(s_{j})}{s_{j}}+pJ(s_{j};\theta_{1})\rho(\frac{s_{j}}{1+a})+p\widehat{H}(s_{j};\theta_{1})+q\tilde{H}(s_{j};\theta_{2})=0,\,j=1,\ldots,3m-1. \end{equation}

Definitely (3.14) provide $3m-1$ equations, but they also introduce the additional unknowns $\rho(\frac{s_{j}}{1+a})$, $j=1,\ldots,3m-1$. However, this problem is overcome by substituting $s=\frac{s_{j}}{1+a}$ in (3.13), so that another equation is derived (for each $j=1,\ldots,3m-1$) where $\rho(\frac{s_{j}}{1+a})$ is given in terms of the unknowns $\rho^{(l)}(\mu)$, $\rho^{(l)}(2\mu)$, $\rho^{(d)}(\nu)$, $\rho^{(d)}(2\nu)$.

3.2. An extension to the GFGM copula

We now assume that the gain interarrivals and the gain sizes are dependent based on a GFGM copula, which belongs to a family of copulas introduced by [Reference Rodrıguez-Lallena and Úbeda-Flores46], and defined by

\begin{equation*} C(u, v) = uv + \theta p(u)g(v), \end{equation*}

where $p(.)$ and $g(.)$ are non-zero real functions with support $[0, 1]$. For other extensions of the classical FGM copula see, for example, [Reference Drouet Mari and Kotz30]. Our motivation on the extensions of FGM is to improve the range of dependence association (as measured by either Kendall’s τ or Spearman’s ρ) between the components of (B, C). In this work, we assume the case where $p(u):=u^{k}(1-u)^{b}$, $g(v):=v^{c}(1-v)^{d}$, with $k,b,c,d\geq 1$. The density associated with the GFGM is given by $c(u,v)=1+\theta p^{\prime}(u)g^{\prime}(v)$. Thus, the joint density of the random vectors $(B_{i},C_{i})$ is given by

\begin{equation*} \begin{array}{rl} f_{B,C}(x,t)=&c(F_{B}(x),F_{C}(t))f_{B}(x)f_{C}(t)= f_{B}(x)f_{C}(t)+\theta p^{\prime}(F_{B}(x))g^{\prime}(F_{C}(t))f_{B}(x)f_{C}(t). \end{array} \end{equation*}

Therefore,

(3.15)\begin{equation} \begin{array}{rl} R(x)=&1-B(x)+\int\limits_{t=0}^{x}\int\limits_{y=0}^{\infty}R((1+a)(x-t)+y)f_{B,C}(t,y)dydt. \end{array} \end{equation}

Assuming that C_i are i.i.d. random variables that follow $exp(\mu)$ (the analysis is still applicable if we consider hyperexponential or Erlang distribution, but to keep the model as simple we can for the shake of readability we assume exponential distribution for the gain sizes), and applying Laplace transforms we come up with the following functional equation:

\begin{equation*} \begin{array}{rl} \rho(s)= &\frac{1-\phi(s)}{s}+ \phi(s)\frac{\mu}{\mu(1+a)-s}(\rho(\frac{s}{1+a})-\rho(\mu))\\ &+\theta\int\limits_{x=0}^{\infty}e^{-sx}\int\limits_{t=0}^{x}\int\limits_{y=0}^{\infty}R((1+a)(x-t)+y)p^{\prime}(F_{B}(t))g^{\prime}(1-e^{-\mu y})f_{B}(t)\mu e^{-\mu y}dydtdx, \end{array} \end{equation*}

where the functions g and h are defined above. In order to rewrite the third term in the RHS of the above equation, let us define the r.v. Z with density given by $g_{Z}(t)=f_{B}(t)-f_{B}(t)p^{\prime}(F_{B}(t))$, and the corresponding Laplace transform $g_{Z}^{*}(s)$. Let also define the function $k_{C}(y)=g^{\prime}(1-e^{-\mu y})\mu e^{-\mu y}$. Thus,

(3.16)\begin{equation} \begin{array}{rl} \rho(s)= &\frac{1-\phi(s)}{s}+ \phi(s)\frac{\mu}{\mu(1+a)-s}(\rho(\frac{s}{1+a})-\rho(\mu))\\ &+\theta\int\limits_{x=0}^{\infty}e^{-sx}\int\limits_{t=0}^{x}\int\limits_{y=0}^{\infty}R((1+a)(x-t)+y)k_{C}(y)[f_{B}(t)-g_{Z}(t)]dydtdx. \end{array} \end{equation}

By the definition of the function g, simple computations imply that

\begin{equation*} \begin{array}{rl} k_{C}(y)= & g^{\prime}(1-e^{-\mu y})\mu e^{-\mu y}=\mu(1-e^{-\mu y})^{c-1}[(c+d)e^{-\mu(d+1)y}-de^{-\mu dy}]\\ =&\mu\sum_{i=0}^{c-1}\left(\!{c-1\atop i}\!\right)(-1)^{i} [(c+d)e^{-\mu(d+c-i)y}-de^{-\mu (c-1-i+d)y}]. \end{array} \end{equation*}

Then, the triple integral in the RHS of (3.16) can be rewritten as follows:

\begin{equation*} \begin{array}{rl} I(s)= &\theta(\phi(s)-g_{Z}^{*}(s))\int\limits_{t=x}^{\infty}e^{-s(x-t)}\int\limits_{z=(1+a)(x-t)}^{\infty}R(z)k_{C}(z-(1+a)(x-t))dzdx \\ =& \theta(\phi(s)-g_{Z}^{*}(s))\mu\sum_{i=0}^{c-1}\left(\!{c-1\atop i}\!\right)(-1)^{i}\left\{(c+d)\int\limits_{w=0}^{\infty}e^{-w(s-\mu(1+a)(c+d-i))}\right.\\ & \times\int\limits_{z=(1+a)w}^{\infty}R(z)e^{-\mu(c+d-i)z}dzdw\\ & \left.-d\int\limits_{w=0}^{\infty}e^{-w(s-\mu(1+a)(c+d-i-1))}\int\limits_{z=(1+a)w}^{\infty}R(z)e^{-\mu(c+d-i-1)z}dzdw\right\} \\ =& \theta(\phi(s)-g_{Z}^{*}(s))\mu\sum_{i=0}^{c-1}\left(\!{c-1\atop i}\!\right)(-1)^{i}\left\{(c+d)\int\limits_{z=0}^{\infty}R(z)e^{-\mu(c+d-i)z}\frac{1-e^{-\frac{z}{1+a}(s-\mu(1+a)(c+d-i))}}{s-\mu(1+a)(c+d-i)}dz\right.\\ & \left.-d\int\limits_{z=0}^{\infty}R(z)e^{-\mu(c+d-i-1)z}\frac{1-e^{-\frac{z}{1+a}(s-\mu(1+a)(c+d-i-1))}}{s-\mu(1+a)(c+d-i-1)}dz\right\}\\ =&\theta(\phi(s)-g_{Z}^{*}(s))\mu\sum_{i=0}^{c-1}\left(\!{c-1\atop i}\!\right)(-1)^{i}\left\{\frac{c+d}{s-\mu(1+a)(c+d-i)}(\rho(\mu(s+d-i))-\rho(\frac{s}{1+a}))\right.\\ &\left.-\frac{d}{s-\mu(1+a)(c+d-i-1)}(\rho(\mu(s+d-i-1))-\rho(\frac{s}{1+a}))\right\}. \end{array} \end{equation*}

Hence, (3.16) is now rewritten as

(3.17)\begin{equation} \rho(s)=J(s;\theta)\rho(\frac{s}{1+a})+H(s;\theta), \end{equation}

where now,

(3.18)\begin{eqnarray} \begin{array}{rl} J(s;\theta):= &\frac{\mu\phi(s)}{\mu(1+a)-s}+\theta(\phi(s)-g_{Z}^{*}(s))\mu\sum_{i=0}^{c-1}\!\left(\!{c-1\atop i}\!\right)(-1)^{i}\!\left[\!\frac{c+d}{\mu(1+a)(c+d-i)-s}-\frac{d}{\mu(1+a)(c+d-i-1)-s}\!\right], \\ H(s;\theta):= & \frac{1-\phi(s)}{s}+\widehat{H}(s;\theta)\\=& \frac{1-\phi(s)}{s}-\frac{\mu\phi(s)\rho(\mu)}{\mu(1+a)-s}\\ &+\theta(\phi(s)-g_{Z}^{*}(s))\mu\sum_{i=0}^{c-1}\left(\!{c-1\atop i}\!\right)(-1)^{i}\left\{\frac{(c+d)\rho(\mu(c+d-i))}{s-\mu(1+a)(c+d-i)}-\frac{d\rho(\mu(c+d-i-1))}{s-\mu(1+a)(c+d-i-1)}\right\}. \end{array}\nonumber\\ \end{eqnarray}

Note that (3.17) has exactly the same form with (3.3), although the functions $J(s;\theta)$, $H(s;\theta)$ are different both because the marginals are now exponentially distributed (when deriving (3.3) we assumed they are Erlang distributed), and especially because we now used the GFGM copula (instead of the standard FGM copula).

Iterating (3.3) N − 1 times results in

\begin{equation*} \rho(s)=\sum_{k=0}^{N-1}\prod_{j=0}^{k-1}J(\frac{s}{(1+a)^{j}} ;\theta)H(\frac{s}{(1+a)^{k}};\theta)+\rho(\frac{s}{(1+a)^{N}})\prod_{j=0}^{N-1}J(\frac{s}{(1+a)^{j}};\theta), \end{equation*}

with an empty product being equal to 1. Observe that for large k, $J(\frac{s}{(1+a)^{k}} ;\theta)$ approaches

\begin{equation*} \begin{array}{l} \frac{\phi(0)}{1+a}+\theta(\phi(0)-g_{Z}^{*}(0))\sum_{i=0}^{c-1}\left(\!{c-1\atop i}\!\right)(-1)^{i}\left[\frac{c+d}{(1+a)(c+d-i)}-\frac{d}{(1+a)(c+d-i-1)}\right]=\frac{\phi(0)}{1+a}=\frac{1}{1+a} \lt 1, \end{array} \end{equation*}

and $H(\frac{s}{(1+a)^{k}};\theta)$ approaches a constant. Thus,

(3.19)\begin{equation} \rho(s)=\sum_{k=0}^{\infty}\prod_{j=0}^{k-1}J(\frac{s}{(1+a)^{j}} ;\theta)H(\frac{s}{(1+a)^{k}};\theta). \end{equation}

The values of $\rho(\mu)$, $\rho(\mu(c+d-i))$, $i=0,\ldots,c$, can be derived by substituting $s=\mu$, $s=\mu(c+d-i)$, $i=0,\ldots,c$, in (3.19), and solving a system of c + 2 equations.

Remark 3.6. It would be interesting to consider the case where with probability p the ith jump has size $au+C_{i}$ (when $U(T_{i}^{-})=u$), and with probability $q:=1-p$, has size $D_{i}\sim\exp(\nu)$, where now the components of $(B_{i},C_{i})$ are dependent through the GFGM with parameters $(\theta_{1},k_{1},b_{1},c_{1},d_{1})$, and the components of $(B_{i},D_{i})$ are dependent through the GFGM with parameters $(\theta_{2},k_{2},b_{2},c_{2},d_{2})$, where $\theta_{1}\neq\theta_{2}, ,k_{1}\neq k_{2},b_{1}\neq b_{2},c_{1}\neq c_{2},d_{1}\neq d_{2}$. Then,

(3.20)\begin{equation} \begin{array}{rl} \rho(s)= &\frac{1-\phi(s)}{s}+ p\left\{\phi(s)\frac{\mu}{\mu(1+a)-s}(\rho(\frac{s}{1+a})-\rho(\mu))\right.\\ &\left.+\theta_{1}\int\limits_{x=0}^{\infty}e^{-sx}\int\limits_{t=0}^{x}\int\limits_{y=0}^{\infty}R((1+a)(x-t)+y)k_{C}(y)[f_{B}(t)-g_{Z}(t)]dydtdx\right\}\\ &+q\left\{\phi(s)\frac{\nu}{\nu-s}(\rho(s)-\rho(\nu))\right.\\ &\left.+\theta_{2}\int\limits_{x=0}^{\infty}e^{-sx}\int\limits_{t=0}^{x}\int\limits_{y=0}^{\infty}R(x-t+y)k_{D}(y)[f_{B}(t)-\tilde{g}_{Z}(t)]dydtdx\right\}, \end{array} \end{equation}

where $\tilde{g}_{Z}(t)=f_{B}(t)-f_{B}(t)p_{*}^{\prime}(F_{B}(t))$, with $p_{*}(u):=u^{k_{2}}(1-u)^{b_{2}}$. Simple but tedious computations yield

(3.21)\begin{equation} \rho(s)=J_{1}(s;\theta_{1},\theta_{2})\rho(\frac{s}{1+a})+H_{1}(s;\theta_{1},\theta_{2}), \end{equation}

where

\begin{equation*} \begin{array}{rl} J_{1}(s;\theta_{1},\theta_{2}):= &\frac{pJ(s;\theta_{1})}{1-qG(s;\theta_{2})}, \\ G(s;\theta_{2}):= & \frac{\nu\phi(s)}{\nu-s}+\theta_{2}(\phi(s)-\tilde{g}_{Z}^{*}(s))\nu\sum_{i=0}^{c_{2}-1}\left(\!{c_{2}-1\atop i}\!\right)(-1)^{i}\left[\frac{c_{2}+d_{2}}{\nu(c_{2}+d_{2}-i)-s}-\frac{d}{\nu(c_{2}+d_{2}-i-1)-s}\right],\\ H_{1}(s;\theta_{1},\theta_{2}):=&\frac{\frac{1-\phi(s)}{s}+p\widehat{H}(s;\theta_{1})+q\tilde{H}(s;\theta_{2})}{1-qG(s;\theta_{2})},\\ \tilde{H}(s;\theta_{2}):=&\frac{\nu\phi(s)\rho(\nu)}{\nu-s}\\ &+\theta_{2}(\phi(s)-\tilde{g}_{Z}^{*}(s))\nu\sum_{i=0}^{c_{2}-1}\left(\!{c_{2}-1\atop i}\!\right)(-1)^{i}\left\{\frac{(c_{2}+d_{2})\rho(\nu(c_{2}+d_{2}-i))}{s-\nu(c_{2}+d_{2}-i)}-\frac{d\rho(\nu(c_{2}+d_{2}-i-1))}{s-\nu(c_{2}+d_{2}-i-1)}\right\}, \end{array} \end{equation*}

and $J(s;\theta_{1})$, $\widehat{H}(s;\theta_{1})$ as given in (3.18), where $\theta_{1}\equiv\theta$, and $c\equiv c_{1}$, $d\equiv d_{1}$, $k\equiv k_{1}$, $b\equiv b_{1}$.

By letting $\nu_{i}:=\nu(d_{2}+i-1)$, $i=1,\ldots,c_{2}+1$, $G(s;\theta_{2})$ is now rewritten as

\begin{equation*}\begin{array}{l} G(s;\theta_{2}):=\frac{\nu\phi(s)\prod_{i=1}^{c_{2}+1}(\nu_{i}-s)+\theta_{2}(\phi(s)-\tilde{g}_{Z}^{*}(s))\nu(\nu-s)\sum_{j=0}^{c_{2}-1}\left(\!{c_{2}-1\atop j}\!\right)(-1)^{j}\Bigl(c(\nu_{c+1-j}-s)-\nu_{c+1}\Bigr)\prod_{l\neq,c-j,c+1-j}(\nu_{l}-s)}{(\nu-s)\prod_{i=1}^{c_{2}+1}(\nu_{i}-s)}.\end{array} \end{equation*}

Note that (3.21) has the same form as the one in (3.17), and its solution will have the form of (3.19), that is,

(3.22)\begin{equation} \rho(s)=\sum_{k=0}^{\infty}\prod_{j=0}^{k-1}J_{1}(\frac{s}{(1+a)^{j}} ;\theta_{1},\theta_{2})H_{1}(\frac{s}{(1+a)^{k}};\theta_{1},\theta_{2}). \end{equation}

However, we still have to obtain $c_{1}+c_{2}+4$ unknown terms, namely, $\rho(\mu)$, $\rho(\mu(d_{1}+i-1))$, $i=1,\ldots,c_{1}+1$, and $\rho(\nu)$, $\rho(\nu(d_{2}+j-1))$, $j=1,\ldots,c_{2}+1$. Clearly, $c_{1}+2$ equations can be derived by simply substituting $s=\mu$, and $s=\mu (d_{1}+i-1)$, $i=1,\ldots,c_{1}+1$ in (3.22). In order to construct other $c_{2}+2$, we need to show that the equation $1-qJ_{2}(s;\theta_{2})=0$ has exactly $c_{2}+2$ roots, say $s_{1},\ldots,s_{c+2}$, such that $Re(s_{k}) \gt 0$, $k=1,\ldots,c_{2}+2$. This task can be accomplished by using Rouché’s theorem [Reference Titchmarsh48] as we did in Proposition 3.5. However, by substituting $s=s_{k}$ in

(3.23)\begin{equation} \rho(s)[1-qG(s;\theta_{2})]=\frac{1-\phi(s)}{s}+p(J(s;\theta_{1})\rho(\frac{s}{1+a})+\widehat{H}(s;\theta_{1}))+q\tilde{H}(s;\theta_{2}), \end{equation}

another unknown is introduced, that is, $\rho(\frac{s_{k}}{1+a})$. However, this unknown can be expressed in terms of $\rho(\mu)$, $\rho(\mu(d_{1}+i-1))$, $i=1,\ldots,c_{1}+1$, and $\rho(\nu)$, $\rho(\nu(d_{2}+j-1))$, $j=1,\ldots,c_{2}+1$, through (3.22).

Proposition 3.7. For $\theta_{2}\neq 0$, the equation $qG(s;\theta_{2})=1$ has exactly $c_{2}+2$ roots, say $s_{1},s_{2},\ldots,s_{c_{2}+2}$ in the right-half complex plane, that is, $Re(s_{j}) \gt 0$, $j=1,2,\ldots,c_{2}+2$.

Proof. The proof is based on Rouché’s theorem [Reference Titchmarsh48] on a contour C_r, consisting of the imaginary axis running from −ir to ir and a semi-circle with radius r running clockwise from ir to −ir. Let $r\to\infty$ and denote by C the limiting contour. It suffices to show that the equation

(3.24)\begin{equation}\begin{array}{l} (\nu-s)\prod_{i=1}^{c_{2}+1}(\nu_{i}-s) =q\nu\left\{\phi(s)\prod_{i=1}^{c_{2}+1}(\nu_{i}-s)\right.\\ \left.+\theta_{2}(\phi(s)-\tilde{g}_{Z}^{*}(s))(\nu-s)\sum_{j=0}^{c_{2}-1}\left(\!{c_{2}-1\atop j}\!\right)(-1)^{j}\Bigl(c(\nu_{c+1-j}-s)-\nu_{c+1}\Bigr)\prod_{l\neq,c-j,c+1-j}(\nu_{l}-s)\right\}, \end{array} \end{equation}

has exactly $c_{2}+2$ roots with positive real parts. Equivalently, we need to show that $|qG(s;\theta_{2})| \lt 1$ on C. Note that $\frac{\nu}{\nu-s}$ and

\begin{equation*} d(s):= \frac{\nu\sum_{j=0}^{c_{2}-1}\left(\!{c_{2}-1\atop j}\!\right)(-1)^{j}\Bigl(c(\nu_{c+1-j}-s)-\nu_{c+1}\Bigr)\prod_{l\neq,c-j,c+1-j}(\nu_{l}-s)}{\prod_{i=1}^{c_{2}+1}(\nu_{i}-s)}, \end{equation*}

are ratios of polynomials with a strictly higher degree at the denominator. Thus, $|qG(s;\theta_{2})|\to 0$ on C (except from $Re(s)=0$). For $Re(s)=0$,

\begin{equation*} \begin{array}{rl} |qG(s;\theta_{2})|\leq &q\left(|\phi(s)\frac{\nu}{\nu-s}|+|\theta_{2}||\phi(s)-\tilde{g}_{Z}(s)||d(s)|\right) \\ \leq &q\left(|\frac{\nu}{\nu-s}|+|\theta_{2}||d(s)|\right)\leq q(1+\theta_{2}|d(0)|)=q \lt 1, \end{array} \end{equation*}

since simple computations imply that $d(0)=0$. Therefore, Rouché’s theorem implies that (3.24) or equivalently $qG(s;\theta_{2})=1$ has exactly $c_{2}+2$ roots, say s_j, $j=1,\ldots,c_{2}+2$ such that $Re(s_{j}) \gt 0$.

3.3. The dual risk model with proportional gains and a linear dependence among gain interarrival times and surplus level

Following the notion of this section, we consider the case where gain interarrivals and gain sizes are dependent based on an FGM copula. On top of that, we further assume that when the surplus level is x the next gain interarrival time is cx, $c\in (0,1)$. Therefore, we assume that the gain interarrival times are linearly dependent on the surplus level. For convenience, we assume that the gain interarrivals are exponentially distributed with rate λ and gain sizes are exponentially distributed with rate µ. Then,

(3.25)\begin{equation} \begin{array}{rl} R(x)=&e^{-\lambda(1-c)x}+\int\limits_{t=cx}^{x}\int\limits_{y=0}^{\infty}R((1+a)(x-t)+y)f_{B,C}(t-cx,y)dydt, \end{array} \end{equation}

where

\begin{equation*} f_{B,C}(t,y)=\lambda\mu e^{-\lambda t}\mu e^{-\mu y}+\theta(2\lambda e^{-2\lambda t}-\lambda e^{-\lambda t})(2\mu e^{-2\mu y}-\mu e^{-\mu y}),\,t,y\geq 0,\theta\in[-1,1]. \end{equation*}

Remark 3.8. Note that for θ = 0, that is, the independent case, and a = 0, our model reduces to the model in [Reference Boxma and Mandjes21, Sect. 3].

Taking Laplace transforms, setting $x-t=w$ and $\bar{c}:=1-c$, we obtain after some algebra:

\begin{equation*} \begin{array}{rl} \rho(s)=&\frac{1}{s+\lambda\bar{c}}+\lambda\mu\int\limits_{x=0}^{\infty}e^{-sx}\int\limits_{w=0}^{\bar{c}x}\int\limits_{y=0}^{\infty}R(w(1+a)+y)\left[(1+\theta) e^{-\lambda (\bar{c}x-w)} e^{-\mu y} +4\theta e^{-2\lambda (\bar{c}x-w)} e^{-2\mu y} \right.\\ &\left.-2\theta e^{-2\lambda (\bar{c}x-w)} e^{-\mu y}-2\theta e^{-\lambda (\bar{c}x-w)} e^{-2\mu y}\right]dydwdx\\ =& \frac{1}{s+\lambda\bar{c}}+\lambda\mu[(1+\theta)I_{1}(s)+4\theta I_{2}(s)-2\theta(I_{3}(s)+I_{4}(s))], \end{array} \end{equation*}

where

\begin{equation*} \begin{array}{rl} I_{1}(s)= & \int\limits_{x=0}^{\infty}e^{-sx}\int\limits_{w=0}^{\bar{c}x}\int\limits_{y=0}^{\infty}R(w(1+a)+y) e^{-\lambda (\bar{c}x-w)} e^{-\mu y}dydwdx \\ =& \int\limits_{x=0}^{\infty}e^{-sx}\int\limits_{w=0}^{\bar{c}x}\int\limits_{z=w(1+a)}^{\infty}R(z) e^{-\lambda (\bar{c}x-w)}\mu e^{-\mu (z-w(1+a))}dydwdx\\ =&\int\limits_{w=0}^{\infty}e^{w(\lambda+\mu(1+a))}\int\limits_{x=w/\bar{c}}^{\infty}e^{-x(s+\lambda\bar{c})}\int\limits_{z=w(1+a)}^{\infty}R(z) e^{-\mu z}dzdxdw\\ =&\frac{1}{s+\lambda\bar{c}}\int\limits_{w=0}^{\infty}e^{-w(\frac{s}{\bar{c}}-\mu(1+a))}\int\limits_{z=w(1+a)}^{\infty}R(z) e^{-\mu z}dzdw\\ =&\frac{\bar{c}}{(s+\lambda\bar{c})(s-\bar{c}\mu(1+a))}(\rho(\mu)-\rho(\frac{s}{\bar{c}(1+a)})). \end{array} \end{equation*}

Similarly,

\begin{equation*} \begin{array}{rl} I_{2}(s)= & \int\limits_{x=0}^{\infty}e^{-sx}\int\limits_{w=0}^{\bar{c}x}\int\limits_{y=0}^{\infty}R(w(1+a)+y) e^{-2\lambda (\bar{c}x-w)} e^{-2\mu y}dydwdx \\ =& \int\limits_{x=0}^{\infty}e^{-sx}\int\limits_{w=0}^{\bar{c}x}\int\limits_{z=w(1+a)}^{\infty}R(z) e^{-2\lambda (\bar{c}x-w)}\mu e^{-2\mu (z-w(1+a))}dydwdx\\ =&\int\limits_{w=0}^{\infty}e^{w(2\lambda+2\mu(1+a))}\int\limits_{x=w/\bar{c}}^{\infty}e^{-x(s+2\lambda\bar{c})}\int\limits_{z=w(1+a)}^{\infty}R(z) e^{-2\mu z}dzdxdw\\ =&\frac{1}{s+2\lambda\bar{c}}\int\limits_{w=0}^{\infty}e^{-w(\frac{s}{\bar{c}}-2\mu(1+a))}\int\limits_{z=w(1+a)}^{\infty}R(z) e^{-2\mu z}dzdw\\ =&\frac{\bar{c}}{(s+2\lambda\bar{c})(s-\bar{c}2\mu(1+a))}(\rho(2\mu)-\rho(\frac{s}{\bar{c}(1+a)})).\\ I_{3}(s)= & \int\limits_{x=0}^{\infty}e^{-sx}\int\limits_{w=0}^{\bar{c}x}\int\limits_{y=0}^{\infty}R(w(1+a)+y) e^{-2\lambda (\bar{c}x-w)} e^{-\mu y}dydwdx\\ =&\frac{\bar{c}}{(s+2\lambda\bar{c})(s-\bar{c}\mu(1+a))}(\rho(\mu)-\rho(\frac{s}{\bar{c}(1+a)})).\\ I_{4}(s)= & \int\limits_{x=0}^{\infty}e^{-sx}\int\limits_{w=0}^{\bar{c}x}\int\limits_{y=0}^{\infty}R(w(1+a)+y) e^{-\lambda (\bar{c}x-w)} e^{-2\mu y}dydwdx\\ =&\frac{\bar{c}}{(s+\lambda\bar{c})(s-\bar{c}2\mu(1+a))}(\rho(2\mu)-\rho(\frac{s}{\bar{c}(1+a)})). \end{array} \end{equation*}

To conclude, the Laplace transform of the ruin probability starting at surplus level x satisfies the following functional equation:

(3.26)\begin{equation} \rho(s)=J(s;\theta)\rho(\zeta(s))+H(s;\theta), \end{equation}

where now, $\zeta(s):=\frac{s}{\bar{c}(1+a)}$, and,

\begin{equation*} \begin{array}{rl} J(s;\theta)=&\frac{\lambda\mu\bar{c}[(s+2\lambda\bar{c})(2\mu\bar{c}(1+a)-s)-\theta s^{2}]}{(s+2\lambda\bar{c})(2\mu\bar{c}(1+a)-s)(s+\lambda\bar{c})(\mu\bar{c}(1+a)-s)}, \\ H(s;\theta)=&\frac{1}{s+\lambda\bar{c}}\left(1-\lambda\mu\bar{c}\left[\frac{s(1-\theta)+2\lambda\bar{c}}{(s+2\lambda\bar{c})(\mu\bar{c}(1+a)-s)}\rho(\mu)+\frac{2\theta s}{(s+2\lambda\bar{c})(2\mu\bar{c}(1+a)-s)}\rho(2\mu)\right]\right). \end{array} \end{equation*}

Note that (3.26) has the same form as the one in [Reference Boxma, Frostig and Palmowski18, Eq. (2.1)]. Iteration of (3.26) results in the following theorem.

Theorem 3.9. The Laplace transform of the ruin probability $\rho(s)$ is given by

(3.27)\begin{equation} \rho(s)=\sum_{k=0}^{\infty}\prod_{j=0}^{k-1}J(\zeta^{(j)}(s);\theta)H(\zeta^{(k)}(s);\theta), \end{equation}

where $\zeta^{(j)}(s):=\zeta(\zeta^{(j-1)}(s))=\frac{s}{(\bar{c}(1+a))^{j}}$, with $\zeta^{(0)}(s)=s$. The unknowns $\rho(\mu)$, $\rho(2\mu)$ (that appear with a prefactor in all $H(\zeta^{(k)}(s);\theta)$) are derived by the system of equations that is constructed by substituting $s=\mu$, $s=2\mu$ in (3.27).

3.4. Numerical example

In the following, we cope with the derivation of some numerical results based on the theoretical findings in Theorem 3.9. In particular, we provide values of the ruin probabilities as functions of the surplus level and investigate the impact of the parameters. All numbers are generated by using the software package Mathematica (Note that no actual time unit is specified here).

Numerical results shown that we need a (relatively) small finite number of terms (depending on the values of the parameters) in order to obtain the values of $\rho(\mu)$, $\rho(2\mu)$ in the expression (3.27) ( $\rho(\mu)$, $\rho(2\mu)$ are unknowns in $H(s;\theta)$ and are derived as stated in Theorem 3.9), as well as to derive the Laplace transform $\rho(s)$. In doing that, we truncate the infinite sum at a specific value of k, such that the absolute difference of two consecutive terms to be smaller than 10⁻⁷, that is, if $u_{1}(k):=\prod_{j=0}^{k-1}J(\zeta^{(j)}(\mu);\theta)$, $u_{2}(k):=\prod_{j=0}^{k-1}J(\zeta^{(j)}(2\mu);\theta)$, find k such that $|u_{m}(k)-u_{m}(k+1)|\leq 10^{-7}$, $m=1,2$. We observed that as we increase a the number of terms needed in the summation in (3.27) decrease very fast (see also Table 1). The ruin probabilities R(x) are obtained by using Euler algorithm described in [Reference Abate and Whitt2, pp. 7–8] to numerically invert $\rho(s)$, that is, $R(x)=\mathcal{L}^{-1}(\rho(s))$. To obtain $\rho(s)$, we need even fewer non-negligible terms in the infinite sum (3.27). In Figure 1, we can see two instances of $\rho(s)$ when we truncate the infinite sum in (40) at k = 2, and k = 30 (when λ = 1, µ = 4, θ = 0.7, a = 3). It is obvious that the two instances of $\rho(s)$ are almost identical. In the numerical examples, we derive below we truncate the infinite sum in (3.27) at k = 30.

In Table 1, we observe the effect of the proportionality parameter a on the ruin probabilities R(x) for given values of the initial surplus x. We observe that as we increase a, R(x) decreases as expected since the surplus level also increases. Note that by increasing a, the number of non-negligible terms in the infinite sum in (3.27) are decreasing. For example, when a = 3, $u_{1}(11)=5.76757*10^{-8}$, with $|u_{1}(11)-u_{1}(12)|=4.16547*10^{-8}$, $u_{2}(10)=7.85483*10^{-8}$, with $|u_{2}(10)-u_{2}(11)|=5.67299*10^{-8}$, thus, $u_{m}(k)\approx 0$ for $k\geq 11$, $m=1,2$. When a = 9, $u_{1}(7)=2.68044*10^{-8}$ with $|u_{1}(7)-u_{1}(8)|=2.38262*10^{-8}$, and $u_{2}(6)=1.04183*10^{-7}$, $|u_{2}(6)-u_{2}(7)|=9.26076*10^{-8}$, thus in that case, $u_{m}(k)\approx 0$ for $k\geq 7$, $m=1,2$. So, the more we increase a, the less non-negligible terms we need. Moreover, when we start with large initial surplus, R(x) decreases, and that decrease becomes more apparent as a increases.

Table 1. The effect of proportionality parameter a on ruin probabilities R(x) when λ = 1, µ = 4, c = 0.1, θ = 0.7.

Figure 2 shows the way the ruin probabilities R(x) are affected when we increase c, for varying values of the initial surplus level. In particular, we can see that as c increases, that is, when the gain interarrival increases proportionally to the surplus level x, the ruin probability also increases as expected, since the time until the next gain interarrival also increases. Moreover, that increase becomes more apparent when the surplus level is closer to zero.

Figure 2. The effect of c on the ruin probabilities for λ = 3, µ = 4, θ = 0.7, a = 6.

In Figure 3, we observe the effect of the dependence parameter θ on $\rho(s)$. As we can see from Figure 3, the dependence parameter θ has a clear impact on the values of $\rho(s)$. In particular, the higher the dependence parameter the higher the value of $\rho(s)$ is. Equivalently, when the dependence relation is positive, the probability of having an important gain increases as the time elapsed since the last gain increases. As a result, the ruin probability will be lower.

Figure 3. The effect of θ on $\rho(s)$ for λ = 8, µ = 4, c = 0.1, a = 6.

To conclude, despite the fact that the explicit expressions of the Laplace transforms are given in terms of infinite sums of products, we can easily derive numerical results due to the fact that numerically, we need only few terms (a finite number) in these sums, thus, making our results fully numerically tractable.

The next step is to focus on the Laplace–Stieltjes transform of the ruin time τ_x, that is, $T(s,x)$:

\begin{equation*}\begin{array}{rl} T(s,x)=&e^{-sx}e^{-\lambda(1-c)x}+\int\limits_{t=cx}^{x}e^{-st}\int\limits_{y=0}^{\infty}T(s,(1+a)(x-t)+y)f_{B,C}(t-cx,y)dydt\\ =& e^{-x(s+\lambda\bar{c})}+\int\limits_{w=0}^{\bar{c}x}e^{-s(x-w)}\int\limits_{y=0}^{\infty}T(s,(1+a)w+y)f_{B,C}(\bar{c}x-w,y)dydw. \end{array} \end{equation*}

Let $\tau(s,\beta)$ the Laplace transform of the ruin time Laplace–Stieltjes transform $T(s,x)$. Then,

\begin{equation*} \begin{array}{rl} \tau(s,\beta)= &\frac{1}{s+\lambda\bar{c}+\beta}+\int\limits_{x=0}^{\infty}e^{-s\beta}\int\limits_{w=0}^{\bar{c}x}e^{-s(x-w)}\int\limits_{z=(1+a)w}^{\infty}R(z)f_{B,C}(\bar{c}x-w,z-w(1+a))dzdwdx. \end{array} \end{equation*}

By repeating similar steps as above, we obtain after some algebra the following functional equation:

(3.28)\begin{equation} \tau(s,\beta)=\tilde{J}(s,\beta;\theta)\tau(s,\psi(s,\beta))+\tilde{H}(s,\beta;\theta), \end{equation}

where now, $\psi(s,\beta):=\frac{r(s,\beta)}{1+a}=\frac{\beta+cs}{\bar{c}(1+a)}$, and,

\begin{equation*} \begin{array}{rl} \tilde{J}(s,\beta;\theta)=&\frac{\lambda\mu[(s+2\lambda\bar{c}+\beta)(2\mu(1+a)-r(s,\beta))-\theta (s+\beta)r(s,\beta)]}{(s+2\lambda\bar{c}+\beta)(2\mu(1+a)-r(s,\beta))(s+\lambda\bar{c}+\beta)(\mu(1+a)-r(s,\beta))}, \\ \tilde{H}(s,\beta;\theta)=&\frac{1}{s+\lambda\bar{c}+\beta}\left(1-\lambda\mu\left[\frac{(s+\beta)(1-\theta)+2\lambda\bar{c}}{(s+2\lambda\bar{c}+\beta)(\mu(1+a)-r(s,\beta))}\tau(s,\mu)+\frac{2\theta (s+\beta)}{(s+2\lambda\bar{c}+\beta)(2\mu(1+a)-r(s,\beta))}\tau(s,2\mu)\right]\right). \end{array} \end{equation*}

Writing $\psi^{(j)}(s,\beta):=\psi(\psi^{(j-1)}(s,\beta))$ with $\psi^{(0)}(s,\beta):=\beta$, we have,

\begin{equation*} \psi^{(j)}(s,\beta)=\frac{\beta}{(\bar{c}(1+a))^{j}}+\sum_{i=1}^{j}\frac{cs}{(\bar{c}(1+a))^{i}},\,j=1,2,\ldots. \end{equation*}

Iterating (3.28), we obtain the following result.

Theorem 3.10. The Laplace transform of the ruin time Laplace–Stieltjes transform $T(s,x)$ is given by

(3.29)\begin{equation} \tau(s,\beta)=\sum_{k=0}^{\infty}\prod_{j=0}^{k-1}\tilde{J}(\psi^{(j)}(s,\beta);\theta)\tilde{H}(\psi^{(k)}(s,\beta);\theta). \end{equation}

The unknowns $\tau(s,\mu)$, $\tau(s,2\mu)$ (that appear with a prefactor in all $\tilde{H}(\psi^{(k)}(s,\beta;\theta)$) are derived by the system of equations that is constructed by substituting $s=\mu$, $s=2\mu$ in (3.29).

Note that the sums of products in (3.27) and (3.29) converge since both $J(.,.)$, $H(.,.)$ (resp. $\tilde{J}(.,.)$, $\tilde{H}(.,.)$) decrease geometrically fast.

Remark 3.11. Note that for c = 0, a = 0, we cope with the standard (no proportional gains) dual risk model with $\exp(\mu)$ initial capital, and for which there is a dependence among gain interarrivals and gain sizes based on the FGM copula. In queueing terms, we are dealing with an M/M/1 queue and studying the busy period starting from an empty system with an $\exp(\mu)$ upward jump, in which the interarrival times and work that enters the system are dependent based on the FGM copula. Note that $\mu\tau(s,\mu)=\int_{0}^{\infty}T(s,x)\mu e^{-\mu x}dx$ should be equal to the to the Laplace–Stieltjes transform of an M/M/1 queue with initial service time to be $\exp(\mu)$, and for which we consider the dependence structure mentioned above. In particular, for $c=a=0$, (3.28) is rewritten as:

(3.30)\begin{equation} \tau(s,\beta)(1-\widehat{J}(s,\beta;\theta))=\widehat{H}(s,\beta;\theta), \end{equation}

where

\begin{equation*} \begin{array}{rl} \widehat{J}(s,\beta;\theta)= &\lambda\mu\frac{(s+2\lambda+\beta)(2\mu-\beta)-\theta(s+\beta)\beta}{(s+\lambda+\beta)(s+2\lambda+\beta)(\mu-\beta)(2\mu-\beta)}, \\ \widehat{H}(s,\beta;\theta)= &\frac{1}{s+\lambda+\beta}\left(1-\lambda\mu\left[\frac{(s+\beta)(1-\theta)+2\lambda}{(s+2\lambda+\beta)(\mu-\beta)}\tau(s,\mu)+\frac{2\theta (s+\beta)}{(s+2\lambda+\beta)(2\mu-\beta)}\tau(s,2\mu)\right]\right). \end{array} \end{equation*}

It is readily seen that for θ = 0, (3.30) reduces to the one in [Reference Boxma and Mandjes21, Eq. (29)]. Moreover, taking $\beta=\mu$, and $\beta=2\mu$, (3.30) gives an identity. We now investigate the roots of the equation $1-\widehat{J}(s,\beta;\theta)=0$.

Proposition 3.12. For $Re(s) \gt 0$, θ ≠ 0, the equation $1-\widehat{J}(s,\beta;\theta)=0$ has exactly two roots, say $u_{1}(s),u_{2}(s)$ with $Re(u_{i}(s)) \gt 0$, $i=1,2$.

Proof. Following Rouché’s theorem [Reference Titchmarsh48], we consider a contour C_r, consisting of the imaginary axis running from −ir to ir and a semi-circle with radius r running clockwise from ir to −ir. Let $r\to\infty$ and denote by C the limiting contour. Note that $1-\widehat{J}(s,\beta;\theta)=0$ is rewritten as

(3.31)\begin{equation} \lambda\mu[(s+2\lambda+\beta)(2\mu-\beta)-\theta(s+\beta)\beta]=(s+\lambda+\beta)(s+2\lambda+\beta)(\mu-\beta)(2\mu-\beta). \end{equation}

So it suffices to show that (3.31) has exactly two roots in the right-half complex plane. It is readily seen that $|\widehat{J}(s,\beta;\theta)|\to 0$ on C (excluding $Re(\beta)=0$), since it is the sum of ratios of polynomials with a strictly higher degree at the denominator. For $Re(\beta)=0$,

\begin{equation*} \begin{array}{rl} |\widehat{J}(s,\beta;\theta)|= &\lambda\mu|\frac{(s+2\lambda+\beta)(2\mu-\beta)-\theta(s+\beta)\beta}{(s+\lambda+\beta)(s+2\lambda+\beta)(\mu-\beta)(2\mu-\beta)}|\leq \lambda\mu|\frac{1}{(s+\lambda+\beta)(\mu-\beta)}|+\lambda\mu|\theta|\frac{(s+\beta)\beta}{(s+\lambda+\beta)(s+2\lambda+\beta)(\mu-\beta)(2\mu-\beta)}|\\ &\leq |\frac{s}{s+\lambda}| \lt 1. \end{array} \end{equation*}

Thus in any case $|\widehat{J}(s,\beta;\theta)| \lt 1$ or equivalently

\begin{equation*} \lambda\mu|(s+2\lambda+\beta)(2\mu-\beta)-\theta(s+\beta)\beta| \lt |(s+\lambda+\beta)(s+2\lambda+\beta)(\mu-\beta)(2\mu-\beta)|, \end{equation*}

and thus by Rouché’s theorem, it follows that (3.31) has the same number of roots as $(s+\lambda+\beta)(s+2\lambda+\beta)(\mu-\beta)(2\mu-\beta)=0$ inside C_r. Since the latter equation has exactly two positive roots inside C_r, we deduce that (3.31) or equivalently $1-\widehat{J}(s,\beta;\theta)=0$ has exactly two roots, say $u_{1}(s),u_{2}(s)$ with positive real parts. Finally, we complete the proof by letting $r\to\infty$.

Substituting in (3.30) $\beta=u_{1}(s)$, and $\beta=u_{2}(s)$ we construct a system of two equations with unknowns $\tau(s,\mu)$, $\tau(s,2\mu)$. By solving this system, we are able to compute $\mu\tau(s,\mu)$ in terms of the $u_{1}(s),u_{2}(s)$. Note that by setting θ = 0 (i.e., the standard M/M/1 queue) $\mu\tau(s,\mu)$ coincides with the expression in [Reference Boxma and Mandjes21, Eq. (29)].

4. A dual risk model with upward and downward jumps and randomly proportional gains

We now focused on the case of two-sided jumps, that is, the upward and downward jumps can be interpreted as company random gains and random losses respectively. Thus, our model is suitable for insurance companies with business in both property and casualty insurance and life annuities. More precisely, we assume that with probability p gains C_i, and with probability q losses $-D_{i}$ ( $i = 1, 2,\ldots $) that arrive according to a renewal process with general interarrival times. Let, C_is and D_is are exponentially distributed (we also consider the case where they follow Erlang distribution).

On top of that, we add the randomly proportional gain feature. More precisely, if the surplus process just before the ith arrival is at level u, then, the capital jumps to:

(4.1)\begin{equation} \left\{\begin{array}{ll} (1+a_{l})u+C_{i},&\text{with probability }p\times k_{l},\,l=1,\ldots,K, \\ \left[(1+\beta_{h})u-D_{i}\right]^{+},& \text{with probability }q\times m_{h},\,h=1,\ldots,L, \end{array}\right. \end{equation}

where $p+q=1$, and $\sum_{l=1}^{L}=1$, $\sum_{h=1}^{M}m_{h}=1$. Thus, with probability p (resp. q), we have an upward (resp. downward) jump of size C_i (resp. $-D_{i}$), with an additional inflow proportional to u, and equal to $a_{l}u$ (resp. $\beta_{h}u$) with probability k_l, $l=1,\ldots,K$ (resp. m_h, $h=1,\ldots,L$). So the type of the jump affects also the values of proportionality coefficient. Note that for p = 1, $k_{1}=1$, our model reduces to the one in [Reference Boxma, Frostig and Palmowski18]. Then, having in mind that $$[(1+\beta_{h})u-D_{i}]^{+}=\left\{\begin{array}{ll} (1+\beta_{h})u-D_{i},&(1+\beta_{h})u\geq D_{i}, \\ 0,&(1+\beta_{h})u \lt D_{i}, \end{array}\right.$$ we have,

(4.2)\begin{eqnarray} \begin{array}{rl} R(x)=&1-B(x)+p\sum_{l=1}^{K}k_{l}\int\limits_{t=0}^{x}\int\limits_{y=0}^{\infty}R((1+a_{l})(x-t)+y)\mu e^{-\mu y}dydB(t)\\ &+\,q\sum_{h=1}^{K}m_{h}\int\limits_{t=0}^{x}\left\{\int\limits_{y=0}^{(1+\beta_{h})(x-t)}R((1+\beta_{h})(x-t)-y)\nu e^{-\nu y}dy+\int\limits_{y=(1+\beta_{h})(x-t)}^{\infty}\nu e^{-\nu y}dy\right\}\\ &\times\, dB(t). \end{array}\nonumber\\ \end{eqnarray}

Applying Laplace transforms, we come up with the following functional equation,

(4.3)\begin{equation} \rho(s)=\frac{1-\phi(s)}{s}+\sum_{l=1}^{K}p_{l}I_{l}(s)+\sum_{h=1}^{L}q_{h}\widehat{I}_{h}(s), \end{equation}

where $p_{l}=p\times k_{l}$, $l=1,\ldots,K$, and $q_{h}=q\times m_{h}$, $h=1,\ldots,L$. After lengthy, but straightforward computations,

(4.4)\begin{equation} \begin{array}{rl} I_{l}(s)= &\frac{\phi(s)\mu}{s-\mu(1+a_{l})}(\rho(\mu)-\rho(\frac{s}{1+a_{l}})),\,l=1,\ldots,K, \\ \widehat{I}_{h}(s)=&\frac{\phi(s)\nu}{s+\nu(1+\beta_{h})}(\rho(\frac{s}{1+\beta_{h}})+\frac{1}{\nu}),\,h=1,\ldots,L. \end{array} \end{equation}

Substituting (4.4) in (4.3) we come up with the following functional equation:

(4.5)\begin{equation} \rho(s)=\phi(s)\sum_{l=1}^{K+L}g_{l}f_{l}(\sigma_{l}(s))\rho(\zeta_{l}(s))+L(s), \end{equation}

where

\begin{equation*} g_{l}=\left\{\begin{array}{ll} \frac{p_{l}}{1+a_{l}}, &\,l=1,\ldots,K, \\ \frac{q_{l-K}}{1+\beta_{l-K}}, &\,l=K+1,\ldots,K+L, \end{array}\right. \end{equation*}

\begin{equation*} f_{l}(s)=\left\{\begin{array}{ll} \frac{\mu}{\mu+s},&l=1,\ldots,K, \\ \frac{\nu}{\nu+s},&l=K+1,\ldots,K+L, \end{array}\right. \end{equation*}

\begin{equation*} \zeta_{l}(s)=\left\{\begin{array}{ll} a_{l}(s)=\frac{s}{1+a_{l}},&l=1,\ldots,K, \\ \widehat{a}_{l-K}(s)=\frac{s}{1+\beta_{l-K}},&l=K+1,\ldots,K+L, \end{array}\right. \end{equation*}

\begin{equation*} \sigma_{l}(s)=\left\{\begin{array}{ll} -a_{l}(s),&l=1,\ldots,K, \\ \widehat{a}_{l-K}(s),&l=K+1,\ldots,K+L, \end{array}\right. \end{equation*}

and

\begin{equation*} L(s)=\frac{1-\phi(s)}{s}+\phi(s)\sum_{l=1}^{K+L}\widehat{g}_{l}f_{l}(\sigma_{l}(s)), \end{equation*}

with $$\widehat{g}_{l}=\left\{\begin{array}{ll} -g_{l}\rho(\mu),&l=1,\ldots,K, \\ \frac{g_{l}}{\nu},&l=K+1,\ldots,K+L. \end{array}\right.$$

Our aim now is to solve (4.5). Note that the form of (4.5) is similar to [Reference Adan, Boxma and Resing3, Eq. (2)], and thus, a similar approach can be employed to solve it. After N − 1 iterations,

\begin{equation*} \begin{array}{rl} \rho(s)=& \sum_{k=0}^{N-1}\sum_{i_{1}+\ldots+i_{K+L}=k}g_{1}^{i_{1}}\ldots g_{K+L}^{i_{K+L}}G_{i_{1},\ldots,i_{K+L}}(s)L(\zeta_{i_{1},\ldots,i_{K+L}}(s)) \\ &+ \sum_{i_{1}+\ldots+i_{K+L}=k}g_{1}^{i_{1}}\ldots g_{K+L}^{i_{K+L}}G_{i_{1},\ldots,i_{K+L}}(s)\rho(\zeta_{i_{1},\ldots,i_{K+L}}(s)), \end{array} \end{equation*}

where $\zeta_{i_{1},\ldots,i_{K+L}}(s):=\zeta_{1}^{i_{1}}(\zeta_{2}^{i_{2}}(\ldots(\zeta_{K+L}^{i_{K+L}}(s))\ldots))$, and $\zeta_{k}^{m}(s)$ is the mth iterate of $\zeta_{k}(s)$, and the functions $G_{i_{1},\ldots,i_{K+L}}(s)$ are recursively obtained as follows:

\begin{equation*} G_{i_{1},\ldots,i_{K+L}}(s)=\sum_{k=1}^{K+L}G_{i_{1},\ldots,i_{k}-1,\ldots,i_{K+L}}(s)L_{0,\ldots,1,\ldots,0}(\zeta_{i_{1},\ldots,i_{k}-1,\ldots,i_{K+L}}(s)), \end{equation*}

with $G_{0,\ldots,0}(s)=1$, $G_{0,\ldots,1,\ldots,0}(s)=L_{0,\ldots,1,\ldots,0}(s)=\phi(s)f_{l}(\sigma_{l}(s))$, where 1 is in the lth position and $l=1,\ldots,K+L$, and $G_{i_{1},\ldots,i_{K+L}}(s)=0$ if one of the indices equals −1. Then, letting $N\to\infty$,

(4.6)\begin{equation}\begin{array}{rl} \rho(s)=&\sum_{k=0}^{\infty}\sum_{i_{1}+\cdots+i_{K+L}=k}g_{1}^{i_{1}}\cdots g_{K+L}^{i_{K+L}}G_{i_{1},\ldots,i_{K+L}}(s)L(\zeta_{i_{1},\ldots,i_{K+L}}(s))\\ &+ \lim_{N\to\infty}\sum_{i_{1}+\cdots+i_{K+L}=N}g_{1}^{i_{1}}\cdots g_{K+L}^{i_{K+L}}G_{i_{1},\ldots,i_{K+L}}(s), \end{array} \end{equation}

since $\rho(\zeta_{i_{1},\ldots,i_{K+L}}(s))\to\rho(0)=1$ since $\zeta_{l}(s)$, $l=1,\ldots,K+L$ are commutative contraction mappings on the closed positive half plane. It is readily seen that $G_{i_{1},\ldots,i_{K+L}}(s)$ can be written as a finite sum of products. As the number of iterations increases, each of these products vanish. Moreover, as the number of iterations increases $L(\zeta_{i_{1},\ldots,i_{K+L}}(s))$ approaches some constant. Thus,

(4.7)\begin{equation}\begin{array}{rl} \rho(s)=&\sum_{k=0}^{\infty}\sum_{i_{1}+\cdots+i_{K+L}=k}g_{1}^{i_{1}}\cdots g_{K+L}^{i_{K+L}}G_{i_{1},\ldots,i_{K+L}}(s)L(\zeta_{i_{1},\ldots,i_{K+L}}(s)).\end{array} \end{equation}

We still need to derive $\rho(\mu)$. This can be achieved by substituting $s=\mu$ in (4.7).

Remark 4.1. Consider the simpler case where $K=1=L$, and $a_{1}=\beta_{1}\equiv a$, so that

(4.8)\begin{equation} \left\{\begin{array}{ll} (1+a)u+C_{i},&\text{with probability }p, \\ \left[(1+a)u-D_{i}\right]^{+},& \text{with probability }q. \end{array}\right. \end{equation}

Then, (4.5) reduces to

(4.9)\begin{equation} \rho(s)=\phi(s)[p\frac{\mu}{\mu(1+a)-s}+q\frac{\nu}{\nu(1+a)+s}]\rho(\frac{s}{1+a})+\frac{1-\phi(s)}{s}+\phi(s)[\frac{q}{\nu(1+a)+s}-\frac{p\mu}{\mu(1+a)-s}\rho(\mu)]. \end{equation}

Note that for q = 0 (so that p = 1) we recover the model in [Reference Boxma, Frostig and Palmowski18]. Equation (4.9) has the same form as in [Reference Boxma, Frostig and Palmowski18, Eq. (2.5)], where now

\begin{equation*} \begin{array}{rl} J(s)= & p\frac{\mu}{\mu(1+a)-s}+q\frac{\nu}{\nu(1+a)+s}, \\ H(s)=& \frac{1-\phi(s)}{s}+\phi(s)[\frac{q}{\nu(1+a)+s}-\frac{p\mu}{\mu(1+a)-s}\rho(\mu)], \end{array} \end{equation*}

and its solution is as given in [Reference Boxma, Frostig and Palmowski18, Eq. (2.7)], with $J(.)$, $H(.)$, as given above.

In case where $a_{1}\equiv a\neq \beta_{1}\equiv \beta$, that is, when (4.8) is given by

(4.10)\begin{equation} \left\{\begin{array}{ll} (1+a)u+C_{i},&\text{with probability }p, \\ \left[(1+\beta)u-D_{i}\right]^{+},& \text{with probability }q, \end{array}\right. \end{equation}

and we now have to solve

(4.11)\begin{equation} \rho(s)=\phi(s)J_{0}(s)\rho(\frac{s}{1+a})+\phi(s)J_{1}(s)\rho(\frac{s}{1+\beta})+H_{1}(s), \end{equation}

where

\begin{equation*} \begin{array}{rl} J_{0}(s)= & p\frac{\mu}{\mu(1+a)-s},\\ J_{1}(s)=&q\frac{\nu}{\nu(1+\beta)+s}, \\ H_{1}(s)=& \frac{1-\phi(s)}{s}+\phi(s)[\frac{q}{\nu(1+\beta)+s}-\frac{p\mu}{\mu(1+a)-s}\rho(\mu)]. \end{array} \end{equation*}

The form of (4.11) is the same as the one in (2.4), and thus, its solution can be derived similarly.

Clearly, the dual risk model with proportional gains and two-sided jumps (i.e., either an upward or a downward jump) can be analyzed similarly, when we consider the FGM copula to describe the dependence among the gain interarrival time and the corresponding size. Indeed, consider the simpler case where $C_{i}\sim \exp(\mu)$, $D_{i}\sim\exp(\nu)$ (the analysis is still applicable when we consider hyperexponential, or Erlang distributions), but there is a dependence with the “gain” interarrival time based on the FGM copula, that is, the joint density function of the random vectors $(B_{i},C_{i})$, $(B_{i},D_{i})$ is as follows:

\begin{equation*} \begin{array}{rl} f_{B,C}(t,y) = &f_{B}(t)\mu e^{-\mu y}+\theta_{1}h(t)(2\mu e^{-2\mu y}-\mu e^{-\mu y}),\,\theta_{1}\in[-1,1], \\ f_{B,D}(t,y) = &f_{B}(t)\nu e^{-\nu y}+\theta_{2}h(t)(2\nu e^{-2\nu y}-\nu e^{-\nu y}),\,\theta_{2}\in[-1,1], \end{array} \end{equation*}

where $h(t)=f_{B}(t)(1-2F_{B}(t))$; see also Section 3. Note that asking $\theta_{1}=0$ and/or $\theta_{2}=0$, we come up with the independent case.

Then, we have

(4.12)\begin{equation} \begin{array}{rl} R(x)=&1-B(x)+p\int\limits_{t=0}^{x}\int\limits_{y=0}^{\infty}R((1+a)(x-t)+y)f_{B,C}(t,y)dydt\\ &+q\int\limits_{t=0}^{x}\int\limits_{y=0}^{(1+\beta)(x-t)}R((1+\beta)(x-t)-y) f_{B,D}(t,y)dydt+q\int\limits_{y=(1+\beta)(x-t)}^{\infty}f_{B,C}(t,y)dydt. \end{array} \end{equation}

Taking Laplace transforms, we come up with the functional equation

(4.13)\begin{equation} \rho(s)=pJ_{1}(s;\theta_{1})\rho(\frac{s}{1+a})+qJ_{2}(s;\theta_{2})\rho(\frac{s}{1+\beta})+H(s;\theta_{1},\theta_{2}), \end{equation}

where now,

\begin{equation*} \begin{array}{rl} J_{1}(s;\theta_{1}):= &\mu\frac{\phi(s)-\theta_{1}h^{*}(s)}{\mu(1+a)-s}+\frac{2\mu\theta_{1}h^{*}(s)}{2\mu(1+a)-s}, \\ J_{2}(s;\theta_{2}):= &\nu\frac{\phi(s)-\theta_{2}h^{*}(s)}{\nu(1+\beta)+s}+\frac{2\nu\theta_{2}h^{*}(s)}{2\nu(1+\beta)+s}, \\ H(s;\theta_{1},\theta_{2}):=&\frac{1-\phi(s)}{s}+q[\frac{\phi(s)-\theta_{2}h^{*}(s)}{\nu(1+\beta)+s}+\frac{\theta_{2}h^{*}(s)}{2\nu(1+\beta)+s}]-p[\mu\rho(\mu)\frac{\phi(s)-\theta_{1}h^{*}(s)}{\mu(1+a)-s}+\rho(2\mu)\frac{2\mu\theta_{1}h^{*}(s)}{2\mu(1+a)-s}]. \end{array} \end{equation*}

Note that (4.12) has the same form as the one in (2.4), so we can apply the same method to solve it. After n − 1 ( $n\geq 1$) iterations, we have

\begin{equation*} \rho(s)=\sum_{k=0}^{n}\sum_{i_{1}+i_{2}=k}p^{i_{1}}q^{i_{2}}L_{i_{1},i_{2}}(s)H(f_{i_{1},i_{2}}(s);\theta_{1},\theta_{2})+\sum_{i_{1}+i_{2}=n}p^{i_{1}}q^{i_{2}}L_{i_{1},i_{2}}(s)\rho(f_{i_{1},i_{2}}(s)), \end{equation*}

where $L_{i_{1},i_{2}}(s)$ are recursively derived as in (2.6), and $f_{i_{1},i_{2}}(s)=a_{i_{1}}(b_{i_{2}}(s))=b_{i_{2}}(a_{i_{1}}(s))=\frac{s}{(1+a)^{i_{1}}(1+\beta)^{i_{2}}}$ with $a_{i}(s)=\frac{s}{(1+a)^{i}}$, $b_{j}(s)=\frac{s}{(1+\beta)^{j}}$. Note that $a_{i}(s)$, $b_{j}(s)$ are commutative contraction mappings on the right half plane, and the methodology presented in [Reference Adan, Boxma and Resing3, Sect. 2] applies. In particular,

(4.14)\begin{equation} \rho(s)=\sum_{k=0}^{\infty}\sum_{i_{1}+i_{2}=k}p^{i_{1}}q^{i_{2}}L_{i_{1},i_{2}}(s)H(f_{i_{1},i_{2}}(s);\theta_{1},\theta_{2})+\lim_{n\to\infty}\sum_{i_{1}+i_{2}=n}p^{i_{1}}q^{i_{2}}L_{i_{1},i_{2}}(s). \end{equation}

We still need to obtain $\rho(\mu)$, $\rho(2\mu)$. This can be done by solving a system of two equations that are constructed by setting $s=\mu$, and $s=2\mu$ in (4.14).

5. The dual risk model with uniformly proportional gains

Consider now the case where the proportional parameter is a random variable, that is, the surplus process U(t) with $U(0)=x \gt 0$ evolves as

\begin{equation*} U(t)=x-t+\sum_{i=1}^{N(t)}(C_{i}+V_{i}U(S_{i}^{-})),\,t\geq 0, \end{equation*}

where $V_{1},V_{2},\ldots$ are i.i.d. uniformly distributed random variables on $[a,b]$, $0 \lt a \lt b$. Moreover, N(t) denotes the number of gains that arrive in $(0,t]$, with $S_{i+1}-S_{i}$ are i.i.d. gain interarrival times having c. distribution function $B(.)$, density $b(.)$, and Laplace–Stieltjes transform $\phi(.)$. By assuming the C_is are i.i.d. exponentially distributed random variables with rate µ, the ruin probability R(x) when starting in x satisfies the following equation:

(5.1)\begin{equation} \begin{array}{rl} R(x)=&1-B(x)+\int\limits_{v=a}^{b}\int\limits_{t=0}^{x}\int\limits_{y=0}^{\infty}R((1+v)(x-t)+y)\mu e^{-\mu y}dB(t)\frac{dv}{b-a}. \end{array} \end{equation}

Applying Laplace transform in (5.1) yields

\begin{equation*} \rho(s)=\frac{1-\phi(s)}{s}+I^{*}(s), \end{equation*}

where

\begin{equation*} \begin{array}{rl} I^{*}(s)=&\frac{\phi(s)}{b-a}\int\limits_{v=a}^{b}\int\limits_{x=t}^{\infty}e^{-(x-t)(s-\mu(1+v))}\int\limits_{z=(1+v)(x-t)}^{\infty}\mu e^{-\mu z}R(z)dzdxdv \\ =&\frac{\phi(s)}{b-a} \int\limits_{v=a}^{b}\int\limits_{z=0}^{\infty}\mu e^{-\mu z}R(z)\frac{e^{-\frac{z}{1+v}(s-\mu(1+v))}-1}{\mu(1+v)-s}dzdv \\ =&\frac{\phi(s)}{b-a} \int\limits_{v=a}^{b}\frac{\mu}{\mu(1+v)-s}(\rho(\frac{s}{1+v})-\rho(\mu))dv \\ =&\frac{\phi(s)}{b-a} \int\limits_{v=a}^{b}\frac{\mu}{\mu(1+v)-s}\rho(\frac{s}{1+v})dv-\frac{\phi(s)}{b-a}\rho(\mu)\ln(\frac{(1+b)\mu-s}{(1+a)\mu-s}). \end{array} \end{equation*}

Thus, we come up with the following functional equation:

(5.2)\begin{equation} \rho(s)=\int\limits_{a}^{b}k(s,v_{1})\rho(\frac{s}{1+v_{1}})dv_{1}+L(s), \end{equation}

where

\begin{equation*} \begin{array}{rl} k(s,v):=& \frac{\phi(s)}{b-a}\frac{\mu}{\mu(1+v)-s}, \\ L(s):=& \frac{1-\phi(s)}{s}-\frac{\phi(s)}{b-a}\rho(\mu)\ln(\frac{(1+b)\mu-s}{(1+a)\mu-s}). \end{array} \end{equation*}

Iterating n times yields

(5.3)\begin{equation} \begin{array}{rl} \rho(s)=&\int\cdots\int\limits_{[a,b]^{n+1}}\prod_{j=1}^{n+1}k(\frac{s}{\prod_{i=1}^{j-1}(1+v_{i})},v_{j})\rho(\frac{s}{\prod_{m=1}^{n+1}(1+v_{m})})dv_1\cdots dv_{n+1} \\ &+L(s)+\sum_{j=1}^{n} \int\cdots\int\limits_{[a,b]^{j}}\prod_{i=1}^{j}k(\frac{s}{\prod_{m=1}^{i-1}(1+v_{m})},v_{i})L(\frac{s}{\prod_{m=1}^{j}(1+v_{m})})dv_1\cdots dv_{j}, \end{array} \end{equation}

with the convention that an empty product equals one. In the following, we will let n tend to $\infty$ to obtain an expression for $\rho(s)$. Thus, we need to estimate the limit of the first term on the RHS of (5.3), as well as to verify the convergence of the summation in the third term. The critical task to cope with this problem is the evaluation of

\begin{equation*} |\int\cdots\int\limits_{[a,b]^{n+1}}\prod_{j=1}^{n+1}k(\frac{s}{\prod_{i=1}^{j-1}(1+v_{i})},v_{j})dv_1\cdots dv_{n+1} |. \end{equation*}

Note that $k(s,v)=\frac{1}{(1+v)(b-a)}\phi(s)A(-\frac{s}{1+v})=\frac{1}{(1+v)(b-a)}E(e^{-s B+\frac{s}{1+v} A})$, where $A(s):=E(e^{-sA})$, $A\sim\exp(\mu)$, $E(e^{-sB})=\phi(s)$. Thus,

\begin{equation*} \begin{array}{l} |\int\cdots\int\limits_{[a,b]^{n+1}}\prod_{j=1}^{n+1}\frac{1}{(1+v_{j})(b-a)}\phi(\frac{s}{\prod_{i=1}^{j-1}(1+v_{i})})A(-\frac{s}{1+v_{j}})dv_1\cdots dv_{n+1}|\\ \leq |\int\cdots\int\limits_{[a,b]^{n+1}}\prod_{j=1}^{n+1}\frac{1}{1+v_{j}}dv_1\cdots dv_{n+1} |=|\int\limits_{[a,b]}\frac{1}{1+v_{1}}dv_1\cdots\int\limits_{[a,b]}\frac{1}{1+v_{n+1}}dv_{n+1}|=(\ln(\frac{b+1}{a+1}))^{n+1}. \end{array} \end{equation*}

Having in mind that $|\rho(s)|\leq 1$, the magnitude of the first term on the right-hand of (5.3) is no more than $(\ln(\frac{b+1}{a+1}))^{n+1}$, which tends to 0 as $n\to\infty$, provided that $\ln(\frac{b+1}{a+1}) \lt 1$, or equivalently $\frac{b+1}{a+1} \lt e$. Therefore,

(5.4)\begin{equation} \begin{array}{rl} \rho(s)=&L(s)+\sum_{j=1}^{\infty} \int\cdots\int\limits_{[a,b]^{j}}\prod_{i=1}^{j}k(\frac{s}{\prod_{m=1}^{i-1}(1+v_{m})},v_{i})L(\frac{s}{\prod_{m=1}^{j}(1+v_{m})})dv_1\cdots dv_{j}. \end{array} \end{equation}

Finally, $\rho(\mu)$ can be derived by substituting $s=\mu$ in (5.4), and solving the derived equation. More precisely,

\begin{equation*} \begin{array}{rl} \rho(\mu)=&\sum_{j=0}^{\infty} \int\cdots\int\limits_{[a,b]^{j}}\prod_{i=1}^{j}k(\frac{\mu}{\prod_{m=1}^{i-1}(1+v_{m})},v_{i})\\ &\times\left[\frac{1-\phi(\frac{\mu}{\prod_{m=1}^{j}(1+v_{m})})}{\frac{\mu}{\prod_{m=1}^{j}(1+v_{m})}}-\phi(\frac{\mu}{\prod_{m=1}^{j}(1+v_{m})})\frac{\rho(\mu)}{b-a}\ln\left(\frac{(1+b)\mu-\frac{\mu}{\prod_{m=1}^{j}(1+v_{m})}}{(1+a)\mu-\frac{\mu}{\prod_{m=1}^{j}(1+v_{m})}}\right)\right] dv_1\cdots dv_{j}.\end{array} \end{equation*}

Hence,

\begin{equation*} \rho(\mu)=\frac{\sum_{j=0}^{\infty} \int\cdots\int\limits_{[a,b]^{j}}\prod_{i=1}^{j}k(\frac{\mu}{\prod_{m=1}^{i-1}(1+v_{m})},v_{i})\frac{1-\phi(\frac{\mu}{\prod_{m=1}^{j}(1+v_{m})})}{\frac{\mu}{\prod_{m=1}^{j}(1+v_{m})}}dv_1\cdots dv_{j}}{1+\frac{1}{b-a}\sum_{j=0}^{\infty} \int\cdots\int\limits_{[a,b]^{j}}\prod_{i=1}^{j}k(\frac{\mu}{\prod_{m=1}^{i-1}(1+v_{m})},v_{i})\phi(\frac{\mu}{\prod_{m=1}^{j}(1+v_{m})})\ln\left(\frac{(1+b)\prod_{m=1}^{j}(1+v_{m})-1}{(1+a){\prod_{m=1}^{j}(1+v_{m})}-1}\right) dv_1\cdots dv_{j}}. \end{equation*}

6. Conclusion and future work

In this work, we cope with several nontrivial generalizations of the dual risk model with proportional gains, for which several independent assumptions among the gain interarrivals and the gain sizes were lifted, but still, we are able to derive explicit expressions for the ruin probability and the time to ruin (in terms of Laplace transforms). Among others, we considered causal dependencies as well as dependencies that are based on the FGM copula. On top of that we also consider cases where the gain size is no longer exponentially distributed, but Erlang or mixed Erlang, which result in additional interesting observations. Moreover, it is well-known that the mixed Erlang distribution belongs to a class of the phase-type distributions, which is dense in the space of distribution functions defined on $[0,\infty)$. Moreover, we cope with the case where the proportional parameter is a uniformly distributed random variable.

Note that the results we obtained so far, that is, the Laplace transform of the ruin probability, and the double Laplace transform of the time to ruin, given the initial surplus to be equal to x, seems to be valid even in the case where the gain sizes are heavy-tailed, although in such a case it is not an easy task to write the Laplace transform of the gain size in terms of elementary functions (as in the case of light-tailed gain sizes). For example, in case of a Pareto (or US-Pareto) gain size (with shape parameter greater than 1), the Laplace transform of its density is written in terms of an incomplete Gamma function; see [Reference Albrecher and Kortschak7, Reference Ramsay43, Reference Ramsay44] for related work on the classical ruin problem. See also in [Reference Bladt, Nielsen and Samorodnitsky16] for an alternative approach to cope with heavy-tailed distributions on problems related to the classical ruin theory, based on infinite-dimensional phase-type distributions. However, in [Reference Abate, Choudhury and Whitt1] a related class of random variables, called Pareto mixture of exponentials, which are defined as products of Pareto and exponential random variables, lead to explicit Laplace–Stieltjes transforms (when the shape parameter of the Pareto mixing distribution is greater than 1). This result would be useful for numerical experiments.

In the future, we plan to investigate several open tasks:

1. (1) To identify the value of the discounted cumulative dividend payments in the presence of dependence among gain interarrivals and gain sizes.

2. (2) To cope with the problem where the jumps up from the level u have a more general form.

3. (3) To consider more sophisticated copulas that describe the dependence structure, but still, to be able to derive explicit results for the values of interest. Finally, it would be interesting to consider a semi-Markovian dual risk model.

4. (4) Other opportunities for future study refer to multidimensional surplus processes, although are anticipated to be highly challenging.