INFORMATICA Informatica 1822-88440868-49520868-4952 Vilnius University INFO1202 10.15388/Informatica.2018.190 Research Article The Gerber–Shiu Discounted Penalty Function for the Bi-Seasonal Discrete Time Risk Model NavickienėOlga

O. Navickienė is a PhD student at Institute of Mathematics of Vilnius University. Her research interests include probability theory, actuarial mathematics and economics.

 SprindysJonas

J. Sprindys is a PhD student at Institute of Mathematics of Vilnius University. His research interests include probability theory, risk theory and actuarial mathematics.

 ŠiaulysJonasjonas.siaulys@mif.vu.lt∗

J. Šiaulys is a professor at Institute of Mathematics of Vilnius University. His fields of interest are number theory, probability theory, risk theory and actuarial mathematics. He has authored or co-authored around 80 papers.

 Institute of Mathematics, Vilnius University, Naugarduko 24, Vilnius LT-03225, Lithuania Corresponding author. 201811201829473375632018102018 © 2018 Vilnius University2018 Open access article under the CC BY license.

In this work, the discrete time risk model with two seasons is considered. In such model, the claims repeat with time periods of two units, i.e. claim distributions coincide at all even instants and at all odd instants. Our purpose is to derive an algorithm for calculating the values of the particular case of the Gerber–Shiu discounted penalty function E ( e − δ T 1 { T < ∞ } ), where T is the time of ruin, and δ is a constant nonnegative force of interest. Theoretical results are illustrated by some numerical examples.

 bi-seasonal model discrete time risk model Gerber–Shiu function penalty function; time of ruin 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.

We say that the insurer’s surplus W u varies according to the bi-seasonal risk model if W u ( n ) = u + n − ∑ i = 1 n Z i for each n ∈ N 0 = { 0 , 1 , 2 , … } and the following assumptions hold:

the initial insurer’s surplus u ∈ N 0 ,

the random claim amounts { Z 1 , Z 2 , … } are nonnegative integer-valued independent r.v.s.,

there exist r.v.s. X and Y such that Z 2 k + 1 = d X, k ∈ N 0 , and Z 2 k = d Y, k ∈ N.

If X = 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 Ψ δ , 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 Ψ δ , w ( u ) = E ( e − δ T u w ( W u ( T u − 1 ) , | W u ( T u ) | ) 1 { T u < ∞ } ) , where force of interest δ ⩾ 0, w ( x , y ) is an arbitrary function of two nonnegative arguments, and T u denotes the time of ruin, i.e. T u = min { n ⩾ 1 : W u ( n ) ⩽ 0 } , ∞ , if W u ( n ) > 0 for all n ∈ N . 0\hspace{2.5pt}\text{for all}\hspace{2.5pt}n\in \mathbb{N}.\end{array}\right.\]]]>

Function w has practical interpretations. For example, if w was interpreted as the benefit amount of reinsurance payable at the time of ruin, then Ψ δ , 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 ψ δ ( u ) = Ψ δ , 1 ( u ) = E ( e − δ T u 1 { T u < ∞ } ) .

If, in addition, force of interest δ = 0, then the Gerber–Shiu discounted penalty function is equal to the ruin probability ψ ( u ) = ψ 0 ( u ) = Ψ 0 , 1 ( u ) = P ( T u < ∞ ) .

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 Ψ δ , w of the homogeneous discrete time risk model can be calculated using the following formulas Ψ δ , w ( 0 ) = e − δ ∑ k = 0 ∞ ∑ l = 0 ∞ ϱ k w ( k , l ) P ( Z = k + l + 1 ) , Ψ δ , w ( u ) = e − δ ∑ k = 0 u − 1 Ψ δ , w ( u − k ) ∑ l = 0 ∞ ϱ l P ( Z = k + l + 1 ) + e − δ ϱ − u ∑ k = u ∞ ϱ k ∑ l = 0 ∞ w ( k , l ) P ( Z = k + l + 1 ) , where Z with E Z < 1 is the integer-valued random variable generating the homogeneous discrete time risk model, and ϱ ∈ ( 0 , 1 ) is the root of equation s e δ = ∑ k = 0 ∞ s k P ( Z = k ) . By arguments provided in Li and Garrido (2002), such a solution exists and is unique for δ > 00$]]>.

If the discrete time risk model is generated by possibly differently distributed random variables Z 1 , Z 2 , … , 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 ψ δ . 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 = P ( X = k ) , y k = P ( Y = k ) , q k = P ( Q = k ) , k ∈ 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. F X ( u ) = P ( X ⩽ u ) = ∑ k = 0 ⌊ u ⌋ x k , F Y ( u ) = P ( Y ⩽ u ) = ∑ k = 0 ⌊ u ⌋ y k , F Q ( u ) = P ( Q ⩽ u ) = ∑ k = 0 ⌊ u ⌋ q k , for each real u. The notation F ‾ is used for the tail of an arbitrary distribution function F, i.e. F ‾ ( u ) = 1 − F ( u ) for each u ∈ R.

The following two assertions enable us to construct an algorithm for calculating values of function ψ δ ( u ) in the bi-seasonal discrete time risk model.

Let the bi-seasonal discrete time risk model be generated by two nonnegative, independent and integer valued random variables X and Y. If E X + E Y < 2, then lim u → ∞ ψ δ ( u ) = 0 for an arbitrary fixed δ ⩾ 0. In addition, if max { E e h X , E e h Y } < ∞ for some positive h, then ∑ l = 0 ∞ ψ δ ( l ) < ∞ for each fixed δ ⩾ 0.

Let all the conditions of Theorem 1 be satisfied. Furthermore, let δ > 00$]]>, and ψ δ denote the Gerber–Shiu function with w ( x , y ) = 1 for all nonnegative x and y. Also denote S δ : = ∑ l = 0 ∞ ψ δ ( l ).

If q 0 = P ( X + Y = 0 ) > 00$]]>, then ψ δ ( n ) = a n ψ δ ( 0 ) + b n S δ + d n

for each n ∈ N 0 , where a n , b n , d n are three sequences of real numbers defined recursively by the following equations: a 0 = 1 , a 1 = − 1 y 0 , a n = 1 q 0 ( e 2 δ a n − 2 − ∑ i = 1 n − 1 q i a n − i − x n − 1 ) , n ∈ { 2 , 3 , … } ; b 0 = 0 , b 1 = − e 2 δ − 1 y 0 , b n = 1 q 0 ( e 2 δ b n − 2 − ∑ i = 1 n − 1 q i b n − i − x n − 1 ( e 2 δ − 1 ) ) , n ∈ { 2 , 3 , … } ; d 0 = 0 , d 1 = e δ E X + y 0 + E Y − 1 y 0 , d n = 1 q 0 ( e 2 δ d n − 2 − ∑ i = 1 n − 1 q i d n − i + x n − 1 y 0 d 1 − e δ F ‾ X ( n − 2 ) − ∑ i = 0 n − 2 x i F ‾ Y ( n − 1 − i ) ) , n ∈ { 2 , 3 , … } .

If x 0 = P ( X = 0 ) = 0 and y 0 = P ( Y = 0 ) ≠ 0, then ψ δ ( n ) = a ˜ n ψ δ ( 0 ) + b ˜ n S δ + d ˜ n

for each n ∈ N 0 , where a ˜ n , b ˜ n , d ˜ n are three sequences of real numbers defined recursively by the following equations: a ˜ 0 = 1 , a ˜ 1 = − 1 y 0 , a ˜ n = 1 q 1 ( e 2 δ a ˜ n − 1 − ∑ i = 1 n − 1 q i + 1 a ˜ n − i − x n ) , n ∈ { 2 , 3 , … } ; b ˜ 0 = 0 , b ˜ 1 = − e 2 δ − 1 y 0 , b ˜ n = 1 q 1 ( e 2 δ b ˜ n − 1 − ∑ i = 1 n − 1 q i + 1 b ˜ n − i − x n ( e 2 δ − 1 ) ) , n ∈ { 2 , 3 , … } ; d ˜ 0 = 0 , d ˜ 1 = e δ E X + y 0 + E Y − 1 y 0 , d ˜ n = 1 q 1 ( e 2 δ d ˜ n − 1 − ∑ i = 1 n − 1 q i + 1 d ˜ n − i + x n y 0 d ˜ 1 − e δ F ‾ X ( n − 1 ) d ˜ n = − ∑ i = 0 n − 2 x i + 1 F ‾ Y ( n − 1 − i ) ) , n ∈ { 2 , 3 , … } .

If x 0 ≠ 0 and y 0 = 0, then ψ δ ( n ) = b ˆ n S δ + d ˆ n

for each n ∈ N 0 , where b ˆ n , d ˆ n are two sequences of real numbers defined recursively by the following equations: b ˆ 0 = − ( e 2 δ − 1 ) , b ˆ n = 1 q 1 ( e 2 δ b ˆ n − 1 − ∑ i = 1 n − 1 q i + 1 b ˆ n − i ) , n ∈ N ; d ˆ 0 = e δ E X + E Y − 1 , d ˆ n = 1 q 1 ( e 2 δ d ˆ n − 1 − ∑ i = 1 n − 1 q i + 1 d ˆ n − i − e δ F ‾ X ( n − 1 ) − ∑ i = 0 n − 1 x i F ‾ Y ( n − i ) ) , n ∈ N .

We observe that case x 0 = y 0 = 0 is impossible due to requirement 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.

 Algorithm for Finding the Values of Function <inline-formula id="j_info1202_ineq_060"><alternatives><mml:math> <mml:msub> <mml:mrow> <mml:mi mathvariant="italic">ψ</mml:mi> </mml:mrow> <mml:mrow> <mml:mi mathvariant="italic">δ</mml:mi> </mml:mrow> </mml:msub></mml:math><tex-math><![CDATA[${\psi _{\delta }}$]]></tex-math></alternatives></inline-formula>

In this section, we describe an algorithm for calculating values of ψ δ ( 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 = P ( X = k ), y k = P ( Y = k ), k ∈ 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 ψ δ ( u ), u ∈ N 0 in the case when x 0 y 0 > 00$]]>. The other possible cases: { x 0 = 0 , y 0 > 0 }0\}$]]>, and { x 0 > 0 , y 0 = 0 }0,{y_{0}}=0\}$]]>, which were described in Theorem 2, can be considered similarly.

Step 1: Select N ∈ { 10 , 20 , 30 , … , 100 } and K ∈ { 1 , … , 5 }.

Step 2: Calculate coefficients a n , b n , d n for all n ∈ { 0 , 1 , … , N } using formulas from Theorem 2.

Step 3: Find ψ ˆ δ ( 0 ) and S ˆ δ satisfying the following system of linear equations a N − K ψ ˆ δ ( 0 ) + b N − K S ˆ δ + d N − K = 0 , a N ψ ˆ δ ( 0 ) + b N S ˆ δ + d N = 0 .

Due to the main formula (1) of Theorem 2 the desired quantity ψ δ ( 0 ) together with sum S δ satisfy the following system a N − K ψ δ ( 0 ) + b N − K S δ + d N − K = ψ δ ( N − K ) , a N ψ δ ( 0 ) + b N S δ + d N = ψ δ ( N ) . However, according to Theorem 1 ψ δ ( N − K ) and ψ δ ( N ) are close to zero for sufficiently large N. We get system (4) from (5) by changing values of ψ δ ( N − K ) and ψ δ ( N ) to zeroes.

Step 4: Test the error | ψ δ ( 0 ) − ψ ˆ δ ( 0 ) |.

Using the Cramer’s rule for both systems of linear equations (4), (5) and the trivial estimate | ψ δ ( n ) | ⩽ 1, n ∈ N 0 , we derive that | ψ δ ( 0 ) − ψ ˆ δ ( 0 ) | ⩽ e − δ ( | b N − K | + | b N | ) | a N − K b N − b N − K a N | . Numerical simulations have shown that the upper estimate of ψ δ ( 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 ψ δ ( 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 ψ δ ( 0 ) when changing system (5) to (4), since ψ δ ( N ) does not converge to zero so quickly. As for parameter K, it should be chosen to minimize the upper estimate of ψ δ ( 0 ) approximation error.

Step 6: Calculate ψ δ ( 1 ) according to the formula (1) by supposing that ψ δ ( 0 ) = ψ ˆ δ ( 0 ) and S δ = S ˆ δ .

Step 7: Calculate values of ψ δ ( u ) for u ⩾ 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, ψ δ ( 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 ∈ 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.

Also, computational errors could arise because by using formula (1) from Theorem 2, we are calculating “small” quantity ψ δ ( u ) as a sum containing “large” in absolute value summands.

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 δ = 0 and w ( x , y ) = 1, was considered. In the paper we have extended the results to the case δ > 00$]]>.

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 ψ δ ( 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.

 Numerical Examples

In this section, we present four numerical examples for calculating the values of ψ δ ( u ), u ∈ N 0 , in the bi-seasonal discrete time risk model. In all examples we consider function ψ δ with three different values of the interest force δ ∈ { 0 ; 0.01 ; 0.1 }. Our algorithm does not allow to compute function values for case δ = 0, so the algorithm and its ψ δ ( 0 ) approximation error upper estimate provided in Damarackas and Šiaulys (2014) were used for this case. Since the function ψ δ ( u ) seems to decay exponentially, all the figures are plotted in log scale (with base 10).

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 { E X , E Y } < 1, and all conditions of Theorem 2 are satisfied. Using the algorithm presented in Section 2 we obtain values of ψ δ ( u ) for u ∈ { 0 , 1 , … , 15 }. These values are presented in Table 1 and are shown in (Fig. 1). The upper estimate of ψ δ ( 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 ψ δ ( 0 ) which is obtained with N = 50 and K = 2.

Values of ψ δ ( u ) in Example 1.

u δ = 0 (<0.000000001) δ = 0.01 (0.083494161) δ = 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

Values of ψ δ ( u ) in Example 1 (log scale).

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 E X < 1, E Y ⩾ 1, but E X + 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 ψ δ ( u ) for u ∈ { 0 , 1 , … , 15 }. Results of this example are based on the value of ψ δ ( 0 ) which is obtained with N = 40 and K = 1. Values of ψ δ ( u ) are also shown in (Fig. 2).

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 ψ δ ( u ) are also changing. The numerical values of ψ δ ( u ) of this model are given in the Table 2 and shown in (Fig. 2) with N = 50 and K = 3.

Values of ψ δ ( u ) in Examples 2 and 3

u Example 2 Example 2 Example 2 Example 3 Example 3 Example 3
δ = 0 δ = 0.01 δ = 0.1 δ = 0 δ = 0.01 δ = 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

Values of ψ δ ( u ) in Examples 2 and 3 (log scale).

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 λ = 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 ψ δ ( 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.

Values of ψ δ ( u ) in Example 4.

u δ = 0 (<0.000000001) δ = 0.01 (0.089541014) δ = 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

Values of ψ δ ( u ) in Example 4 (log scale),

 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 ψ δ ( 0 ) approximation error could be derived.

In the bi-seasonal discrete time risk model, claims with distributions satisfying E X + E Y ⩾ 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 ψ δ ( 0 ) and ψ δ ( 1 ) should be derived.

Proofs of the Main Results

According to Theorem 2.3 of (Damarackas and Šiaulys, 2014), we have that lim u → ∞ ψ ( u ) = 0 . This implies that lim u → ∞ ψ δ ( u ) = 0 for an arbitrary fixed δ ⩾ 0, because 0 ⩽ ψ δ ( u ) ⩽ ψ ( u ) for all δ , u ⩾ 0.

Furthermore, for i ∈ N denote η i = Z i − 1. The conditions of Theorem 1 imply that sup i ∈ N E ( e h η i ) = max { E ( e h X ) , E ( e h Y ) } < ∞ , lim u → ∞ sup i ∈ N E ( | η i | 1 { η i ⩽ − u } ) = lim u → ∞ max { E ( ( 1 − X ) 1 { X ⩽ 1 − u } ) , E ( ( 1 − Y ) 1 { Y ⩽ 1 − u } ) } = 0 , lim sup n → ∞ 1 n ∑ i = 1 n E η i = E X + E Y − 2 2 < 0 .

Hence, according to Lemma 1 by (Andrulytė et al., 2015), we have ψ δ ( u ) ⩽ ψ ( u ) = P ( sup k ⩾ 1 ∑ i = 1 k η i > u ) ⩽ c 1 e − c 2 u , u ⩾ 0 , 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 ∑ l = 0 ∞ ψ δ ( l ) < ∞ for each fixed δ ⩾ 0. □

Suppose that δ > 00$]]> and u ∈ N 0 . According to definition of function ψ δ we have ψ δ ( u ) = ∑ m = 1 ∞ E ( e − δ m 1 { T u = m } ) = ∑ m = 1 ∞ e − δ m P ( ∑ i = 1 j Z i < j + u for j ∈ { 1 , 2 , … , m − 1 } and ∑ i = 1 m Z i ⩾ m + u ) = e − δ P ( Z 1 ⩾ 1 + u ) + e − 2 δ P ( Z 1 < 1 + u , Z 1 + Z 2 ⩾ 2 + u ) + ∑ m = 3 ∞ e − δ m P ( ∑ i = 1 j Z i < j + u for j ∈ { 1 , 2 , … , m − 1 } and ∑ i = 1 m Z i ⩾ m + u ) .

Since X = d Z 1 = d Z 3 = d Z 5 = d … and Y = d Z 2 = d Z 4 = d Z 6 = d … we get that ψ δ ( u ) = e − δ ∑ l ⩾ 1 + u x l + e − 2 δ ∑ l ⩽ u ∑ k ⩾ 2 + u − l x l y k + ∑ m = 3 ∞ e − δ m P ( Z 1 ⩽ u , Z 1 + Z 2 ⩽ 1 + u , Z 1 + Z 2 + ∑ i = 3 j Z i < j + u for j ∈ { 3 , … , m − 1 } and Z 1 + Z 2 + ∑ i = 3 m Z i ⩾ m + u ) = e − δ F ‾ X ( u ) + e − 2 δ ∑ l = 0 u x l F ‾ Y ( 1 + u − l ) + ∑ l = 0 u ∑ k = 0 1 + u − l x l y k ∑ m = 3 ∞ e − δ m P ( ∑ i = 3 j Z i < j + u − l − k for j ∈ { 3 , … , m − 1 } and ∑ i = 3 m Z i ⩾ m + u − k − l ) = e − δ F ‾ X ( u ) + e − 2 δ ∑ l = 0 u x l F ‾ Y ( 1 + u − l ) + e − 2 δ ∑ l = 0 u ∑ k = 0 1 + u − l x l y k ∑ m = 3 ∞ e − δ ( m − 2 ) P ( ∑ i = 1 j Z i < j + u − l − k for j ∈ { 1 , … , m − 3 } and ∑ i = 1 m − 2 Z i ⩾ m + u − k − l ) = e − δ F ‾ X ( u ) + e − 2 δ ∑ l = 0 u x l F ‾ Y ( 1 + u − l ) + e − 2 δ ∑ l = 0 u ∑ k = 0 1 + u − l x l y k ψ δ ( u + 2 − k − l ) .

For each m ∈ N 0 q m = P ( Q = m ) = ∑ k = 0 m x k y m − k .

Therefore the last sum in equality (6) can be expressed by ∑ l = 0 1 + u ∑ k = 0 1 + u − l x l y k ψ δ ( u + 2 − ( k + l ) ) − x u + 1 y 0 ψ δ ( 1 ) = ∑ l = 0 1 + u q l ψ δ ( u + 2 − l ) − x u + 1 y 0 ψ δ ( 1 ) = ∑ l = 0 u q l ψ δ ( u + 2 − l ) + ( q u + 1 − x u + 1 y 0 ) ψ δ ( 1 ) .

Substituting this expression into equality (6) we obtain that ψ δ ( u ) = e − δ F ‾ X ( u ) + e − 2 δ ∑ l = 0 u x l F ‾ Y ( 1 + u − l ) + e − 2 δ ( ∑ l = 0 u q l ψ δ ( u + 2 − l ) + ( q u + 1 − x u + 1 y 0 ) ψ δ ( 1 ) ) for each u ∈ N 0 .

By summing these last equalities from u = 0 to u = N ∈ N we get that ∑ u = 0 N ψ δ ( u ) = e − δ ∑ u = 0 N F ‾ X ( u ) + e − 2 δ ∑ u = 0 N ∑ l = 0 u x l F ‾ Y ( 1 + u − l ) + e − 2 δ ( ∑ u = 0 N ∑ l = 0 u q l ψ δ ( u + 2 − l ) + ψ δ ( 1 ) ∑ u = 0 N ( q u + 1 − x u + 1 y 0 ) ) .

We observe that ∑ u = 0 N ∑ l = 0 u x l F ‾ Y ( 1 + u − l ) = ∑ u = 1 N + 1 F ‾ Y ( u ) F X ( N + 1 − u ) and, similarly, ∑ u = 0 N ∑ l = 0 u q l ψ δ ( u + 2 − l ) = ∑ u = 2 N + 2 ψ δ ( u ) F Q ( N + 2 − u ) .

Hence, it follows from (8) that ∑ u = 0 N + 2 ψ δ ( u ) ( 1 − e − 2 δ F Q ( N + 2 − u ) ) = e − δ ∑ u = 0 N F ‾ X ( u ) + e − 2 δ ∑ u = 1 N + 1 F ‾ Y ( u ) F X ( N + 1 − u ) + ψ δ ( N + 1 ) + ψ δ ( N + 2 ) + e − 2 δ ψ δ ( 1 ) ∑ u = 0 N ( q u + 1 − x u + 1 y 0 ) − e − 2 δ ( ψ δ ( 0 ) F Q ( N + 2 ) + ψ δ ( 1 ) F Q ( N + 1 ) ) for each N ∈ N.

Now we are in a position to let N → ∞. It is obvious that: lim N → ∞ ∑ u = 0 N F ‾ X ( u ) = E X , lim N → ∞ F Q ( N + 1 ) = lim N → ∞ F Q ( N + 2 ) = 1 , lim N → ∞ ∑ u = 0 N q u + 1 = 1 − q 0 , lim N → ∞ ∑ u = 0 N + 1 x u + 1 = 1 − x 0 . Theorem 1 implies that lim N → ∞ ψ δ ( N + 1 ) = lim N → ∞ ψ δ ( N + 2 ) = 0 . Consider the second term in the right side of equality (9). Obviously lim N → ∞ ∑ u = 1 N + 1 F ‾ Y ( u ) F X ( N + 1 − u ) ⩽ ∑ u = 1 ∞ F ‾ Y ( u ) . On the other hand, for an arbitrary M ∈ N lim N → ∞ ∑ u = 1 N + 1 F ‾ Y ( u ) F X ( N + 1 − u ) ⩾ lim N → ∞ F X ( N + 1 − M ) ∑ u = 1 M F ‾ Y ( u ) = ∑ u = 1 M F ‾ Y ( u ) . Consequently, lim N → ∞ ∑ u = 1 N + 1 F ‾ Y ( u ) F X ( N + 1 − u ) = ∑ u = 1 ∞ F ‾ Y ( u ) = y 2 + 2 y 3 + 3 y 4 + … = y 0 + E Y − 1 . Only the left side of equality (9) is left for consideration. Due to Theorem 1 lim N → ∞ ∑ u = 0 N + 2 ψ δ ( u ) = S δ < ∞ . In addition, lim N → ∞ ∑ u = 0 N + 2 ψ δ ( u ) F Q ( N + 2 − u ) ⩽ ∑ u = 0 ∞ ψ δ ( u ) = S δ , and, for an arbitrary chosen M ∈ N, lim N → ∞ ∑ u = 0 N + 2 ψ δ ( u ) F Q ( N + 2 − u ) ⩾ lim N → ∞ F Q ( N + 2 − M ) ∑ u = 0 M ψ δ ( u ) . Hence, lim N → ∞ ∑ u = 0 N + 2 ψ δ ( u ) ( 1 − e − 2 δ F Q ( N + 2 − u ) ) = ( 1 − e − 2 δ ) S δ .

Substituting all limiting relations (10)–(14) into equality (9) we get ( 1 − e − 2 δ ) S δ = e − δ E X + e − 2 δ ( y 0 + E Y − 1 ) − e − 2 δ ψ δ ( 1 ) y 0 − e − 2 δ ψ δ ( 0 ) .

From this point we consider the three cases described in the formulation of Theorem separately.

(I) If q 0 > 00$]]> then equality (15) implies that ψ δ ( 1 ) = a 1 ψ δ ( 0 ) + b 1 S δ + d 1 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 ∈ { 0 , 1 }.

Now we need to prove this equality for all n ∈ N. For this we use induction. Suppose that equality (1) holds for all n ∈ { 0 , 1 , … , K } for the defined sequences a n , b n and d n .

The induction hypothesis and equality (7) with u = K − 1 imply that e 2 δ ψ δ ( K − 1 ) = e 2 δ ( a K − 1 ψ δ ( 0 ) + b K − 1 S δ + d K − 1 ) = e δ F ‾ X ( K − 1 ) + ∑ l = 0 K − 1 x l F ‾ Y ( K − l ) + q 0 ψ δ ( K + 1 ) + ∑ l = 1 K − 1 q l ( a K + 1 − l ψ δ ( 0 ) + b K + 1 − l S δ + d K + 1 − l ) + ( q K − x K y 0 ) ( a 1 ψ δ ( 0 ) + b 1 S δ + d 1 ) . Therefore, q 0 ψ δ ( K + 1 ) = ψ δ ( 0 ) ( e 2 δ a K − 1 − ∑ l = 1 K q l a K + 1 − l + x K y 0 a 1 ) + S δ ( e 2 δ b K − 1 − ∑ l = 1 K q l b K + 1 − l + x K y 0 b 1 ) + ( e 2 δ d K − 1 − ∑ l = 1 K q l d K + 1 − l + x K y 0 d 1 − e δ F ‾ X ( K − 1 ) − ∑ l = 0 K − 1 x l F ‾ Y ( K − l ) ) , or ψ δ ( K + 1 ) = a K + 1 ψ δ ( 0 ) + b K + 1 S δ + 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 ∈ N 0 . The first part of Theorem 2 is proved.

(II) If x 0 = 0, y 0 ≠ 0, then equality (15) implies that ψ δ ( 1 ) = a ˜ 1 ψ δ ( 0 ) + b ˜ 1 S δ + d ˜ 1 where a ˜ 1 , b ˜ 1 and d ˜ 1 are as defined in formulation of Theorem. So, we have that the main equality (2) holds if n ∈ { 0 , 1 }. Similarly as in case (I), we finish the proof using induction method and equality (7).

(III) If x 0 ≠ 0, y 0 = 0, then equality (15) implies that ψ δ ( 0 ) = b ˆ 0 S δ + d ˆ 0 , where b ˆ 0 and 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. □

 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 T u ˆ ( ( X + Y ) / 2 ) denotes the ruin time of this model, i.e. T u ˆ ( ( X + Y ) / 2 ) = min { n ⩾ 1 : W ˆ u ( n ) ⩽ 0 } , ∞ , if W ˆ u ( n ) > 0 for all n ∈ N , 0\hspace{2.5pt}\text{for all}\hspace{2.5pt}\hspace{2.5pt}n\in \mathbb{N},\end{array}\right.\]]]> where W u ˆ ( n ) = u + n − ∑ i = 1 n Z ˆ i with independent r.v.s. Z ˆ 1 , Z ˆ 2 , … distributed as ( X + Y ) / 2.

Then we have that W u ( n | X , Y ) ⩾ W u ( n | X + Y , 0 ), and hence T u ( X , Y ) = min { n ∈ N : W u ( n | X , Y ) ⩽ 0 } ⩾ min { n ∈ N : W u ( n | X + Y , 0 ) ⩽ 0 } = min { n ∈ 2 N − 1 : u + n − ∑ i = 1 [ ( n + 1 ) / 2 ] ( X i + Y i ) ⩽ 0 } = 2 min { k ∈ N : u + 2 k − 1 − ∑ i = 1 k ( X i + Y i ) ⩽ 0 } − 1 = 2 T ˆ ( u − 1 ) / 2 ( ( X + Y ) / 2 ) − 1 . Thus we get ψ δ ( u ) = ψ ˆ δ ( T u ( X , Y ) ) ⩽ ψ ˆ δ ( 2 T ˆ ( u − 1 ) / 2 ( ( X + Y ) / 2 ) − 1 ) = e δ ψ ˆ 2 δ ( T ˆ ( u − 1 ) / 2 ( ( X + Y ) / 2 ) ) , where ψ ˆ Δ ( T ) = E ( e − Δ T 1 { T < ∞ } ) for a r.v. T and arbitrary Δ > 00$]]>. We also have W u ( n | X , Y ) ⩽ W u ( n | 0 , X + Y ) . Hence, similarly as above we obtain ψ δ ( u ) ⩾ ψ ˆ 2 δ ( T ˆ u / 2 ( ( X + Y ) / 2 ) ) . To summarize, upper and lower bounds for ψ δ ( u ) were obtained: ψ ˆ 2 δ ( T ˆ u / 2 ( ( X + Y ) / 2 ) ) ⩽ ψ δ ( u ) ⩽ e δ ψ ˆ 2 δ ( T ˆ ( u − 1 ) / 2 ( ( X + Y ) / 2 ) ) for all δ > 00$]]> and u ⩾ 1.

 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)

Acknowledgements

We are grateful to the referee for his/her useful comments, and especially for the proposed lower and upper bounds of the Gerber–Shiu function.

The second and the third authors are supported by a grant (No.S-MIP-17-72) from the Research Council of Lithuania.

 References Andrulytė, I.M., Bernackaitė, E., Kievinaitė, D., Šiaulys, J. (2015). A Lundberg-type inequality for an inhomogeneous renewal risk model. Modern Stochastics: Theory and Applications, 2, 173–184. Avram, F., Palmowski, Z., Pistorius, M.R. (2015). On Gerber–Shiu functions and optimal dividend distribution for a Levy risk process in the presence of a penalty function. The Annals of Applied Probability, 25, 1868–1935. Bao, Z., Liu, Y. (2016). A discrete-time risk model with dependence between interclaim arrivals and claim sizes. Advances in Difference Equations, 2016, 188. Bieliauskienė, E., Šiaulys, J. (2012). Gerber–Shiu function for the discrete inhomogeneous claim case. International Journal of Computer Mathematics, 89, 1617–1630. Bischoff-Ferrari, H.A., Orav, J.E., Barrett, J.A., Baron, J.A. (2007). Effect of seasonality and weather on fracture risk in individuals 65 years and older. Osteoporosis International, 18(9), 1225–1233. Bratiichuk, M. (2012). On the Gerber–Shiu function for a risk model with multi-layer dividend strategy. Statistics & Probability Letters, 82, 496–504. Castañer, A., Claramunt, M.M., Gathy, M., Lefèvre, Cl., Mármol, M. (2013). Ruin problems for a discrete time risk model with non-homogeneous conditions. Scandinavian Actuarial Journal, 2013, 83–102. Cheng, S., Gerber, H.U., Shiu, E.S.W. (2000). Discounted probabilities and ruin theory in the compound binomial model. Insurance: Mathematics and Economics, 26, 239–250. Cheung, E.C.K., Liu, H. (2016). On the joint analysis of the total discounted payments to policyholders and shareholders: threshold dividend strategy. Annals of Actuarial Science, 10, 236–269. Cheung, E.C.K., Landriault, D., Badescu, A.L. (2011). On a generalization of the risk model with Markovian claim arrivals. Stochastic Models, 27, 407–430. Cheung, E.C.K., Liu, H., Woo, J-K. (2015). On the joint analysis of the total discounted payments to policyholders and shareholders: dividend barrier strategy. Risks, 3, 491–514. Chi, Y., Jaimungal, S., Lin, X.S. (2010). An insurance risk model with stochastic volatility. Insurance: Mathematics and Economics, 46, 52–66. Cossette, H., Marceau, E., Marri, F. (2011). Constant dividend barrier in a risk model with generalized Farlie-Gumbel-Morgenstern copula. Methodology and Computing in Applied Probability, 13, 487–510. Damarackas, J., Šiaulys, J. (2014). Bi-seasonal discrete time risk model. Applied Mathematics and Computation, 247, 930–940. De La Sen, M. (2016). About some controllability properties of linear discrete-time systems in probabilistic metric spaces. Informatica, 27, 503–526. De Vylder, F.E., Goovaerts, M.J. (1988). Recursive calculation of finite-time ruin probabilities. Insurance: Mathematics and Economics, 7, 1–8. Deng, X., Barnett, B.J., Vedenov, D.V., West, J.W. (2007). Hedging dairy production losses using weather based index insurance. Agricultural Economics, 36, 271–280. Dickson, D.C.M., Waters, H.R. (1991). Recursive calculation of survival probabilities. ASTIN Bulletin, 21, 199–221. Dickson, D.C.M., Qazvini, M. (2016). Gerber–Shiu analysis of a risk model with capital injections. The European Actuarial Journal, 6, 409–440. Dong, Y., Wang, G., Yuen, K.C. (2009). On the renewal risk model under a threshold strategy. Journal of Computational and Applied Mathematics, 230, 22–33. Ferreira, F., Pacheco, A., Ribeiro, H. (2017). Moments of losses during busy-periods of regular and nonpreemptive oscillating M-X/G/1/n systems. Annals of Operations Research, 252, 191–211. Gerber, H.U., Shiu, E.S.W. (1998). On the time value of ruin. North American Actuarial Journal, 2, 48–78. Landriault, D., Willmot, G.E. (2008). On the Gerber–Shiu discounted penalty function in the Sparre Andersen model with an arbitrary inter-claim time distribution. Insurance: Mathematics and Economics, 42, 600–608. Li, S. (2005). On a class of discrete time renewal risk models. Scandinavian Actuarial Journal, 2005, 241–260. Li, S., Garrido, J. (2002). On the time value of ruin in the discrete time risk model. Working Paper 02-18, Business Economics, University Carlos III of Madrid. Li, S., Garrido, J. (2004). On ruin for the Erlang(n) risk process. Insurance: Mathematics and Economics, 34, 391–408. Li, S., Lu, Y. (2013). On the generalized Gerber–Shiu function for surplus processes with interest. Insurance: Mathematics and Economics, 52, 127–134. Li, J.-Z., Wu, R. (2015). The Gerber–Shiu discounted penalty function for a compound binomial risk model with by-claims. Acta Mathematicae Applicatae Sinica, English Series, 31, 181–190. Li, Z., Sendova, K.P. (2015). On a ruin model with both interclaim times and premium depending on claim sizes. Scandinavian Actuarial Journal, 2015, 245–265. Li, S., Lu, Y., Garrido, J. (2009). A review of discrete-time risk models. Revista de la Real Academia de Ciencias Naturales. Serie A, Matemáticas, 103, 321–337. Lin, X.S., Pavlova, K.P. (2006). The compound Poisson risk model with a threshold dividend strategy. Insurance: Mathematics and Economics, 42, 617–627. Lin, X.S., Sendova, K.P. (2008). The compound Poisson risk model with multiple thresholds. Insurance: Mathematics and Economics, 38, 57–80. Lin, X.S., Willmot, G.E., Drekic, S. (2003). The classical risk model with a constant dividend barrier: analysis of the Gerber–Shiu discounted penalty function. Insurance: Mathematics and Economics, 33, 551–566. Liu, S.X., Guo, J.Y. (2006). Discrete risk model revisited. Methodology and Computing in Applied Probability, 8, 303–313. Liu, X., Xiong, J., Zhang, S. (2015). The Gerber–Shiu discounted penalty function in the classical risk model with impulsive dividend policy. Statistics & Probability Letters, 107, 183–190. Liu, C., Zhang, Z., Yang, H. (2017). A note on a discrete time MAP risk model. Journal of Computational and Applied Mathematics, 309, 111–121. Marceau, E. (2009). On the discrete-time compound renewal risk model with dependence. Insurance: Mathematics and Economics, 44, 245–259. Marciniak, E., Palmowski, Z. (2016). On the optimal dividend problem for insurance risks models with surplus-dependent premiums. Journal of Optimization Theory and Applications, 168, 723–742. Mihálýko, É.O., Mihálýko, C. (2011). Mathematical investigation of the Gerber–Shiu function in the case dependent inter-claim time and claim size. Insurance: Mathematics and Economics, 48, 378–383. Pavlova, K.P., Willmot, G.E. (2004). The discrete stationary renewal risk model and the Gerber–Shiu discounted penalty function. Insurance: Mathematics and Economics, 35, 267–277. Răducan, A.M., Vernic, R., Zbăganu, G. (2015a). Recursive calculation of ruin probabilities at or before claim instants for non-identically distributed claims. ASTIN Bulletin, 45, 421–443. Răducan, A.M., Vernic, R., Zbăganu, G. (2015b). On the ruin probability for nonhomogeneous claims and arbitrary inter-claim revenues. Journal of Computational and Applied Mathematics, 290, 319–333. Schmidli, H. (1999). On the distribution of the surplus prior and at ruin. ASTIN Bulletin, 29, 227–244. Schmidli, H. (2015). Extended Gerber–Shiu functions in a risk model with interest. Insurance: Mathematics and Economics, 61, 271–275. Shi, Y., Liu, P., Zhang, C. (2013). On the compound Poisson risk model with dependence and a threshold dividend strategy. Statistics & Probability Letters, 83, 1998–2006. Shiraishi, H. (2016). Review of statistical actuarial risk modelling. Cogent Mathematics & Statistics, 3, 1123945. Tsai, C.C.L. (2003). On the expectations of the present values of the time of ruin perturbed by diffusion. Insurance: Mathematics and Economics, 32, 413–429. Tsai, C.C.L., Willmot, G.E. (2002). A generalized defective renewal equation for the surplus process perturbed by diffusion. Insurance: Mathematics and Economics, 30, 51–66. Willmot, G.E., Dickson, D.C.M. (2003). The Gerber–Shiu discounted penalty function in the stationary renewal risk model. Insurance: Mathematics and Economics, 32, 403–411. Woo, J.-K., Xu, R., Yang, H. (2017). Gerber–Shiu analysis with two-sided acceptable levels. Journal of Computational and Applied Mathematics, 321, 185–210. Xu, L., Yang, H., Wang, R. (2014). Cox risk model with variable premium rate and stochastic return on investment. Journal of Computational and Applied Mathematics, 256, 52–64. Zhang, Z., Cheung, E.C.K. (2016). The Markov additive risk process under an Erlangized dividend barrier strategy. Methodology and Computing in Applied Probability, 18, 275–306. Zhang, Z., Yang, H., Yang, H. (2012). On a Sparre Andersen risk model with time-dependent claim sizes and jump-diffusion perturbation. Methodology and Computing in Applied Probability, 14, 973–995. Zhang, Z., Wu, X., Yang, H. (2014). On a perturbed Sparre Andersen risk model with dividend barrier and dependence. Journal of the Korean Statistical Society, 43, 585–598. Zhang, Z., Cheung, E.C.K., Yang, H. (2017a). Levy insurance risks process with Poissonian taxation. Scandinavian Actuarial Journal, 2017, 51–87. Zhang, Z., Yang, Y., Liu, C. (2017b). On a perturbed compound Poisson model with varying premium rates. Journal of Industrial and Management Optimization, 13, 721–736. Zhou, Z., Xiao, H., Deng, Y. (2015). Markov-dependent risk model with multi-layer dividend strategy. Applied Mathematics and Computation, 252, 273–286.