Numerical Method for a Markov-Modulated Risk Model with Two-Sided Jumps

Abstract

This paper considers a perturbed Markov-modulated risk model with two-sided jumps, where both the upward and downward jumps follow arbitrary distribution. We first derive a system of differential equations for the Gerber-Shiu function. Furthermore, a numerical result is given based on Chebyshev polynomial approximation. Finally, an example is provided to illustrate the method.

1. Introduction

The risk model with two-sided jumps was first proposed by Boucherie et al. [1] and has been further investigated by many authors during the last few years. For example, Kou and Wang [2] studied the Laplace transform of the first passage time and the overshoot for a perturbed compound Poisson model with double exponential jumps. Xing et al. [3] extended the results of Kou and Wang [2] to the case that the surplus process with phase-type downward and arbitrary upward jumps. Zhang et al. [4] assumed that the downward jumps follow arbitrary distribution and the upward jumps have a rational Laplace transform. They derived the Laplace transform of the Gerber-Shiu function by using the roots of the generalized Lundberg equation. Under the assumption that the upward jumps follow Laplace distribution and arbitrary downward jumps, Chi [5] obtained a closed-form expression for the Gerber-Shiu function by applying Wiener-Hopf factorization technique. The applications of the model in finance were also discussed. Jacobsen [6] studied a perturbed renewal risk model with phase-type interclaim times and two-sided jumps, where both the jumps have rational Laplace transforms. Based on the roots of the Cramér-Lundberg equation, the joint Laplace transform on the time to ruin, and the undershoot at ruin were given. However, in all the aforemetioned papers, the topic that the jumps in both directions are arbitrary distributions is still not discussed. The Markov-modulated risk model (Markovian regime switching model) was first proposed by Asmussen [7] to extend the classical risk model. Since then, it has received remarkable attention in actuarial mathematics, see, for example, Zhu and Yang [8, 9], Zhang et al. [4], Ng and Yang [10], Li and Lu [11], Lu and Tsai [12], and references therein. Motivated by the papers mentioned above, in this paper, we will study the Markov-modulated risk model with two-sided jumps.

Let {𝐽(𝑡),𝑡≥0} be a homogenous, irreducible, and recurrent Markov process with finite state space 𝐸={1,2,…,𝑛}. Denote the intensity matrix of {𝐽(𝑡),𝑡≥0} by 𝐀=(𝛼_𝑖𝑗)^𝑛_𝑖,𝑗=1 with ∑^𝑛_𝑗=1𝛼_𝑖𝑗=0 and 𝛼_𝑖𝑖∶=−𝛼_𝑖=−∑^𝑛_𝑗≠𝑖𝛼_𝑖𝑗 for 𝑖∈𝐸. Let {𝑋_𝑖,𝑖=1,2,…} be a sequence of independent random variables representing the jumps, and 𝐵(𝑡) be a standard Brownian motion with 𝐵(0)=0. Here we assume that the premium rates, claim interarrival times, the distributions of the jumps, and the diffusion parameter are all influenced by the environment process {𝐽(𝑡),𝑡≥0}. When 𝐽(𝑡)=𝑖, the premium rate is 𝑐_𝑖, jumps arrive according to a Poisson process with intensity 𝜆_𝑖, the diffusion parameter is 𝜎_𝑖>0, and the size of the jumps which arrives at time 𝑡 follows the distribution 𝐹_𝑖 with density 𝑓_𝑖 and finite mean 𝜇_𝑖. Then the Markov-modulated diffusion risk model {𝑈(𝑡),𝑡≥0} is defined by 𝑈(𝑡)=𝑢+^𝑡₀𝑐_𝐽(𝑠)d𝑠−^𝑁(𝑡)_𝑗=1𝑋_𝑗+^𝑡₀𝜎_𝐽(𝑠)d𝐵(𝑠), (1.1) where 𝑢≥0 is the initial surplus. If we denote the stationary distribution of {𝐽(𝑡),𝑡≥0} by 𝜋=(𝜋₁,𝜋₂,…,𝜋_𝑛), then the positive security loading condition is given by ^𝑛_𝑖=1𝜋_𝑖𝑐_𝑖−𝜆_𝑖𝜇_𝑖>0. (1.2)

In this paper, we further assume that the jumps in (1.1) are two-sided. The upward jumps can be explained as the random income (premium or investment), while the downward jumps are interpreted as the random loss. In this case, the density function is given by 𝑓_𝑖(𝑥)=𝑝_𝑖𝑓_𝑖,𝑑(𝑥)𝐼(𝑥≥0)+𝑞_𝑖𝑓_𝑖,𝑢(−𝑥)𝐼(𝑥<0),for𝐽(𝑡)=𝑖,𝑖=1,2,…,𝑛, (1.3) where 0<𝑝_𝑖≤1, 𝑝_𝑖+𝑞_𝑖=1, 𝐼(⋅) is the indicator function, 𝑓_𝑖,𝑑 and 𝑓_𝑖,𝑢 are two arbitrary functions on [0,∞).

Let 𝑇=inf{𝑡≥0∶𝑈(𝑡)≤0} (∞ otherwise) be the time to ruin. For 𝛿≥0, let 𝜙_𝑖(𝑢)=𝐸𝑒^−𝛿𝑇𝜔(𝑈(𝑇−),|𝑈(T)|)𝐼(𝑇<∞)∣𝐽(0)=𝑖,𝑈(0)=𝑢,𝑢≥0, (1.4) be the Gerber-Shiu function at ruin given that the initial state is 𝑖, where 𝜔(𝑥₁,𝑥₂) is a nonnegative penalty function, 𝑈(𝑇−) is the surplus immediately prior to ruin, and |𝑈(𝑇)| is the deficit at ruin. Without loss of generality, we assume that 𝜔(0,0)=1. Thus 𝜙_𝑖(0)=1 for 𝑖=1,2,…,𝑛. When 𝜔=1, (1.4) reduces to the Laplace transform of the time to ruin 𝜓_𝛿,𝑖(𝑢)=𝐸𝑒^−𝛿𝑇𝐼(𝑇<∞)∣𝐽(0)=𝑖,𝑈(0)=𝑢,𝑢≥0, (1.5) when 𝜔=1 and 𝛿=0, (1.4) reduces to the probability of ruin 𝜓_𝑖(𝑢)=𝑃(𝑇<∞∣𝐽(0)=𝑖,𝑈(0)=𝑢),𝑢≥0. (1.6)

The purpose of this paper is to present some numerical results on the Gerber-Shiu function for the Markov-modulated diffusion risk model with arbitrary upward and downward jumps. In Section 2 we derive a system of integrodifferential equations and approximate solutions for 𝜙_𝑖(𝑢). Numerical example is given in the last section.

2. Integrodifferential Equations and Approximate Solution

Theorem 2.1. For 𝑢≥0, 𝜙_𝑖(𝑢)(𝑖=1,2,…,𝑛) satisfies the following integrodifferential equation 𝜎²_𝑖2𝜙_𝑖^′′(𝑢)+𝑐_𝑖𝜙^′_𝑖(𝑢)−𝜆_𝑖+𝛿𝜙_𝑖(𝑢)+^𝑛_𝑘=1𝛼_𝑖𝑘𝜙_𝑘(𝑢)=−𝜆_𝑖^∞_−∞𝜙_𝑖𝑢−𝑦𝑓𝑦d𝑦=−𝜆_𝑖𝑝_𝑖^𝑢₀𝜙_𝑖𝑢−𝑦𝑓_𝑖,𝑑𝑦d𝑦−𝜆_𝑖𝑞_𝑖^∞₀𝜙_𝑖𝑢+𝑦𝑓_𝑖,𝑢𝑦d𝑦−𝜔_𝑖(𝑢), (2.1) where 𝜔_𝑖(𝑢)=𝜆_𝑖𝑝_𝑖^∞_𝑢𝜔𝑢,𝑦−𝑢𝑓_𝑖,𝑑𝑦d𝑦, (2.2) with boundary conditions 𝜙_𝑖(0)𝜙=1,_𝑖(∞)=0. (2.3)

Proof. Similar to Ng and Yang [10].

Remark 2.2. When 𝑛=1, (2.1) is identical to (3.2) in Zhang et al. [4].

Clearly, (2.1) is a system of second order linear integrodifferential equations of Fredholm-Volterra type. As is well known, it is very difficult to find analytical solution of this system. Motivated by Akyüz-Dascioglu [13], we will study an alternative system defined on [0,1] by Chebyshev collocation method. First, we transform the interval [0,∞) to [0,1]. Following Diko and Usábel [14], we set 𝑢=ℎ(𝑥), that is, ℎ∶[0,1]→[0,∞). Furthermore, we assume that ℎ is an arbitrary strictly monotone, twice continuously differentiable function throughout the paper.

Theorem 2.3. Let ℎ(𝑥) be a monotone increase function and 𝜒_𝑖(𝑥)=𝜙_𝑖(ℎ(𝑥)) for 𝑥∈[0,1]. Then 𝜒_𝑖(𝑥) satisfies the following integrodifferential equation 𝜎²_𝑖2ℎ^′(𝑥)²𝜒_𝑖^′′(𝑥)+𝑐_𝑖ℎ^′(𝑥)−𝜎²_𝑖ℎ^′′(𝑥)2ℎ^′(𝑥)³𝜒^′_𝑖(𝑥)−𝜆_𝑖+𝛿𝜒_𝑖(𝑥)+^𝑛_𝑘=1𝛼_𝑖𝑘𝜒_𝑘(𝑥)+^𝑥₀𝐾_𝑖(𝑥,𝑡)𝜒_𝑖(𝑡)d𝑡+¹₀𝐿_𝑖(𝑥,𝑡)𝜒_𝑖(𝑡)d𝑡+𝑊_𝑖(𝑥)=0, (2.4) where 𝐾_𝑖(𝑥,𝑡)=𝜆_𝑖𝑝_𝑖𝑓_𝑖,𝑑(ℎ(𝑥)−ℎ(𝑡))−𝑞_𝑖𝑓_𝑖,𝑢(ℎ(𝑡)−ℎ(𝑥))ℎ^′(𝑡),𝐿_𝑖(𝑥,𝑡)=𝜆_𝑖𝑞_𝑖𝑓_𝑖,𝑢(ℎ(𝑡)−ℎ(𝑥))ℎ^′(𝑡),𝑊_𝑖(𝑥)=𝜔_𝑖(ℎ(𝑥)), (2.5) with boundary conditions 𝜒_𝑖(0)𝜒=1,_𝑖(1)=0. (2.6)

Proof. By the definitions of function ℎ and 𝜒_𝑖, we have 𝑝_𝑖^𝑢₀𝜙_𝑖𝑢−𝑦𝑓_𝑖,𝑑𝑦d𝑦+𝑞_𝑖^∞₀𝜙_𝑖𝑢+𝑦𝑓_𝑖,𝑢𝑦d𝑦=𝑞_𝑖¹₀𝜙_𝑖(ℎ(𝑡))𝑓_𝑖,𝑢(ℎ(𝑡)−ℎ(𝑥))ℎ^′(𝑡)+d𝑡_𝑥0𝜙_𝑖(ℎ(𝑡))𝑝_𝑖𝑓_𝑖,𝑑(ℎ(𝑥)−ℎ(𝑡))−𝑞_𝑖𝑓_𝑖,𝑢(ℎ(𝑡)−ℎ(𝑥))ℎ^′(𝑡)d𝑡=𝑞_𝑖¹₀𝜒_𝑖(𝑡)𝑓_𝑖,𝑢(ℎ(𝑡)−ℎ(𝑥))ℎ^′(𝑡)+d𝑡_𝑥0𝜒_𝑖(𝑡)𝑝_𝑖𝑓_𝑖,𝑑(ℎ(𝑥)−ℎ(𝑡))−𝑞_𝑖𝑓_𝑖,𝑢(ℎ(𝑡)−ℎ(𝑥))ℎ^′(𝑡)d𝑡. (2.7)
Substituting (2.7) and 𝜒_𝑖(𝑥)=𝜙_𝑖(ℎ(𝑥)) into (2.1) and simplifying lead to (2.4). The boundary conditions are direct result of the boundary conditions in Theorem 2.1. This completes the proof.

Remark 2.4. The existence of the solution for the system of integrodifferential equations (2.4) can be found in Fariborzi and Behzadi [15].

According to Akyüz-Dascioglu [13], 𝜒_𝑖(𝑥) and its derivatives have truncated Chebyshev series expression 𝜒_𝑖^(𝑗)(𝑥)=^𝑁_𝑟=0𝑎^(𝑗)_𝑖𝑟𝑇^∗_𝑟(𝑥),𝑖=1,2,…,𝑛,𝑗=0,1,2,…,𝑥∈[0,1], (2.8) where 𝜒_𝑖⁽⁰⁾(𝑥)=𝜒_𝑖(𝑥), 𝑎⁽⁰⁾_𝑖𝑟=𝑎_𝑖𝑟,𝑇^∗_𝑟(𝑥) are shifted Chebyshev polynomials of the first kind and 𝑎^(𝑗)_𝑖𝑟 are the unknown coefficients to be determined.

Let 𝐓^∗(𝑥)=(𝑇^∗₁(𝑥),𝑇^∗₂(𝑥),…,𝑇^∗_𝑁(𝑥))^⊤, 𝐀_𝑖=(𝑎_𝑖0,𝑎_𝑖1,…,𝑎_𝑖𝑁)^⊤, 𝐀_𝑖^(𝑗)=(𝑎^(𝑗)_𝑖0,𝑎^(𝑗)_𝑖1,…,𝑎^(𝑗)_𝑖𝑁)^⊤. Then (2.8) can be written in the matrix form 𝜒_𝑖(𝑥)=𝐓^∗(𝑥)𝐀_𝑖,𝜒_𝑖^(𝑗)(𝑥)=𝐓^∗(𝑥)𝐀_𝑖^(𝑗)=4^𝑗𝐓^∗_𝑟(𝑥)𝐌^𝑛𝐀_𝑖^(𝑗),𝑗=1,2,…, (2.9) where 01𝐌=2032052⋯𝑁2002040⋯0000305⋯𝑁⋮⋮⋮⋮⋮⋮⋱⋮000000⋯𝑁000000⋯0_{(𝑁+1)×(𝑁+1)}, (2.10) for odd 𝑁, and 01𝐌=2032052⋯0002040⋯𝑁000305⋯0⋮⋮⋮⋮⋮⋮⋱⋮000000⋯0000000⋯𝑁_{(𝑁+1)×(𝑁+1)}, (2.11) for even 𝑁.

Similarly, the kernel functions 𝐾_𝑖(𝑥,𝑡) and 𝐿_𝑖(𝑥,𝑡) can be expanded to univariate Chebyshev series 𝐾_𝑖(𝑥,𝑡)=𝜿_𝑖(𝑥)𝑇^∗_𝑟(𝑡),𝐿_𝑖(𝑥,𝑡)=𝓵_𝑖(𝑥)𝑇^∗_𝑟(𝑡), (2.12) where 𝜿_𝑖(𝑥)=12𝜅_𝑖,1(𝑥),𝜅_𝑖,2(𝑥)1,…,2𝜅_𝑖,𝑛(𝑥,𝓵_𝑖)(𝑥)=12ℓ_𝑖,1(𝑥),ℓ_𝑖,2(𝑥)1,…,2ℓ_𝑖,𝑛(𝑥), (2.13) with 𝜅_𝑖,𝑟 and ℓ_𝑖,𝑟 are Chebyshev coefficients determined by Clenshaw and Curtis [16].

Theorem 2.5. For 0≤𝑥≤1, an approximate expression for 𝜒_𝑖(𝑥) is given by 𝜒_𝑖(𝑥)=𝐓^∗(𝑥)𝐀_𝑖,𝑖=1,2,…,𝑛, (2.14) where the column vector 𝐀_𝑖 can be determined by the following systems 2𝜎²_𝑖ℎ^′𝑥_𝑗²𝐓^∗𝑥_𝑗𝐌𝐀_𝑖+𝑐_𝑖ℎ^′𝑥_𝑗−𝜎²_𝑖ℎ^′′𝑥_𝑗2ℎ^′𝑥_𝑗³𝐓^∗𝑥_𝑗𝐀_𝑖−𝜆_𝑖+𝛿𝐓^∗𝑥_𝑗𝐀_𝑖+^𝑛_𝑘=1𝛼_𝑖𝑘𝐓^∗𝑥_𝑗𝐀_𝑘+12𝜿_𝑖𝑥_𝑗𝐙2𝑥_𝑗−1𝐀_𝑖+12𝓵_𝑖𝑥_𝑗𝐙𝐀_𝑖+𝑊_𝑖𝑥_𝑗𝐓=0,_∗(0)𝐀_𝑖𝐓=1,_∗(1)𝐀_𝑖=0, (2.15) where matrix 𝐙=(𝑧_𝑖𝑗) with elements 𝑧_𝑖𝑗=11−𝑖−𝑗²+11−𝑖+𝑗²,𝑓𝑜𝑟𝑒𝑣𝑒𝑛𝑖+𝑗,0,𝑓𝑜𝑟𝑜𝑑𝑑𝑖+𝑗, (2.16) matrix 𝐙(𝑥)=(𝑧_𝑖𝑗(𝑥)) with elements 𝑧_𝑖𝑗(𝑥)=142𝑥²𝑇−2,𝑖+𝑗=1,_𝑖+𝑗+1−𝑇𝑖+𝑗+1^{𝑖+𝑗−1}−1𝑖+𝑗−1+1𝑖+𝑗+1𝑖+𝑗−1+𝑥²𝑇−1,𝑖−𝑗=1,_𝑖+𝑗+1+𝑇𝑖+𝑗+1^{1−𝑖−𝑗}+𝑇1−𝑖−𝑗^{1+𝑖−𝑗}+𝑇1+𝑖−𝑗^{1−𝑖+𝑗}11−𝑖+𝑗+21−𝑖+𝑗²+11−𝑖−𝑗²𝑇,𝑓𝑜𝑟𝑒𝑣𝑒𝑛𝑖+𝑗,_𝑖+𝑗+1+𝑇𝑖+𝑗+1^{1−𝑖−𝑗}+𝑇1−𝑖−𝑗^{1+𝑖−𝑗}+𝑇1+𝑖−𝑗^{1−𝑖+𝑗}11−𝑖+𝑗−21−𝑖+𝑗²+11−𝑖−𝑗²,𝑓𝑜𝑟𝑜𝑑𝑑𝑖+𝑗 (2.17) and 𝑥_𝑗(𝑗=0,1,2,…,𝑁−1) are collocations.

Proof. Using (2.8) and (2.12), one obtains ^𝑥₀𝐾_𝑖(𝑥,𝑡)𝜒_𝑖(𝑡)d𝑡=^𝑥₀𝜿_𝑖(𝑥)𝐓^∗(𝑡)^⊤𝐓^∗(𝑡)𝐀_𝑖=1d𝑡2𝜿_𝑖(𝑥)𝐙(2𝑥−1)𝐀_𝑖,₁₀𝐿_𝑖(𝑥,𝑡)𝜒_𝑖(𝑡)d𝑡=¹₀𝓵_𝑖(𝑥)𝐓^∗(𝑡)^⊤𝐓^∗(𝑡)𝐀_𝑖=1d𝑡2𝓵_𝑖(𝑥)𝐙𝐀_𝑖, (2.18) Substituting (2.18) into (2.4), we have 2𝜎²_𝑖ℎ^′(𝑥)²𝐓^∗(𝑥)𝐌𝐀_𝑖+𝑐_𝑖ℎ^′(𝑥)−𝜎²_𝑖ℎ^′′(𝑥)2ℎ^′(𝑥)³𝐓^∗(𝑥)𝐀_𝑖−𝜆_𝑖+𝛿𝐓^∗(𝑥)𝐀_𝑖+^𝑛_𝑘=1𝛼_𝑖𝑘𝐓^∗(𝑥)𝐀_𝑘+12𝜿_𝑖(𝑥)𝐙(2𝑥−1)𝐀_𝑖+12𝓵_𝑖(𝑥)𝐙𝐀_𝑖+𝑊_𝑖(𝑥)=0, (2.19) which is identical to (2.15) in form. Substituting the collocations 𝑥_𝑗(𝑗=0,1,…,𝑁−1) into (2.19) leads to (2.15). 𝐓^∗(0)𝐀_𝑖=1 and 𝐓^∗(1)𝐀_𝑖=0 can be obtained by (2.6).

Example 2.6. To illustration our method, we use the example of Zhang et al. [4]. Let 𝑛=1, 𝑐=2, 𝜎²=2, 𝜆=1, 𝑝₁=0.6, 𝑞₁=0.4, 𝛿=0.3, the downward jumps are exponentially distributed with parameter 0.3, and the upward jump density is given by 𝑓_1,𝑢(𝑥)=0.08𝑒^−0.4𝑥+0.64𝑒^−0.8𝑥. We set 𝑢=ℎ(𝑥)=−ln(1−𝑥) and the collocation points are 𝑥_𝑗=((1+cos(𝑠𝜋/𝑁))/2)(𝑖=1,2,…,𝑁).

Figure 1 shows that the approximate solution is very near to the exact solution for any initial surplus 𝑢. We remark that the horizontal axis in Figure 1 is cos(𝑠𝜋/𝑁)(𝑖=1,2,…,𝑁) and 𝑢=−ln((1−cos(𝑠𝜋/𝑁))/2).

Figure 1 

Laplace transform of the time to ruin.

From Table 1 we can see that the errors between the approximate solutions and the exact solutions decrease when 𝑁 increases. The initial surplus 𝑢 can also influence the approximate solution: the bigger 𝑢 need a bigger 𝑁 to decrease the error.

Table 1 

Laplace transform of the time to ruin.

Acknowledgments

This work is supported by the Natural Science Foundation of Shandong (no. ZR2010AQ015), the Tianyuan fund for Mathematics (no. 11226251), and the Natural Science Foundation of Qufu Normal University (no. 2012ZRB01473).