On dual risk models with proportional gains and dependencies

Abstract

In this work, we consider extensions of the dual risk model with proportional gains by introducing dependence structures among gain sizes and gain interarrival times. Among others, we further consider the case where the proportionality parameter is randomly chosen, the case where it is a uniformly random variable, as well as the case where we may have upward as well as downward jumps. Moreover, we consider the case with causal dependence structure, as well as the case where the dependence is based on the generalized Farlie–Gumbel–Morgenstern copula. The ruin probability and the distribution of the time to ruin are investigated.

1. Introduction

In this work, we investigate the dual risk model with constant expense rate normalized to 1, by considering several nontrivial generalizations of the model considered in [18]. Among others, our primary aim is to lift several independence assumptions among gain interarrival times and gain sizes, but still, to be able to obtain explicit results (in terms of Laplace transforms) regarding some major metrics of interest such as ruin probability and the time to ruin.

For a detailed study on the fundamentals of ruin probabilities in the conventional, permanently inspected, Cramer–Lundberg context see [9, 36]. It is well-known that there is a duality property among the ruin theory and the queueing theory. In particular, the Cramer–Lundberg model is dual to the M/G/1 queueing model with the same arrival rate and with a service time distribution that equals the claim size distribution in the Cramer–Lundberg model; see [9, 36]. Quite recently, in the seminal book [36], the authors presented the main results and the most important probabilistic methods related to the Cramer–Lundberg model and exploited connections with the related model in queueing theory.

Its dual process has also attracted interest in the insurance risk literature. As pointed out in [12], the dual risk model describes the surplus or equity of a company with fixed expense rate and occasional income inflows of random size, called innovations or gains. These gains arise due to some contingent events (e.g., discoveries, sales). Examples where the dual risk model applies are in pharmaceutical, petroleum, or R&D companies. Other examples are commission-based businesses, such as real estate agents or brokerage firms that sell mutual funds or insurance products with a front-end load. As stated in [18, see Figure 1], an illustrative realistic example of the dual risk model with proportional gains refers to start-ups or e-companies, where their gains strongly depend on the amount of investments, which is often proportional to the value of the company; for example, the CD Projekt (see https://www.cdprojekt.com/en/), one of the biggest Polish companies producing computer games.

Figure 1. Instances of $\rho(s)$ when we truncate the infinite sum in (3.27) at k = 2 and k = 30 (when λ = 1, µ = 4, θ = 0.7, a = 3).

1.1. Related work

The standard dual risk model (without proportional gains) has been treated in the past in [32] where it was called as the negative claims model; see also [13]. Moreover, since the dual risk model corresponds to the first busy period in a single-server queue with initial workload x, we further refer [25, 42]. The ruin probability in the dual risk model under a taxation system was analyzed in [5], while in [41], the authors studied the finite-time ruin probability under a discrete-time setup. The vast majority of works related to the dual risk model focused on dividend barriers. The authors in [10, 12] considered cases where gains follow an exponential distribution or a mixture of exponential distributions, and they derived explicit formulas for the expected discounted dividend value; see also [4]. A Laplace transform method to investigate the dual model perturbed by a diffusion was applied in [11], while the authors in [14, 15] proved the optimality of a barrier strategy for all spectrally positive Lévy processes by employing fluctuation theory. A general dual risk process where the rate of the costs depends on the present amount of reserves was analyzed in [37]. In [17], the authors considered the time to ruin and the expected discounted dividends for a special dividend policy. Recently in [33], the authors studied a mixed dividend strategy, which is the combination of a threshold dividend and a Parisian implementation delays dividend under periodic observation. In [45], the dual risk model with Erlang distributed gain interarrival time was analyzed, and expressions for the ruin probability and the Laplace transform of the time of ruin for an arbitrary single gain distribution were derived. We further mention [40, 49, 50] that considered terminal costs and dividends that are paid continuously at a constant rate when the surplus is above that barrier, and [31] for the asymptotic analysis of optimal dividends in a dual risk model. Quite recently, the authors in [18] considered a variant of a dual risk model by adding the proportional gain feature, that is, if the surplus process just before a gain arrival is u, then for a > 0 the capital jumps up to the level $(1 + a)u + C$, where C is the size of the gain. The authors dealt with the ruin probability and the distribution of the time to ruin, that is, the first time the surplus level hits zero. They further identified the value of discounted cumulative dividend payments and also considered a random perturbation of the basic risk process modeled by an independent Brownian motion with drift.

1.2. Our contribution and motivation

The major contribution of this work relies on the investigation of dual risk models with a proportional gain mechanism where in addition: (a) several assumptions of independence among gain interarrivals and gain sizes are lifted, (b) the proportional parameter can be a random variable and/or taking specific values according to a probabilistic rule (and not only a fixed constant as in [18]), and still, to be able to derive explicit expressions (in terms of Laplace Transforms) for the ruin probability and the distribution of the time to ruin (given the initial surplus level). For example, among others, we investigate the case where the distribution of the gain size depends on the gain interarrival time, and this type of dependence it also affects the proportional parameter. Moreover, we focus on the case where the dependence structure is based on a copula, and the case where the proportional parameter is a random variable with a finite support. To our best knowledge, it is the first time in the related literature that the dual risk model with proportional gains and a general dependence structure is analyzed.

Let U(t) be the surplus process with $U(0) = x \gt 0$. Our interest relies on the derivation of the ruin probability, given that the initial surplus equals x, that is, $R(x): = P(\tau_{x} \lt \infty)$, where $\tau_{x} = \inf\{t\geq 0:U(t)=0|U(0)=x\}$, as well as, of the distribution of τ_x, that is, the time to ruin, which is defined as the first time the surplus equals zero. In particular, we focus on the derivation of $\tau(s,\alpha):=\int_{x=0}^{\infty}e^{-sx}E(e^{-\alpha\tau_{x}}1(\tau_{x} \lt \infty)|U(0)=x)dx$, where $1(A)$ be the indicator function of the event A. To investigate R(x) and $\tau(s,\alpha)$, we use a one-step analysis where the process under study is viewed at successive gain times. The approach we follow bears similarities to the method developed in [19, 20, 22, 29, 34] to study reflected autoregressive processes. To our best knowledge, this work provides for the first time explicit expressions for the Laplace transform of the ruin probability and the double Laplace transform of the distribution of the time to ruin, given that the $U(0)=x \gt 0$, in a general dual risk model with proportional gains under a dependent setting, and thus, should be viewed as a starting point for obtaining analytical results in more general dependent scenarios. In all models considered in this paper, we assume constant expense rate normalized to 1.

Our motivation in this work is twofold: First, by mathematical point of view, to consider several nontrivial generalizations of the model developed in [18], (a) by lifting independence assumptions among gain interarrivals and gain sizes. Among others, we introduce a dependence structure based on the Farlie–Gumbel–Morgenstern (FGM) copula that has both a mathematical and a practical interest, (b) by allowing the proportionality parameter to depend on gain interarrivals, as well as to be a constant or a random variable. Moreover, note that in queueing terms, the surplus process describes the workload in an M/G/1 queue with a constant demand rate and occasional inflow that depends proportionally on the current amount of work in the system, and where there is a dependence structure among interarrival times and the work that brought in the system. To our best knowledge, this is the first time that such a queueing model is investigated. Moreover, related works on queueing models with dependence structure based on copulas are very rare. Thus, the derived results have merit both in insurance mathematics and in queueing theory. Second, as mentioned above, by practical point of view, the concept of proportional gains, that is, gains that are proportional to the value of the company, appears lately in the budgets of start-ups and e-companies. By further considering dependence structures make our model even more realistic.

The paper is organized as follows. In Section 2, we consider the dual risk model with causal dependence structure that relates the gain size, the gain interarrival, and the proportional parameter, and obtain the Laplace transform of the ruin probability as well as the Laplace transform of the ruin time Laplace–Stieltjes transform (Laplace–Stieltjes transform). We present result where the gain size follows exponential, Erlang, or even mixed Erlang distribution. In Section 3, we focus on the dual risk model with proportional gains and where the gain size and the gain interarrival times are dependent based on the (generalized) FGM copula. We also consider the case where, in addition, there is a linear dependence among gain interarrival times and surplus level. In the latter case, we also performed some numerical examples to illustrate the ability of the derived expressions regarding the numerical calculations of ruin measures. Section 4 is devoted to the analysis of the dual risk model with randomly proportional gains with upward and downward jumps. We further consider the case of dependence based on the FGM copula among gain sizes and gain interarrivals. Finally, in Section 5, we consider the case where the proportional parameter is uniformly distributed.

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

In this section, we generalize the model in [18] 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 [6], the authors considered the classical ruin model with causal dependencies (motivated by the work in [23]), 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 [18] where the authors assumed exponentially distributed gain sizes (they just mentioned in [18, 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 [3, 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 [47]. 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 [18] 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 [39], 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., [27]), 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 [27], 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 [24], the authors generalized the work in [27], by considering an extension to the renewal or Sparre Andersen risk process, where the claim interarrivals are Erlang distributed. Later, in [26], 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 [21], where some risk models with a semi-linear dependent structure (no copula dependent) were discussed.

As stated in [39], “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., [39]) 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 [35, 39] 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, [8, 28, 38] for applications of copulas in actuarial science and financial risk management.

In this section, we generalize the model in [18], 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 [18, 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$.

Remark 3.1. Setting θ = 0, that is, by assuming independence among the gain interarrival and the gain size, and letting n = 1, that is, by considering exponentially distributed gain sizes, we recover the functional equation (2.5) in [18].

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

Theorem 3.2. The Laplace transform of the ruin probability $\rho(s)$ is given in (4.9), where the terms $\rho^{(l)}(\mu)$, $l=0,1,\ldots,n-1$, $\rho^{(l)}(2\mu)$, $l=0,1,\ldots,2n-2$, are obtained as the solution of a system of $3n-1$ equations that is derived by differentiating l times (4.9) with respect to s, and setting $s=\mu$, and $s=2\mu$, respectively.

Remark 3.3. Note that both $J(s;\theta)$ and $H(s,\theta)$ have a singularity at $s=\mu(1+a)$, and $s=2\mu(1+a)$, and thus, one expects that the expression $\rho(s)$ in (3.5) has a singularity at $s=\mu(1+a)^{j+1}$, and $s=2\mu(1+a)^{j+1}$, $j=0,1,\ldots$. However, such a singularity is removable (as in the independent case that was treated in [18, Remark 2.2]).

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 [48]. 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 [48], 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 [48] 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 [46], 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, [30]. 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 [48] 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 [48] 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 [21, 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 [18, 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 [2, 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 [21, 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 [48], 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 [21, 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 [18]. 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 [3, 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 [18]. Equation (4.9) has the same form as in [18, 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 [18, 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 [3, 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*}

Remark 5.1. Clearly, the approach we followed in this section can be used to incorporate dependencies among the gain interarrivals and the gain sizes, as well as the case where gain sizes follow Erlang, or mixed Erlang distribution.

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 [7, 43, 44] for related work on the classical ruin problem. See also in [16] 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 [1] 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.