1 Introduction and Main Results
In this paper, we consider the so called bi-seasonal discrete time risk model, which is the direct generalization of the classical discrete time risk model.
Definition 1.
We say that the insurer’s surplus ${W_{u}}$ varies according to the bi-seasonal risk model if
for each $n\in {\mathbb{N}_{0}}=\{0,1,2,\dots \}$ and the following assumptions hold:
-
• the initial insurer’s surplus $u\in {\mathbb{N}_{0}}$,
-
• the random claim amounts $\{{Z_{1}},{Z_{2}},\dots \}$ are nonnegative integer-valued independent r.v.s.,
-
• there exist r.v.s. X and Y such that ${Z_{2k+1}}\stackrel{d}{=}X$, $k\in {\mathbb{N}_{0}}$, and ${Z_{2k}}\stackrel{d}{=}Y$, $k\in \mathbb{N}$.
If $X\stackrel{d}{=}Y$, then the bi-seasonal discrete time risk model becomes the classical discrete time risk model.
There exists practical motivation for seasonal risk models in different spheres of insurance risks. In Bischoff-Ferrari et al. (2007) the effect of seasonality on fracture risk is found to be statistically significant. Another example of risk influenced by seasonality is dairy production loss risk, as found by Deng et al. (2007).
The Gerber–Shiu discounted penalty function ${\Psi _{\delta ,w}}$ is one of the main critical characteristics for risk models of any types. According to the definition presented in Gerber and Shiu (1998) for the discrete time risk model
\[ {\Psi _{\delta ,w}}(u)=\mathbb{E}\big({\mathrm{e}^{-\delta {T_{u}}}}\hspace{0.1667em}w\big({W_{u}}({T_{u}}-1),|{W_{u}}({T_{u}})|\big){\mathbb{1}_{\{{T_{u}}<\infty \}}}\big),\]
where force of interest $\delta \geqslant 0$, $w(x,y)$ is an arbitrary function of two nonnegative arguments, and ${T_{u}}$ denotes the time of ruin, i.e.
Function w has practical interpretations. For example, if w was interpreted as the benefit amount of reinsurance payable at the time of ruin, then ${\Psi _{\delta ,w}}(u)$ is the single premium of the reinsurance.
In the particular case considered in this paper when $w(x,y)=1$ for all nonnegative x and y, the discounted penalty function is equal to the following expression
If, in addition, force of interest $\delta =0$, then the Gerber–Shiu discounted penalty function is equal to the ruin probability
After Gerber and Shiu (1998) presented the concept of function named on their behalf, various properties of this function were considered by many authors. The main part of the known results on the Gerber–Shiu function is related with the Sparre Andersen model and various generalizations of this model. For instance, several cases of the Sparre Andersen model were considered by Dickson and Qazvini (2016), Landriault and Willmot (2008), Li and Garrido (2004), Li and Sendova (2015), Lin et al. (2003), Schmidli (1999), Willmot and Dickson (2003). Properties of the Gerber–Shiu function in the risk renewal models perturbed by diffusion were investigated by Chi et al. (2010), Tsai (2003), Tsai and Willmot (2002), Xu et al. (2014), Zhang and Cheung (2016), Zhang et al. (2012, 2017b, 2014). The Gerber–Shiu function of the risk models with various special strategies were considered by Avram et al. (2015), Bratiichuk (2012), Cheung and Liu (2016), Cheung et al. (2015), Dong et al. (2009), Lin and Pavlova (2006), Lin and Sendova (2008), Liu et al. (2015), Marciniak and Palmowski (2016), Shi et al. (2013), Shiraishi (2016), Woo et al. (2017), Zhang et al. (2017a), Zhou et al. (2015). This function for the risk models with various dependence structures or for risk models with investment strategies was considered by Cheung et al. (2011), Cossette et al. (2011), Li and Lu (2013), Mihálýko and Mihálýko (2011), Schmidli (2015), among others.
In the above articles, the general risk renewal models of continuous time were considered. In such a case, the defective renewal equation is the main tool to obtain a suitable information about the exact values or the asymptotic behaviour of the Gerber–Shiu function. If we consider the discrete time risk model, then the recursive relations between values of the Gerber–Shiu function play role of the defective renewal equation. Recursive methods were successfully analysed in many diverse fields, ranging from queuing models (Ferreira et al., 2017) to dynamical systems (De La Sen, 2016). Various properties of the Gerber–Shiu function in the discrete time risk models were considered by Bao and Liu (2016), Cheng et al. (2000), Li and Wu (2015), Li (2005), Li and Garrido (2002), Li et al. (2009), Liu et al. (2017), Liu and Guo (2006), Marceau (2009), Pavlova and Willmot (2004). For instance, in Li and Garrido (2002), it is shown that values of function ${\Psi _{\delta ,w}}$ of the homogeneous discrete time risk model can be calculated using the following formulas
\[\begin{aligned}{}{\Psi _{\delta ,w}}(0)=& {\mathrm{e}^{-\delta }}{\sum \limits_{k=0}^{\infty }}{\sum \limits_{l=0}^{\infty }}{\varrho ^{k}}w(k,l)\mathbb{P}(Z=k+l+1),\\ {} {\Psi _{\delta ,w}}(u)=& {\mathrm{e}^{-\delta }}{\sum \limits_{k=0}^{u-1}}{\Psi _{\delta ,w}}(u-k){\sum \limits_{l=0}^{\infty }}{\varrho ^{l}}\hspace{0.1667em}\mathbb{P}(Z=k+l+1)\\ {} & +{\mathrm{e}^{-\delta }}{\varrho ^{-u}}{\sum \limits_{k=u}^{\infty }}{\varrho ^{k}}{\sum \limits_{l=0}^{\infty }}w(k,l)\mathbb{P}(Z=k+l+1),\end{aligned}\]
where Z with $\mathbb{E}Z<1$ is the integer-valued random variable generating the homogeneous discrete time risk model, and $\varrho \in (0,1)$ is the root of equation
By arguments provided in Li and Garrido (2002), such a solution exists and is unique for $\delta >0$.If the discrete time risk model is generated by possibly differently distributed random variables ${Z_{1}},{Z_{2}},\dots $ , then the above formulas do not hold anymore. The situation in the nonhomogeneous discrete time risk model is much more complicated.
In this paper, we consider the behaviour of the special case of Gerber–Shiu penalty function for the bi-seasonal discrete time risk model which is a particular case of nonhomogeneous discrete time risk models. Our results supplement the results of Castañer et al. (2013), Răducan et al. (2015a) and Răducan et al. (2015b). We derive the specific recursive equality for function ${\psi _{\delta }}$. Using the derived formula we construct an algorithm to calculate approximate values of this function. The running of the algorithm is illustrated by several examples. The ideas from Bieliauskienė and Šiaulys (2012), Damarackas and Šiaulys (2014), De Vylder and Goovaerts (1988), Dickson and Waters (1991) were used to get the main results of this paper.
We consider the bi-seasonal discrete time risk model generated by two nonnegative, independent and integer valued random variables X and Y. By
\[ {x_{k}}=\mathbb{P}(X=k),\hspace{2em}{y_{k}}=\mathbb{P}(Y=k),\hspace{2em}{q_{k}}=\mathbb{P}(Q=k),\hspace{1em}k\in {\mathbb{N}_{0}}\]
we denote the local probabilities of random variables X, Y and $Q=X+Y$ respectively. Distribution functions of these random variables we denote by ${F_{X}}$, ${F_{Y}}$ and ${F_{Q}}$, i.e.
\[\begin{array}{l}\displaystyle {F_{X}}(u)=\mathbb{P}(X\leqslant u)={\sum \limits_{k=0}^{\lfloor u\rfloor }}{x_{k}},\\ {} \displaystyle {F_{Y}}(u)=\mathbb{P}(Y\leqslant u)={\sum \limits_{k=0}^{\lfloor u\rfloor }}{y_{k}},\\ {} \displaystyle {F_{Q}}(u)=\mathbb{P}(Q\leqslant u)={\sum \limits_{k=0}^{\lfloor u\rfloor }}{q_{k}},\end{array}\]
for each real u. The notation $\overline{F}$ is used for the tail of an arbitrary distribution function F, i.e. $\overline{F}(u)=1-F(u)$ for each $u\in \mathbb{R}$.The following two assertions enable us to construct an algorithm for calculating values of function ${\psi _{\delta }}(u)$ in the bi-seasonal discrete time risk model.
Theorem 1.
Let the bi-seasonal discrete time risk model be generated by two nonnegative, independent and integer valued random variables X and Y. If $\mathbb{E}X+\mathbb{E}Y<2$, then ${\lim \nolimits_{u\to \infty }}{\psi _{\delta }}(u)=0$ for an arbitrary fixed $\delta \geqslant 0$. In addition, if $\max \{\mathbb{E}{\mathrm{e}^{hX}},\mathbb{E}{\mathrm{e}^{hY}}\}<\infty $ for some positive h, then ${\textstyle\sum _{l=0}^{\infty }}{\psi _{\delta }}(l)<\infty $ for each fixed $\delta \geqslant 0$.
Theorem 2.
Let all the conditions of Theorem 1 be satisfied. Furthermore, let $\delta >0$, and ${\psi _{\delta }}$ denote the Gerber–Shiu function with $w(x,y)=1$ for all nonnegative x and y. Also denote ${\mathcal{S}_{\delta }}:={\textstyle\sum _{l=0}^{\infty }}{\psi _{\delta }}(l)$.
for each $n\in {\mathbb{N}_{0}}$, where ${a_{n}}$, ${b_{n}}$, ${d_{n}}$ are three sequences of real numbers defined recursively by the following equations:
\[\begin{array}{l}\displaystyle {a_{0}}=1,\hspace{2em}{a_{1}}=-\frac{1}{{y_{0}}},\hspace{2em}{a_{n}}=\frac{1}{{q_{0}}}\bigg({\mathrm{e}^{2\delta }}{a_{n-2}}-{\sum \limits_{i=1}^{n-1}}{q_{i}}{a_{n-i}}-{x_{n-1}}\bigg),\\ {} \displaystyle \hspace{1em}n\in \{2,3,\dots \};\\ {} \displaystyle {b_{0}}=0,\hspace{2em}{b_{1}}=-\frac{{\mathrm{e}^{2\delta }}-1}{{y_{0}}},\\ {} \displaystyle {b_{n}}=\frac{1}{{q_{0}}}\bigg({\mathrm{e}^{2\delta }}{b_{n-2}}-{\sum \limits_{i=1}^{n-1}}{q_{i}}{b_{n-i}}-{x_{n-1}}({\mathrm{e}^{2\delta }}-1)\bigg),\hspace{1em}n\in \{2,3,\dots \};\\ {} \displaystyle {d_{0}}=0,\hspace{1em}{d_{1}}=\frac{{\mathrm{e}^{\delta }}\mathbb{E}X+{y_{0}}+\mathbb{E}Y-1}{{y_{0}}},\\ {} \displaystyle {d_{n}}=\frac{1}{{q_{0}}}\bigg({\mathrm{e}^{2\delta }}{d_{n-2}}-{\sum \limits_{i=1}^{n-1}}{q_{i}}{d_{n-i}}+{x_{n-1}}{y_{0}}{d_{1}}-{\mathrm{e}^{\delta }}{\overline{F}_{X}}(n-2)-{\sum \limits_{i=0}^{n-2}}{x_{i}}{\overline{F}_{Y}}(n-1-i)\bigg),\\ {} \displaystyle \hspace{1em}n\in \{2,3,\dots \}.\end{array}\]
for each $n\in {\mathbb{N}_{0}}$, where ${\tilde{a}_{n}}$, ${\tilde{b}_{n}}$, ${\tilde{d}_{n}}$ are three sequences of real numbers defined recursively by the following equations:
\[\begin{array}{l}\displaystyle {\tilde{a}_{0}}=1,\hspace{2em}{\tilde{a}_{1}}=-\frac{1}{{y_{0}}},\hspace{2em}{\tilde{a}_{n}}=\frac{1}{{q_{1}}}\bigg({\mathrm{e}^{2\delta }}{\tilde{a}_{n-1}}-{\sum \limits_{i=1}^{n-1}}{q_{i+1}}{\tilde{a}_{n-i}}-{x_{n}}\bigg),\\ {} \displaystyle n\in \{2,3,\dots \};\\ {} \displaystyle {\tilde{b}_{0}}=0,\hspace{2em}{\tilde{b}_{1}}=-\frac{{\mathrm{e}^{2\delta }}-1}{{y_{0}}},\\ {} \displaystyle {\tilde{b}_{n}}=\frac{1}{{q_{1}}}\bigg({\mathrm{e}^{2\delta }}{\tilde{b}_{n-1}}-{\sum \limits_{i=1}^{n-1}}{q_{i+1}}{\tilde{b}_{n-i}}-{x_{n}}({\mathrm{e}^{2\delta }}-1)\bigg),\hspace{1em}n\in \{2,3,\dots \};\\ {} \displaystyle {\tilde{d}_{0}}=0,\hspace{2em}{\tilde{d}_{1}}=\frac{{\mathrm{e}^{\delta }}\mathbb{E}X+{y_{0}}+\mathbb{E}Y-1}{{y_{0}}},\\ {} \displaystyle {\tilde{d}_{n}}=\frac{1}{{q_{1}}}\bigg({\mathrm{e}^{2\delta }}{\tilde{d}_{n-1}}-{\sum \limits_{i=1}^{n-1}}{q_{i+1}}{\tilde{d}_{n-i}}+{x_{n}}{y_{0}}{\tilde{d}_{1}}-{\mathrm{e}^{\delta }}{\overline{F}_{X}}(n-1)\\ {} \displaystyle \phantom{{\tilde{d}_{n}}=}-{\sum \limits_{i=0}^{n-2}}{x_{i+1}}{\overline{F}_{Y}}(n-1-i)\bigg),\hspace{1em}n\in \{2,3,\dots \}.\end{array}\]
for each $n\in {\mathbb{N}_{0}}$, where ${\hat{b}_{n}}$, ${\hat{d}_{n}}$ are two sequences of real numbers defined recursively by the following equations:
\[\begin{array}{l}\displaystyle {\hat{b}_{0}}=-({\mathrm{e}^{2\delta }}-1),\\ {} \displaystyle {\hat{b}_{n}}=\frac{1}{{q_{1}}}\bigg({\mathrm{e}^{2\delta }}{\hat{b}_{n-1}}-{\sum \limits_{i=1}^{n-1}}{q_{i+1}}{\hat{b}_{n-i}}\bigg),\hspace{1em}n\in \mathbb{N};\\ {} \displaystyle {\hat{d}_{0}}={\mathrm{e}^{\delta }}\mathbb{E}X+\mathbb{E}Y-1,\\ {} \displaystyle {\hat{d}_{n}}=\frac{1}{{q_{1}}}\bigg({\mathrm{e}^{2\delta }}{\hat{d}_{n-1}}-{\sum \limits_{i=1}^{n-1}}{q_{i+1}}{\hat{d}_{n-i}}-{\mathrm{e}^{\delta }}{\overline{F}_{X}}(n-1)-{\sum \limits_{i=0}^{n-1}}{x_{i}}{\overline{F}_{Y}}(n-i)\bigg),\hspace{1em}n\in \mathbb{N}.\end{array}\]
Remark 1.
We observe that case ${x_{0}}={y_{0}}=0$ is impossible due to requirement $\mathbb{E}(X+Y)<2$. This observation shows that all possible cases of the discrete r.v.s. X and Y are considered in Theorem 2.
The rest of the paper is organized in the following way: in Section 2 we describe an algorithm for calculating values of Gerber–Shiu function; next, in Section 3 we present a few numerical examples which illustrate the applicability of our algorithm; in Section 4 some concluding remarks and directions for future work are provided; Section 5 deals with proofs of the main results; in Section 6 some lower and upper bounds for Gerber–Shiu function are derived; finally, in Section 7 the algorithm code in R language is provided.
2 Algorithm for Finding the Values of Function ${\psi _{\delta }}$
In this section, we describe an algorithm for calculating values of ${\psi _{\delta }}(u)$ in the case of the bi-seasonal risk model. The algorithm was implemented with R language, using increased numerical precision package Rmpfr. Our algorithm is based on formula (1) from Theorem 2 and the results of Theorem 1. As usual, it is assumed that we have a positive force of interest δ, and the bi-seasonal discrete time risk model is generated by two nonnegative, integer-valued and differently distributed r.v.s. X, Y with local probabilities ${x_{k}}=\mathbb{P}(X=k)$, ${y_{k}}=\mathbb{P}(Y=k)$, $k\in {\mathbb{N}_{0}}$. Of course, these two r.v.s. should satisfy all requirements of Theorem 2. Below we present the detailed, step by step algorithm for calculating ${\psi _{\delta }}(u)$, $u\in {\mathbb{N}_{0}}$ in the case when ${x_{0}}{y_{0}}>0$. The other possible cases:$\{{x_{0}}=0,{y_{0}}>0\}$, and $\{{x_{0}}>0,{y_{0}}=0\}$, which were described in Theorem 2, can be considered similarly.
Step 1: Select $N\in \{10,20,30,\dots ,100\}$ and $K\in \{1,\dots ,5\}$.
Step 2: Calculate coefficients ${a_{n}}$, ${b_{n}}$, ${d_{n}}$ for all $n\in \{0,1,\dots ,N\}$ using formulas from Theorem 2.
Step 3: Find ${\hat{\psi }_{\delta }}(0)$ and ${\hat{\mathcal{S}}_{\delta }}$ satisfying the following system of linear equations
Due to the main formula (1) of Theorem 2 the desired quantity ${\psi _{\delta }}(0)$ together with sum ${\mathcal{S}_{\delta }}$ satisfy the following system
However, according to Theorem 1 ${\psi _{\delta }}(N-K)$ and ${\psi _{\delta }}(N)$ are close to zero for sufficiently large N. We get system (4) from (5) by changing values of ${\psi _{\delta }}(N-K)$ and ${\psi _{\delta }}(N)$ to zeroes.
(5)
\[ \left\{\begin{array}{l}{a_{N-K}}{\psi _{\delta }}(0)+{b_{N-K}}{\mathcal{S}_{\delta }}+{d_{N-K}}={\psi _{\delta }}(N-K),\\ {} {a_{N}}{\psi _{\delta }}(0)+{b_{N}}{\mathcal{S}_{\delta }}+{d_{N}}={\psi _{\delta }}(N).\end{array}\right.\]
Step 4: Test the error $|{\psi _{\delta }}(0)-{\hat{\psi }_{\delta }}(0)|$.
Using the Cramer’s rule for both systems of linear equations (4), (5) and the trivial estimate $|{\psi _{\delta }}(n)|\leqslant 1$, $n\in {\mathbb{N}_{0}}$, we derive that
\[ \big|{\psi _{\delta }}(0)-{\hat{\psi }_{\delta }}(0)\big|\leqslant \frac{{e^{-\delta }}\big(|{b_{N-K}}|+|{b_{N}}|\big)}{|{a_{N-K}}{b_{N}}-{b_{N-K}}{a_{N}}|}.\]
Numerical simulations have shown that the upper estimate of ${\psi _{\delta }}(0)$ approximation error tends to 0 as N grows. This is consistent with the behaviour of the approximation error itself. As for parameter K, its choice does not have a clear effect on the upper estimate of ${\psi _{\delta }}(0)$ approximation error.
Step 5: If the size of error in Step 4 is suitable, then pass to Step 6. If the size of error is not suitable, then return to Step 1 choosing different parameters N and K.
We remark only that the sets provided in Step 1 for choosing these parameters are not strictly defined, and different sets can be used successfully. However, choosing N much larger than 100 would result in very large coefficients ${a_{N}}$, ${b_{N}}$ and ${d_{N}}$, and owing to that some computational difficulties may arise. Besides that, in this case computational speed would be reduced. And conversely, choosing N too small would result in big approximation error of ${\psi _{\delta }}(0)$ when changing system (5) to (4), since ${\psi _{\delta }}(N)$ does not converge to zero so quickly. As for parameter K, it should be chosen to minimize the upper estimate of ${\psi _{\delta }}(0)$ approximation error.
Step 6: Calculate ${\psi _{\delta }}(1)$ according to the formula (1) by supposing that ${\psi _{\delta }}(0)={\hat{\psi }_{\delta }}(0)$ and ${\mathcal{S}_{\delta }}={\hat{\mathcal{S}}_{\delta }}$.
Step 7: Calculate values of ${\psi _{\delta }}(u)$ for $u\geqslant 2$ while the algorithm works correctly, applying either formula (1) from Theorem 2 or the main recursive formula (7) from the proof of Theorem 2.
By saying that the algorithm works correctly, we mean that its results do not conflict with mathematical properties. Namely, ${\psi _{\delta }}(u)$ is a function taking values between 0 and 1, nonincreasing with respect to u and decreasing with respect to δ. However, sometimes algorithm produces results that are not compatible with these properties. This could happen due to the following reasons:
-
• In some particular cases of X, Y and δ, coefficients ${a_{n}}$, ${b_{n}}$, ${d_{n}}$, $n\in {\mathbb{N}_{0}}$ in the main equality of Theorem 2 are rapidly growing and fluctuating. Consequently, it is quite difficult to get precise values of these coefficients.
Remark 2.
Many ideas for constructing a recursive algorithm were taken from Damarackas and Šiaulys (2014). In this article infinite time ruin probability, which is a special case of Gerber–Shiu function with $\delta =0$ and $w(x,y)=1$, was considered. In the paper we have extended the results to the case $\delta >0$.
Remark 3.
In Bieliauskienė and Šiaulys (2012), an analogous problem to ours is considered. While we analyse a less general model than the one provided in their paper, there are some advantages in our algorithm. Namely, our approach of finding ${\psi _{\delta }}(0)$ is more efficient. The formula provided in Theorem 3 of Bieliauskienė and Šiaulys (2012) is applicable to all numerical examples of Section 3 except the last one, which deals with random variables having infinite support. But the problem with this formula is its combinatorial form, and even for relatively simple distributions it is not easy to implement. The computational speed is also reduced for the same reason. Furthermore, our proposed algorithm is less prone to computational errors, because we do not use multiple way recursion.
3 Numerical Examples
In this section, we present four numerical examples for calculating the values of ${\psi _{\delta }}(u)$, $u\in {\mathbb{N}_{0}}$, in the bi-seasonal discrete time risk model. In all examples we consider function ${\psi _{\delta }}$ with three different values of the interest force $\delta \in \{0;0.01;0.1\}$. Our algorithm does not allow to compute function values for case $\delta =0$, so the algorithm and its ${\psi _{\delta }}(0)$ approximation error upper estimate provided in Damarackas and Šiaulys (2014) were used for this case. Since the function ${\psi _{\delta }}(u)$ seems to decay exponentially, all the figures are plotted in log scale (with base 10).
Example 1.
Let us assume that the bi-seasonal discrete time risk model is generated by the following independent random claim amounts X and Y 
In this example, both claim amounts are “good” because $\max \{\mathbb{E}X,\mathbb{E}Y\}<1$, and all conditions of Theorem 2 are satisfied. Using the algorithm presented in Section 2 we obtain values of ${\psi _{\delta }}(u)$ for $u\in \{0,1,\dots ,15\}$. These values are presented in Table 1 and are shown in (Fig. 1). The upper estimate of ${\psi _{\delta }}(0)$ approximation error, described in Step 4 of algorithm, is provided in the parenthesis near the value of δ. The results of this example are based on the value of ${\psi _{\delta }}(0)$ which is obtained with $N=50$ and $K=2$.
Table 1
Values of ${\psi _{\delta }}(u)$ in Example 1.
| u | $\delta =0$ (<0.000000001) | $\delta =0.01$ (0.083494161) | $\delta =0.1$ (0.000001786) |
| 0 | 0.735808540 | 0.715289725 | 0.588111815 |
| 1 | 0.528382921 | 0.505099453 | 0.379732449 |
| 2 | 0.308008652 | 0.283691781 | 0.168950439 |
| 3 | 0.186932507 | 0.166883336 | 0.082819297 |
| 4 | 0.109425467 | 0.094115383 | 0.036822099 |
| 5 | 0.064774209 | 0.053789118 | 0.016949434 |
| 6 | 0.038352631 | 0.030752904 | 0.007818717 |
| 7 | 0.022665488 | 0.017539770 | 0.003572849 |
| 8 | 0.013406572 | 0.010015276 | 0.001640920 |
| 9 | 0.007928948 | 0.005717783 | 0.000753055 |
| 10 | 0.004688946 | 0.003263965 | 0.000345342 |
| 11 | 0.002773172 | 0.001863371 | 0.000158466 |
| 12 | 0.001639884 | 0.001063758 | 0.000072701 |
| 13 | 0.000970174 | 0.000607275 | 0.000033353 |
| 14 | 0.000573054 | 0.000346681 | 0.000015302 |
| 15 | 0.000340345 | 0.000197913 | 0.000007020 |
Example 2.
Suppose now that the bi-seasonal discrete time risk model is generated by r.v.s. X and Y having the following distributions 
We observe that $\mathbb{E}X<1$, $\mathbb{E}Y\geqslant 1$, but $\mathbb{E}X+\mathbb{E}Y<2$ in this case. Consequently, the model is “good” only on average and all conditions of Theorem 2 are satisfied. Using the algorithm from Section 2, Table 2 is filled out with values of ${\psi _{\delta }}(u)$ for $u\in \{0,1,\dots ,15\}$. Results of this example are based on the value of ${\psi _{\delta }}(0)$ which is obtained with $N=40$ and $K=1$. Values of ${\psi _{\delta }}(u)$ are also shown in (Fig. 2).
Example 3.
Let us consider the mirror reflection of the bi-seasonal discrete time risk model from Example 2, i.e. the order of claims appearance is reversed.
From the obtained calculations we can easily see that when the positions of claims are changed, the values of ${\psi _{\delta }}(u)$ are also changing. The numerical values of ${\psi _{\delta }}(u)$ of this model are given in the Table 2 and shown in (Fig. 2) with $N=50$ and $K=3$.
| u | Example 2 | Example 2 | Example 2 | Example 3 | Example 3 | Example 3 |
| $\delta =0$ | $\delta =0.01$ | $\delta =0.1$ | $\delta =0$ | $\delta =0.01$ | $\delta =0.1$ | |
| (<0.000000001) | $(0.006459348)$ | (<0.000000001) | (<0.000000001) | $(0.001104494)$ | (<0.000000001) | |
| 0 | 0.850000000 | 0.826902130 | 0.697524567 | 0.950000000 | 0.936126346 | 0.839178292 |
| 1 | 0.500000000 | 0.455345718 | 0.274354439 | 0.625000000 | 0.588031587 | 0.427209666 |
| 2 | 0.250000000 | 0.207339723 | 0.075270358 | 0.312500000 | 0.267757665 | 0.117206868 |
| 3 | 0.125000000 | 0.094411255 | 0.020650757 | 0.156250000 | 0.121922306 | 0.032156225 |
| 4 | 0.062500010 | 0.042989761 | 0.005665627 | 0.078125010 | 0.055516800 | 0.008822203 |
| 5 | 0.031250000 | 0.019575203 | 0.001554390 | 0.039062510 | 0.025279337 | 0.002420411 |
| 6 | 0.015625000 | 0.008913485 | 0.000426454 | 0.019531260 | 0.011510838 | 0.000664050 |
| 7 | 0.007812502 | 0.004058717 | 0.000116999 | 0.009765629 | 0.005241411 | 0.000182185 |
| 8 | 0.003906251 | 0.001848120 | 0.000032099 | 0.004882816 | 0.002386654 | 0.000049983 |
| 9 | 0.001953125 | 0.000841533 | 0.000008807 | 0.002441409 | 0.001086753 | 0.000013713 |
| 10 | 0.000976563 | 0.000383189 | 0.000002416 | 0.001220706 | 0.000494848 | 0.000003762 |
| 11 | 0.000488281 | 0.000174483 | 0.000000663 | 0.000610354 | 0.000225327 | 0.000001032 |
| 12 | 0.000244141 | 0.000079450 | 0.000000182 | 0.000305178 | 0.000102602 | 0.000000283 |
| 13 | 0.000122070 | 0.000036177 | 0.000000050 | 0.000152590 | 0.000046719 | 0.000000078 |
| 14 | 0.000061035 | 0.000016473 | 0.000000014 | 0.000076296 | 0.000021273 | 0.000000021 |
| 15 | 0.000030518 | 0.000007501 | 0.000000004 | 0.000038149 | 0.000009687 | 0.000000006 |
Example 4.
Suppose that the bi-seasonal discrete time risk model is generated by r.v.s. X and Y, where X has Poisson distribution with parameter $\lambda =0.8$ and Y has geometric distribution with parameter $p=0.7$.
In this case, the model generators have infinite supports, but all requirements of Theorem 2 are satisfied. So we can use the algorithm from Section 2 to calculate values of ${\psi _{\delta }}(u)$. These values are given in Table 3 and shown in (Fig. 3). The results are obtained by choosing $N=60$ and $K=4$ in the first step of the algorithm.
Table 3
Values of ${\psi _{\delta }}(u)$ in Example 4.
| u | $\delta =0$ (<0.000000001) | $\delta =0.01$ (0.089541014) | $\delta =0.1$ (0.000002568) |
| 0 | 0.678504300 | 0.667146224 | 0.582922968 |
| 1 | 0.357239100 | 0.346815995 | 0.278446415 |
| 2 | 0.170682700 | 0.162951735 | 0.116632815 |
| 3 | 0.080801850 | 0.075772347 | 0.047817117 |
| 4 | 0.038827470 | 0.035788750 | 0.020007214 |
| 5 | 0.018862780 | 0.017104346 | 0.008536891 |
| 6 | 0.009203741 | 0.008213946 | 0.003676915 |
| 7 | 0.004496317 | 0.003949953 | 0.001588588 |
| 8 | 0.002197207 | 0.001900018 | 0.000686862 |
| 9 | 0.001073798 | 0.000913991 | 0.000297021 |
| 10 | 0.000524834 | 0.000439670 | 0.000128443 |
| 11 | 0.000256585 | 0.000211501 | 0.000055544 |
| 12 | 0.000125498 | 0.000101741 | 0.000024019 |
| 13 | 0.000061448 | 0.000048942 | 0.000010387 |
| 14 | 0.000030139 | 0.000023543 | 0.000004492 |
| 15 | 0.000014871 | 0.000011325 | 0.000001942 |
4 Concluding Remarks
In this work, the bi-seasonal discrete time risk model is considered. We derived a recursive algorithm for calculating the values of a special case of Gerber–Shiu discounted penalty function. Theoretical results are illustrated by some numerical examples.
The results obtained in this paper could be improved in the following directions:
-
• Instead of taking $w(x,y)=1$ in the Gerber–Shiu function, arbitrary function $w(x,y)$ could be taken. This would allow to reflect insurer’s economic costs at the time of ruin in a more realistic way.
-
• Our results could be generalized to the models with more complex structure of claims’ non-homogeneity. For instance, models with cyclically distributed claims with an arbitrary cycle length could be considered. In this paper, the model with cycle length equal to 2 is considered.
-
• Other model extensions may be useful to consider. For example, a model may include investment strategies on premiums or claims following some dependence structure.
-
• In the Step 4 of our presented algorithm, more subtle estimation of ${\psi _{\delta }}(0)$ approximation error could be derived.
-
• In the bi-seasonal discrete time risk model, claims with distributions satisfying $\mathbb{E}X+\mathbb{E}Y\geqslant 2$ could be considered. The difficulty arises here because limiting relations in Theorem 1 and Theorem 2 do not hold anymore. Therefore an alternate way of finding ${\psi _{\delta }}(0)$ and ${\psi _{\delta }}(1)$ should be derived.
5 Proofs of the Main Results
Proof of Theorem 1.
According to Theorem 2.3 of (Damarackas and Šiaulys, 2014), we have that
This implies that ${\lim \nolimits_{u\to \infty }}{\psi _{\delta }}(u)=0$ for an arbitrary fixed $\delta \geqslant 0$, because $0\leqslant {\psi _{\delta }}(u)\leqslant \psi (u)$ for all $\delta ,u\geqslant 0$.
Furthermore, for $i\in \mathbb{N}$ denote ${\eta _{i}}={Z_{i}}-1$. The conditions of Theorem 1 imply that
\[\begin{array}{l}\displaystyle \underset{i\in \mathbb{N}}{\sup }\mathbb{E}\big({\mathrm{e}^{h{\eta _{i}}}}\big)=\max \big\{\mathbb{E}\big({\mathrm{e}^{hX}}\big),\mathbb{E}\big({\mathrm{e}^{hY}}\big)\big\}<\infty ,\\ {} \displaystyle \underset{u\to \infty }{\lim }\underset{i\in \mathbb{N}}{\sup }\mathbb{E}\big(|{\eta _{i}}|{\mathbb{1}_{\{{\eta _{i}}\leqslant -u\}}}\big)\\ {} \displaystyle \hspace{1em}=\underset{u\to \infty }{\lim }\max \big\{\mathbb{E}\big((1-X){\mathbb{1}_{\{X\leqslant 1-u\}}}\big),\mathbb{E}\big((1-Y){\mathbb{1}_{\{Y\leqslant 1-u\}}}\big)\big\}=0,\\ {} \displaystyle \underset{n\to \infty }{\limsup }\frac{1}{n}{\sum \limits_{i=1}^{n}}\mathbb{E}{\eta _{i}}=\frac{\mathbb{E}X+\mathbb{E}Y-2}{2}<0.\end{array}\]
Hence, according to Lemma 1 by (Andrulytė et al., 2015), we have
\[ {\psi _{\delta }}(u)\leqslant \psi (u)=\mathbb{P}\bigg(\underset{k\geqslant 1}{\sup }{\sum \limits_{i=1}^{k}}{\eta _{i}}>u\bigg)\leqslant {c_{1}}{\mathrm{e}^{-{c_{2}}u}},\hspace{1em}u\geqslant 0,\]
for some positive constants ${c_{1}},{c_{2}}$.Therefore, it follows immediately that ${\textstyle\sum _{l=0}^{\infty }}{\psi _{\delta }}(l)<\infty $ for each fixed $\delta \geqslant 0$. □
Proof of Theorem 2.
Suppose that $\delta >0$ and $u\in {\mathbb{N}_{0}}$. According to definition of function ${\psi _{\delta }}$ we have
\[\begin{aligned}{}{\psi _{\delta }}(u)=& {\sum \limits_{m=1}^{\infty }}\mathbb{E}\big({\mathrm{e}^{-\delta m}}{\mathbb{1}_{\{{T_{u}}=m\}}}\big)\\ {} =& {\sum \limits_{m=1}^{\infty }}{\mathrm{e}^{-\delta m}}\hspace{0.1667em}\mathbb{P}\bigg({\sum \limits_{i=1}^{j}}{Z_{i}}<j+u\hspace{2.5pt}\mathrm{for}\hspace{2.5pt}j\in \{1,2,\dots ,m-1\}\\ {} & \mathrm{and}\hspace{2.5pt}{\sum \limits_{i=1}^{m}}{Z_{i}}\geqslant m+u\bigg)\\ {} =& {\mathrm{e}^{-\delta }}\hspace{0.1667em}\mathbb{P}({Z_{1}}\geqslant 1+u)+{\mathrm{e}^{-2\delta }}\hspace{0.1667em}\mathbb{P}({Z_{1}}<1+u,{Z_{1}}+{Z_{2}}\geqslant 2+u)\\ {} & +{\sum \limits_{m=3}^{\infty }}{\mathrm{e}^{-\delta m}}\hspace{0.1667em}\mathbb{P}\bigg({\sum \limits_{i=1}^{j}}{Z_{i}}<j+u\hspace{2.5pt}\mathrm{for}\hspace{2.5pt}j\in \{1,2,\dots ,m-1\}\\ {} & \mathrm{and}\hspace{2.5pt}{\sum \limits_{i=1}^{m}}{Z_{i}}\geqslant m+u\bigg).\end{aligned}\]
Since $X\stackrel{d}{=}{Z_{1}}\stackrel{d}{=}{Z_{3}}\stackrel{d}{=}{Z_{5}}\stackrel{d}{=}\dots \hspace{0.1667em}$ and $Y\stackrel{d}{=}{Z_{2}}\stackrel{d}{=}{Z_{4}}\stackrel{d}{=}{Z_{6}}\stackrel{d}{=}\dots \hspace{0.1667em}$ we get that
(6)
\[\begin{aligned}{}{\psi _{\delta }}(u)=& {\mathrm{e}^{-\delta }}\sum \limits_{l\geqslant 1+u}{x_{l}}+{\mathrm{e}^{-2\delta }}\sum \limits_{l\leqslant u}\sum \limits_{k\geqslant 2+u-l}{x_{l}}{y_{k}}\\ {} & +{\sum \limits_{m=3}^{\infty }}{\mathrm{e}^{-\delta m}}\hspace{0.1667em}\mathbb{P}\bigg({Z_{1}}\leqslant u,{Z_{1}}+{Z_{2}}\leqslant 1+u,{Z_{1}}+{Z_{2}}+{\sum \limits_{i=3}^{j}}{Z_{i}}<j+u\\ {} & \mathrm{for}\hspace{2.5pt}j\in \{3,\dots ,m-1\}\hspace{2.5pt}\mathrm{and}\hspace{2.5pt}{Z_{1}}+{Z_{2}}+{\sum \limits_{i=3}^{m}}{Z_{i}}\geqslant m+u\bigg)\\ {} =& {\mathrm{e}^{-\delta }}{\overline{F}_{X}}(u)+{\mathrm{e}^{-2\delta }}{\sum \limits_{l=0}^{u}}{x_{l}}{\overline{F}_{Y}}(1+u-l)\\ {} & +{\sum \limits_{l=0}^{u}}{\sum \limits_{k=0}^{1+u-l}}{x_{l}}{y_{k}}{\sum \limits_{m=3}^{\infty }}{\mathrm{e}^{-\delta m}}\hspace{0.1667em}\mathbb{P}\bigg({\sum \limits_{i=3}^{j}}{Z_{i}}<j+u-l-k\\ {} & \mathrm{for}\hspace{2.5pt}j\in \{3,\dots ,m-1\}\hspace{2.5pt}\mathrm{and}\hspace{2.5pt}{\sum \limits_{i=3}^{m}}{Z_{i}}\geqslant m+u-k-l\bigg)\\ {} =& {\mathrm{e}^{-\delta }}{\overline{F}_{X}}(u)+{\mathrm{e}^{-2\delta }}{\sum \limits_{l=0}^{u}}{x_{l}}{\overline{F}_{Y}}(1+u-l)\\ {} & +{\mathrm{e}^{-2\delta }}{\sum \limits_{l=0}^{u}}{\sum \limits_{k=0}^{1+u-l}}{x_{l}}{y_{k}}{\sum \limits_{m=3}^{\infty }}{\mathrm{e}^{-\delta (m-2)}}\hspace{0.1667em}\mathbb{P}\bigg({\sum \limits_{i=1}^{j}}{Z_{i}}<j+u-l-k\\ {} & \mathrm{for}\hspace{2.5pt}j\in \{1,\dots ,m-3\}\hspace{2.5pt}\mathrm{and}\hspace{2.5pt}{\sum \limits_{i=1}^{m-2}}{Z_{i}}\geqslant m+u-k-l\bigg)\\ {} =& {\mathrm{e}^{-\delta }}{\overline{F}_{X}}(u)+{\mathrm{e}^{-2\delta }}{\sum \limits_{l=0}^{u}}{x_{l}}{\overline{F}_{Y}}(1+u-l)\\ {} & +{\mathrm{e}^{-2\delta }}{\sum \limits_{l=0}^{u}}{\sum \limits_{k=0}^{1+u-l}}{x_{l}}{y_{k}}\hspace{0.1667em}{\psi _{\delta }}(u+2-k-l).\end{aligned}\]For each $m\in {\mathbb{N}_{0}}$
Therefore the last sum in equality (6) can be expressed by
\[\begin{array}{l}\displaystyle {\sum \limits_{l=0}^{1+u}}{\sum \limits_{k=0}^{1+u-l}}{x_{l}}{y_{k}}\hspace{0.1667em}{\psi _{\delta }}(u+2-(k+l))-{x_{u+1}}{y_{0}}{\psi _{\delta }}(1)\\ {} \displaystyle \hspace{1em}={\sum \limits_{l=0}^{1+u}}{q_{l}}{\psi _{\delta }}(u+2-l)-{x_{u+1}}{y_{0}}{\psi _{\delta }}(1)\\ {} \displaystyle \hspace{1em}={\sum \limits_{l=0}^{u}}{q_{l}}{\psi _{\delta }}(u+2-l)+({q_{u+1}}-{x_{u+1}}{y_{0}}){\psi _{\delta }}(1).\end{array}\]
Substituting this expression into equality (6) we obtain that
for each $u\in {\mathbb{N}_{0}}$.
(7)
\[\begin{aligned}{}{\psi _{\delta }}(u)=& {\mathrm{e}^{-\delta }}{\overline{F}_{X}}(u)+{\mathrm{e}^{-2\delta }}{\sum \limits_{l=0}^{u}}{x_{l}}{\overline{F}_{Y}}(1+u-l)\\ {} & +{\mathrm{e}^{-2\delta }}\bigg({\sum \limits_{l=0}^{u}}{q_{l}}{\psi _{\delta }}(u+2-l)+({q_{u+1}}-{x_{u+1}}{y_{0}}){\psi _{\delta }}(1)\bigg)\end{aligned}\]By summing these last equalities from $u=0$ to $u=N\in \mathbb{N}$ we get that
(8)
\[\begin{aligned}{}{\sum \limits_{u=0}^{N}}{\psi _{\delta }}(u)=& {\mathrm{e}^{-\delta }}{\sum \limits_{u=0}^{N}}{\overline{F}_{X}}(u)+{\mathrm{e}^{-2\delta }}{\sum \limits_{u=0}^{N}}{\sum \limits_{l=0}^{u}}{x_{l}}{\overline{F}_{Y}}(1+u-l)\\ {} & +{\mathrm{e}^{-2\delta }}\bigg({\sum \limits_{u=0}^{N}}{\sum \limits_{l=0}^{u}}{q_{l}}{\psi _{\delta }}(u+2-l)+{\psi _{\delta }}(1){\sum \limits_{u=0}^{N}}({q_{u+1}}-{x_{u+1}}{y_{0}})\bigg).\end{aligned}\]We observe that
and, similarly,
Hence, it follows from (8) that
for each $N\in \mathbb{N}$.
(9)
\[\begin{array}{l}\displaystyle {\sum \limits_{u=0}^{N+2}}{\psi _{\delta }}(u)\big(1-{\mathrm{e}^{-2\delta }}{F_{Q}}(N+2-u)\big)\\ {} \displaystyle \hspace{1em}={\mathrm{e}^{-\delta }}{\sum \limits_{u=0}^{N}}{\overline{F}_{X}}(u)+{\mathrm{e}^{-2\delta }}{\sum \limits_{u=1}^{N+1}}{\overline{F}_{Y}}(u){F_{X}}(N+1-u)\\ {} \displaystyle \hspace{2em}+{\psi _{\delta }}(N+1)+{\psi _{\delta }}(N+2)+{\mathrm{e}^{-2\delta }}{\psi _{\delta }}(1){\sum \limits_{u=0}^{N}}({q_{u+1}}-{x_{u+1}}{y_{0}})\\ {} \displaystyle \hspace{2em}-{\mathrm{e}^{-2\delta }}\big({\psi _{\delta }}(0){F_{Q}}(N+2)+{\psi _{\delta }}(1){F_{Q}}(N+1)\big)\end{array}\]Now we are in a position to let $N\to \infty $. It is obvious that:
Theorem 1 implies that
Consider the second term in the right side of equality (9). Obviously
Only the left side of equality (9) is left for consideration. Due to Theorem 1
(12)
\[ \underset{N\to \infty }{\lim }{\psi _{\delta }}(N+1)=\underset{N\to \infty }{\lim }{\psi _{\delta }}(N+2)=0.\]
\[ \underset{N\to \infty }{\lim }{\sum \limits_{u=1}^{N+1}}{\overline{F}_{Y}}(u){F_{X}}(N+1-u)\leqslant {\sum \limits_{u=1}^{\infty }}{\overline{F}_{Y}}(u).\]
On the other hand, for an arbitrary $M\in \mathbb{N}$
\[ \underset{N\to \infty }{\lim }{\sum \limits_{u=1}^{N+1}}{\overline{F}_{Y}}(u){F_{X}}(N+1-u)\geqslant \underset{N\to \infty }{\lim }{F_{X}}(N+1-M){\sum \limits_{u=1}^{M}}{\overline{F}_{Y}}(u)={\sum \limits_{u=1}^{M}}{\overline{F}_{Y}}(u).\]
Consequently,
(13)
\[\begin{aligned}{}\underset{N\to \infty }{\lim }{\sum \limits_{u=1}^{N+1}}{\overline{F}_{Y}}(u){F_{X}}(N+1-u)=& {\sum \limits_{u=1}^{\infty }}{\overline{F}_{Y}}(u)={y_{2}}+2{y_{3}}+3{y_{4}}+\dots \\ {} =& {y_{0}}+\mathbb{E}Y-1.\end{aligned}\]
\[ \underset{N\to \infty }{\lim }{\sum \limits_{u=0}^{N+2}}{\psi _{\delta }}(u)={\mathcal{S}_{\delta }}<\infty .\]
In addition,
\[ \underset{N\to \infty }{\lim }{\sum \limits_{u=0}^{N+2}}{\psi _{\delta }}(u){F_{Q}}(N+2-u)\leqslant {\sum \limits_{u=0}^{\infty }}{\psi _{\delta }}(u)={\mathcal{S}_{\delta }},\]
and, for an arbitrary chosen $M\in \mathbb{N}$,
Hence,
From this point we consider the three cases described in the formulation of Theorem separately.
(I) If ${q_{0}}>0$ then equality (15) implies that
where ${a_{1}}$, ${b_{1}}$ and ${d_{1}}$ are as defined in formulation of Theorem. So, we have that the main equality (1) holds if $n\in \{0,1\}$.
Now we need to prove this equality for all $n\in \mathbb{N}$. For this we use induction. Suppose that equality (1) holds for all $n\in \{0,1,\dots ,K\}$ for the defined sequences ${a_{n}}$, ${b_{n}}$ and ${d_{n}}$.
The induction hypothesis and equality (7) with $u=K-1$ imply that
\[\begin{aligned}{}{\mathrm{e}^{2\delta }}{\psi _{\delta }}(K-1)=& {\mathrm{e}^{2\delta }}\big({a_{K-1}}\hspace{0.1667em}{\psi _{\delta }}(0)+{b_{K-1}}\hspace{0.1667em}{\mathcal{S}_{\delta }}+{d_{K-1}}\big)\\ {} =& {\mathrm{e}^{\delta }}{\overline{F}_{X}}(K-1)+{\sum \limits_{l=0}^{K-1}}{x_{l}}{\overline{F}_{Y}}(K-l)+{q_{0}}\hspace{0.1667em}{\psi _{\delta }}(K+1)\\ {} & +{\sum \limits_{l=1}^{K-1}}{q_{l}}\big({a_{K+1-l}}\hspace{0.1667em}{\psi _{\delta }}(0)+{b_{K+1-l}}\hspace{0.1667em}{\mathcal{S}_{\delta }}+{d_{K+1-l}}\big)\\ {} & +({q_{K}}-{x_{K}}{y_{0}})\big({a_{1}}{\psi _{\delta }}(0)+{b_{1}}{\mathcal{S}_{\delta }}+{d_{1}}\big).\end{aligned}\]
Therefore,
\[\begin{aligned}{}{q_{0}}\hspace{0.1667em}{\psi _{\delta }}(K+1)=& {\psi _{\delta }}(0)\big({\mathrm{e}^{2\delta }}{a_{K-1}}-{\sum \limits_{l=1}^{K}}{q_{l}}{a_{K+1-l}}+{x_{K}}{y_{0}}{a_{1}}\big)\\ {} & +{\mathcal{S}_{\delta }}\big({\mathrm{e}^{2\delta }}{b_{K-1}}-{\sum \limits_{l=1}^{K}}{q_{l}}{b_{K+1-l}}+{x_{K}}{y_{0}}{b_{1}}\big)\\ {} & +\bigg({\mathrm{e}^{2\delta }}{d_{K-1}}-{\sum \limits_{l=1}^{K}}{q_{l}}{d_{K+1-l}}+{x_{K}}{y_{0}}{d_{1}}-{\mathrm{e}^{\delta }}{\overline{F}_{X}}(K-1)\\ {} & -{\sum \limits_{l=0}^{K-1}}{x_{l}}{\overline{F}_{Y}}(K-l)\bigg),\end{aligned}\]
or
\[ {\psi _{\delta }}(K+1)={a_{K+1}}{\psi _{\delta }}(0)+{b_{K+1}}{\mathcal{S}_{\delta }}+{d_{K+1}}\]
due to the definition of sequences ${a_{n}}$, ${b_{n}}$ and ${d_{n}}$.The induction principle implies that equality (1) holds for all $n\in {\mathbb{N}_{0}}$. The first part of Theorem 2 is proved.
(II) If ${x_{0}}=0$, ${y_{0}}\ne 0$, then equality (15) implies that
\[ {\psi _{\delta }}(1)={\tilde{a}_{1}}{\psi _{\delta }}(0)+{\tilde{b}_{1}}{\mathcal{S}_{\delta }}+{\tilde{d}_{1}}\]
where ${\tilde{a}_{1}}$, ${\tilde{b}_{1}}$ and ${\tilde{d}_{1}}$ are as defined in formulation of Theorem. So, we have that the main equality (2) holds if $n\in \{0,1\}$. Similarly as in case (I), we finish the proof using induction method and equality (7).
(III) If ${x_{0}}\ne 0$, ${y_{0}}=0$, then equality (15) implies that
where ${\hat{b}_{0}}$ and ${\hat{d}_{0}}$ are as defined in formulation of Theorem. So, we have that the main equality (3) holds if $n=0$. Similarly as in case (I), we finish the proof using induction method and equality (7).
Now Theorem 2 is proved. □
6 Bounds for the Gerber–Shiu Discounted Penalty Function
Let us consider the bi-seasonal risk model ${W_{u}}(n)={W_{u}}(n|X,Y)$ defined in Section 1. Let ${T_{u}}={T_{u}}(X,Y)$ be the time of ruin for the model.
Also, let us consider homogeneous discrete time risk model generated by claim amount $(X+Y)/2$. Suppose that $\hat{{T_{u}}}((X+Y)/2)$ denotes the ruin time of this model, i.e.
\[ \hat{{T_{u}}}((X+Y)/2)=\left\{\begin{array}{l}\min \{n\geqslant 1:{\hat{W}_{u}}(n)\leqslant 0\},\\ {} \infty ,\hspace{1em}\text{if}\hspace{2.5pt}{\hat{W}_{u}}(n)>0\hspace{2.5pt}\text{for all}\hspace{2.5pt}\hspace{2.5pt}n\in \mathbb{N},\end{array}\right.\]
where $\hat{{W_{u}}}(n)=u+n-{\textstyle\sum _{i=1}^{n}}{\hat{Z}_{i}}$ with independent r.v.s. ${\hat{Z}_{1}}$, ${\hat{Z}_{2}}$, … distributed as $(X+Y)/2$.Then we have that ${W_{u}}(n|X,Y)\geqslant {W_{u}}(n|X+Y,0)$, and hence
\[\begin{aligned}{}{T_{u}}(X,Y)=& \min \big\{n\in \mathbb{N}:{W_{u}}(n|X,Y)\leqslant 0\big\}\\ {} \geqslant & \min \big\{n\in \mathbb{N}:{W_{u}}(n|X+Y,0)\leqslant 0\big\}\\ {} =& \min \bigg\{n\in 2\mathbb{N}-1:u+n-{\sum \limits_{i=1}^{[(n+1)/2]}}({X_{i}}+{Y_{i}})\leqslant 0\bigg\}\\ {} =& 2\min \bigg\{k\in \mathbb{N}:u+2k-1-{\sum \limits_{i=1}^{k}}({X_{i}}+{Y_{i}})\leqslant 0\bigg\}-1\\ {} =& 2{\hat{T}_{(u-1)/2}}\big((X+Y)/2\big)-1.\end{aligned}\]
Thus we get
\[\begin{aligned}{}{\psi _{\delta }}(u)=& {\hat{\psi }_{\delta }}\big({T_{u}}(X,Y)\big)\leqslant {\hat{\psi }_{\delta }}\big(2{\hat{T}_{(u-1)/2}}\big((X+Y)/2\big)-1\big)\\ {} =& {\mathrm{e}^{\delta }}{\hat{\psi }_{2\delta }}\big({\hat{T}_{(u-1)/2}}\big((X+Y)/2\big)\big),\end{aligned}\]
where
\[ {\hat{\psi }_{\Delta }}(T)=\mathbb{E}\big({\mathrm{e}^{-\Delta T}}{\mathbb{1}_{\{T<\infty \}}}\big)\]
for a r.v. T and arbitrary $\Delta >0$. We also have ${W_{u}}(n|X,Y)\leqslant {W_{u}}(n|0,X+Y).$ Hence, similarly as above we obtain
\[ {\psi _{\delta }}(u)\geqslant {\hat{\psi }_{2\delta }}\big({\hat{T}_{u/2}}\big((X+Y)/2\big)\big).\]
To summarize, upper and lower bounds for ${\psi _{\delta }}(u)$ were obtained:
\[ {\hat{\psi }_{2\delta }}\big({\hat{T}_{u/2}}((X+Y)/2)\big)\leqslant {\psi _{\delta }}(u)\leqslant {\mathrm{e}^{\delta }}{\hat{\psi }_{2\delta }}\big({\hat{T}_{(u-1)/2}}\big((X+Y)/2\big)\big)\]
for all $\delta >0$ and $u\geqslant 1$.7 Algorithm Code in R
library(Rmpfr)
# set the values of parameters
delta = 0.1; N = 30; K = 2; umax = 20
# initialise vectors
q = numeric(N)
a = mpfrArray(0, precBits = 1024, dim = c(N,1))
b = mpfrArray(0, precBits = 1024, dim = c(N,1))
d = mpfrArray(0, precBits = 1024, dim = c(N,1))
psi = mpfrArray(0, precBits = 1024, dim = c((umax+1),1))
FX = numeric(N)
FY = numeric(N)
# choose the distributions of claims (4 different distributions are
considered as described in Numerical examples section)
x = c(0.6, 0.2, 0.2); y = c(0.5, 0.2, 0.2, 0.1)
# x = c(0.4,0.6); y = c(0.1,0.6,0.3)
# y = c(0.4,0.6); x = c(0.1,0.6,0.3)
# lambda = 0.8; prob = 0.7; x = dpois(c(0:N),lambda);
y = dgeom(c(0:N),prob)
# compute quantities related with claims’ distributions
Xmax <- length(x)-1; Ymax <- length(y)-1
X = 0:Xmax; Y = 0:Ymax
EX = sum(X * x); EY = sum(Y * y)
x[(Xmax+2):N] = 0; y[(Ymax+2):N] = 0
for (i in 0:(Xmax+Ymax)) {
for (k in 1:(i+1))
q[i+1] = q[i+1] + x[k] * y[(i+2) - k] }
FX[1] = x[1]
for (u in 1:(N-1)) {FX[u+1] = FX[u] + x[u+1]}
F_X = 1 - FX
FY[1] = y[1]
for (u in 1:(N-1)) {FY[u+1] = FY[u] + y[u+1]}
F_Y = 1 - FY
# calculate the coefficients of algorithm
a[1] = mpfr(1,1024)
a[2] = mpfr(-1,1024) / mpfr(y[1],1024)
for (n in 2:(N-1)) {
a[n + 1] = mpfr(1,1024) / mpfr(q[1],1024) * (mpfr(exp(2 * delta),1024)
* mpfr(a[(n + 1) - 2],1024) + mpfr(x[n],1024) * mpfr(y[1],1024)
* mpfr(a[2],1024))
for (i in 2:n)
a[n + 1] = mpfr(a[n + 1],1024) - (mpfr(1,1024) / mpfr(q[1],1024))
* (mpfr(q[i],1024) * mpfr(a[n - i + 2],1024)) }
b[1] = 0
b[2] = -(exp(2 * delta) - 1) / mpfr(y[1],1024)
for (n in 2:(N-1)) {
b[n + 1] = 1 / mpfr(q[1],1024) * (exp(2 * delta)
* mpfr(b[(n + 1) - 2],1024) + mpfr(x[n],1024) * mpfr(y[1],1024)
* mpfr(b[2],1024))
for (i in 2:n)
b[n + 1] = mpfr(b[n + 1],1024) - (1 / mpfr(q[1],1024))
* (mpfr(q[i],1024) * mpfr(b[n - i + 2],1024)) }
d[1] = 0
d[2] = (exp(delta) * mpfr(EX,1024) + mpfr(y[1],1024) + mpfr(EY,1024) - 1)
/ mpfr(y[1],1024)
for (n in 2:(N-1)) {
d[n + 1] = 1 / mpfr(q[1],1024) * (exp(2 * delta)
* mpfr(d[(n + 1) - 2],1024)
+ mpfr(x[n],1024) * mpfr(y[1],1024) * mpfr(d[2],1024)
- exp(delta) * mpfr(F_X[n - 1],1024))
for (i in 2:n)
d[n + 1] = mpfr(d[n + 1],1024) - (1 / mpfr(q[1],1024))
* (mpfr(q[i],1024) * mpfr(d[n - i + 2],1024)
+ mpfr(x[i - 1],1024) * mpfr(F_Y[n - i + 2],1024)) }
# solve the system of linear equations
eqA = array(c(mpfr(a[N-K],1024), mpfr(a[N],1024), mpfr(b[N-K],1024),
mpfr(b[N],1024)),
dim = c(2, 2))
eqb = array(c(mpfr(-d[N-K],1024), mpfr(-d[N],1024)))
detA = mpfr(eqA[1,1],1024) * mpfr(eqA[2,2],1024)
- mpfr(eqA[1,2],1024) * mpfr(eqA[2,1],1024)
eqA_inv = 1/detA * array(c(mpfr(eqA[2,2],1024), mpfr(-eqA[2,1],1024),
mpfr(-eqA[1,2],1024), mpfr(eqA[1,1],1024)), dim = c(2, 2))
eqx = mpfr(eqA_inv,1024) %*% mpfr(eqb,1024)
id_mat = mpfr(eqA_inv,1024) %*% mpfr(eqA,1024)
psi[1] = eqx[1]
S = eqx[2]
# check the accuracy of solutions
acc_psi0 = mpfr(exp(-delta),1024) * (abs(mpfr(b[N-K],1024))
+ abs(mpfr(b[N],1024))) / abs(mpfr(detA,1024))
acc_S = exp(-delta) * (abs(a[N-K]) + abs(a[N])) / abs(detA)
# calculate the values of Gerber--Shiu function
psi[2] = a[2] * psi[1] + b[2] * S + d[2]
psi[3:(umax+1)] = a[3:(umax+1)] * psi[1] + b[3:(umax+1)] * S
+ d[3:(umax+1)]
psi2 = asNumeric(psi)
