The Perturbed Dual Risk Model with Constant Interest and a Threshold Dividend Strategy

Abstract

We consider the perturbed dual risk model with constant interest and a threshold dividend strategy. Firstly, we investigate the moment-generation function of the present value of total dividends until ruin. Integrodifferential equations with certain boundary conditions are derived for the present value of total dividends. Furthermore, using techniques of sinc numerical methods, we obtain the approximation results to the expected present value of total dividends. Finally, numerical examples are presented to show the impact of interest on the expected present value of total dividends and the absolute ruin probability.

1. Introduction

In insurance mathematics, the classical risk model has been the center of focus for decades. The surplus in the classical model at time 𝑡 can be expressed as 𝑈(𝑡)=𝑢+𝑐𝑡−𝑆(𝑡),𝑡≥0, (1) where 𝑢≥0 is the initial surplus, 𝑐>0 is the premium rate, and 𝑆(𝑡) is the aggregate claims by time 𝑡, usually modeled by a compound Poisson process. In recent years, quite a few interesting papers have been written on a model which is dual to the classical insurance risk model. See, for example, Avanzi et al. [1], Avanzi and Gerber [2]. The surplus of a dual classical risk model at time 𝑡 is 𝑈(𝑡)=𝑢−𝑐𝑡+𝑆(𝑡),𝑡≥0. (2)

In this model, the premium rate should be viewed as an expense rate and claims should be viewed as profits or gains. While not very popular in insurance mathematics, this model has appeared in various literature (see Cramér [3, Section  5.13] and Seal [4, pages 116–119] and Takács [5, pages 152–154]). There were many possible interpretations for this model. For example, we can treat the surplus as the amount of the capital of a business engaged in research and development. The company paid expenses for research, and occasional profit of random amounts (such as the award of a patent or a sudden increase in sales) arises according to a Poisson process. A similar model was used by Bayraktar and Egami [6] to model the capital of a venture capital investment. Another model was an annuity business. The company issues payments continuously to annuitants, while the gross reserve of an annuitant was released as emerging profit when he died. Yang and Zhu [7] generalized the dual model into a regime-switching setting and calculated bounds for ruin probabilities. One of the current topics of interest in insurance mathematics is the problem of maximizing the expected total discounted dividends until ruin, which goes back to de Finetti [8], and it had also been studied by Bühlmann [9, Section  6.4], Gerber [10, Sections  7 and 8], and Gerber [11, Section 10.1]. In Avanzi et al. [1], the authors studied the expected total discounted dividends until ruin for the dual model under the barrier strategy by means of integrodifferential equations. They derived explicit formulas when profits or gains followed an exponential or a mixture of exponential distributions and showed that the optimal value of the dividend barrier under the dual model was independent of the initial surplus. Avanzi and Gerber [2] studied a dual model perturbed by diffusion and discussed how the optimal value of the dividend barrier can be determined. Albrecher et al. [12] studied a dual model that also paid taxes when the surplus was at a running maximum and calculated the expected total discounted dividends before ruin for exponentially distributed profits.

Considering the perturbed dual risk model, the surplus of an insurer has the following form 𝑈(𝑡)=𝑢−𝑐𝑡+^𝑁(𝑡)_𝑖=1𝑋_𝑖+𝜎𝐵(𝑡),𝑡≥0, (3) where 𝑢≥0 is the initial surplus and 𝑐>0 is the constant rate. {𝑋_𝑖,𝑖=1,2,…} is a sequence of independent income size random variables with a common distribution function 𝐹(𝑥) with 𝐹(0)=0 and density function 𝑓(𝑥). {𝑁(𝑡),𝑡≥0} is the Poisson income-number process with an intensity 𝜆>0 and is defined as 𝑁(𝑡)=sup{𝑘∶𝑇₁+⋯+𝑇_𝑘≤𝑡}, where the i.i.d. interincome times {𝑇_𝑖}^∞_𝑖=1 have a common exponential distribution with a parameter 𝜆. {𝐵(𝑡),𝑡≥0} is a standard Brownian motion with 𝐵(0)=0 and 𝜎>0 is a constant, representing the diffusion volatility parameter. In addition, {𝑋_𝑖,𝑖=1,2,…}, {𝑁(𝑡),𝑡≥0} and {𝐵(𝑡),𝑡≥0}, are mutually independent.

Dividend strategy for insurance risk models was first proposed by de Finetti [8] to reflect more realistically the surplus cash flowed in an insurance portfolio. Various threshold strategies have been studied by many authors, including Gerber and Landry [13], Ng [14], and Lin and Pavlova [15]. Among them, Ng [14] showed that the threshold strategy was optimal when the dividend rate was bounded and the individual claim distribution was exponential. Recently, the threshold dividend strategy has been considered in the class of compound Poisson process perturbed by diffusion and its generalizations; readers may refer to Avanzi and Gerber [2], Gao and Liu [16], Wan [17], and the references therein.

Then, we consider the modification of the surplus process by a threshold strategy with a threshold level 𝑏. When 𝑈(𝑡) is below 𝑏, no dividends are paid and the surplus decreases at the original rate 𝑐₁. But when 𝑈(𝑡) is above 𝑏, the surplus would decrease at a different rate 𝑐₂>𝑐₁ and dividends are paid at rate 𝑐₂−𝑐₁. Incorporating the threshold strategy into (3) yields the surplus process {𝑈_𝑏(𝑡),𝑡≥0} which can be expressed by 𝑑𝑈_𝑏(𝑡)=−𝑐₁𝑑𝑡+𝑑𝑆(𝑡)+𝜎𝑑𝐵(𝑡),𝑈_𝑏(𝑡)<𝑏,−𝑐₂𝑑𝑡+𝑑𝑆(𝑡)+𝜎𝑑𝐵(𝑡),𝑈_𝑏(𝑡)≥𝑏, (4) where 𝑈_𝑏(0)=𝑢, 0<𝑏<∞, and 𝑆(𝑡)=∑^𝑁(𝑡)_𝑖=1𝑋_𝑖.

In this paper, we consider that the insurance company earns credit interest with a constant force 𝑟>0. In the mean time, the insurer will receive interest. Incorporating interest into (4) yields the surplus process {𝑈_𝑏(𝑡),𝑡≥0} which can be expressed by 𝑑𝑈_𝑏(𝑡)=−𝑐₁𝑑𝑡+𝑟𝑈_𝑏(𝑡)𝑑𝑡+𝑑𝑆(𝑡)+𝜎𝑑𝐵(𝑡),𝑈_𝑏(𝑡)<𝑏,−𝑐₂𝑑𝑡+𝑟𝑈_𝑏(𝑡)𝑑𝑡+𝑑𝑆(𝑡)+𝜎𝑑𝐵(𝑡),𝑈_𝑏(𝑡)≥𝑏. (5)

Let 𝐷(𝑡) denote the cumulative amount of dividends paid out up to time 𝑡 and 𝛿>0 the force of interest; then 𝐷_𝑢,𝑏=^𝑇𝑏₀𝑒^−𝛿𝑡𝑑𝐷(𝑡) (6) is the present value of all dividends until 𝑇_𝑏, where 𝑇_𝑏 denoted by 𝑇_𝑏=inf{𝑡≥0∶𝑈_𝑏(𝑡)≤0} is the time of ruin. An alternative expression for 𝐷_𝑢,𝑏 is 𝐷_𝑢,𝑏=𝑐₂−𝑐₁^𝑇𝑏₀𝑒^−𝛿𝑡𝐼𝑈_𝑏(𝑡)>𝑏𝑑𝑡, (7) with 𝐼(⋅) denoting the indicator function. It is obvious that 0<𝐷_𝑢,𝑏≤(𝑐₂−𝑐₁)/𝛿.

In the sequel, we will be interested in the following moment generating function: 𝑀𝑢,𝑦;𝑏=𝐸𝑒^{𝑦𝐷𝑢,𝑏} (8) (for those values of 𝑦 where it exist) and the expected total discounted dividends until ruin is 𝑉(𝑢;𝑏)=𝐸𝐷_𝑢,𝑏∣𝑈_𝑏(0)=𝑢. (9) Throughout this paper, we assume that 𝑀(𝑢,𝑦;𝑏) and 𝑉(𝑢;𝑏) are sufficiently smooth functions in 𝑢 and 𝑦, respectively.

2. Integrodifferential Equations

2.1. Integrodifferential Equations for 𝑀(𝑢,𝑦;𝑏)

In this section, we will give the integrodifferential equations satisfied by the moment generating function 𝑀(𝑢,𝑦;𝑏), respectively. It is easy to see that 𝑀(𝑢,𝑦;𝑏) behaves differently with different initial surplus. Hence, for notation convenience, we set 𝑀=𝑀𝑢,𝑦;𝑏¹𝑀𝑢,𝑦;𝑏,if0<𝑢<𝑏,₂𝑢,𝑦;𝑏,if𝑏≤𝑢<∞. (10) In the case of 𝜎=0, we write 𝑀(𝑢,𝑦;𝑏) as 𝑀⁰(𝑢,𝑦;𝑏), and the other one is similar.

In the following, we firstly derive the integrodifferential equations satisfied by 𝑀⁰(𝑢,𝑦;𝑏).

Theorem 1. When 0<𝑢<𝑏, 𝑟𝑢−𝑐₁𝜕𝑀⁰₁𝑢,𝑦;𝑏𝜕𝑢=𝑦𝛿𝜕𝑀⁰₁𝑢,𝑦;𝑏𝜕𝑦+𝜆𝑀⁰₁𝑢,𝑦;𝑏−𝜆₀^𝑏−𝑢𝑀⁰₁𝑢+𝑥,𝑦;𝑏𝑓(𝑥)𝑑𝑥−𝜆^∞_𝑏−𝑢𝑀⁰₂𝑢+𝑥,𝑦;𝑏𝑓(𝑥)𝑑𝑥, (11) and  for 𝑏≤𝑢<∞, 𝑟𝑢−𝑐₂𝜕𝑀⁰₂𝑢,𝑦;𝑏𝜕𝑢=𝑦𝛿𝜕𝑀⁰₂𝑢,𝑦;𝑏𝜕𝑦+𝜆𝑀⁰₂𝑢,𝑦;𝑏−𝜆^∞₀𝑀⁰₁𝑢+𝑥,𝑦;𝑏𝑓(𝑥)𝑑𝑥, (12) with the following boundary conditions: 𝑀⁰₁𝑏−,𝑦;𝑏=𝑀⁰₂𝑏,𝑦;𝑏,lim_𝑢→∞𝑀⁰₂𝑢,𝑦;𝑏=𝑒^{((𝑐2−𝑐1)/𝑦)𝛿},𝑟𝑏−𝑐₁𝜕𝑀⁰₁𝑢,𝑦;𝑏𝜕𝑢_𝑢=𝑏−=𝑟𝑏−𝑐₂𝜕𝑀⁰₂𝑢,𝑦;𝑏𝜕𝑢_𝑢=𝑏. (13)(14)(15)

Proof. When 0<𝑢<𝑏, consider 𝑡>0 such that the surplus cannot reach level 𝑏 by time; that is, ℎ_𝛿(𝑢,𝜏)=𝑒^𝑟𝜏(𝑢−(𝑐₁/𝑟))+(𝑐₁/𝑟)>0. By conditioning on the time and amount of first claim, if it occurs by 𝑡, and whether the claim causes ruin, one gets 𝑀⁰₁𝑢,𝑦;𝑏=(1−𝜆𝜏)𝑀⁰₁ℎ_𝛿(𝑢,𝜏),𝑦𝑒^−𝛿𝜏𝑀;𝑏+𝜆𝜏𝐸⁰₁ℎ_𝛿(𝑢,𝜏)+𝑥,𝑦𝑒^−𝛿𝜏;𝑏. (16) By Taylor’s expansion, 𝑀⁰₁ℎ_𝛿(𝑢,𝜏),𝑦𝑒^−𝛿𝜏;𝑏=𝑀⁰₁𝑢,𝑦;𝑏+𝑟𝑢−𝑐₁𝜏𝜕𝑀⁰₁𝑢,𝑦;𝑏𝜕𝑢−𝑦𝛿𝜏𝜕𝑀⁰₁𝑢,𝑦;𝑏(𝜏).𝜕𝑦+𝑜 (17) Also, due to 𝐸𝑀⁰₁ℎ_𝛿(𝑢,𝜏),𝑦𝑒^−𝛿𝜏=;𝑏_∞0𝑀⁰₁ℎ_𝛿(𝑢,𝜏)+𝑥,𝑦𝑒^−𝛿𝜏;𝑏𝑓(𝑥)𝑑𝑥. (18) Substituting the above expansion into (16), dividing both sides of (16) by 𝑡, and letting 𝑡→0, we can get (11). Using the same way as (11), we can easily prove (12) when 𝑢→∞, 𝑇_𝑏=∞ so the condition (12) is correct, by the method used in Gao and Liu [16], and it is easy to check that the boundary conditions hold. This ends the proof of Theorem 1.

Now we consider the case of 𝜎≠0. For 𝑡≥0, define 𝑊(𝑡)=𝜎∫^𝑡₀𝑒^−𝑟𝑠𝑑𝐵(𝑠). Thus, 𝑊(𝑡) is an Ito stochastic integral. Denote by {⟨𝑊⟩(𝑡),𝑡≥0} the variance process of {𝑊(𝑡),𝑡≥0}. We have ⟨𝑊⟩(𝑡)=𝜎²∫^𝑡₀𝑒^−2𝑟𝑠𝑑𝑠=(2𝑟)⁻¹𝜎²(1−𝑒^−2𝑟𝑠) for 𝑡≥0. Let 𝑣(𝑠)=inf{𝑡∶⟨𝑊⟩(𝑡)>𝑠}; then 𝑣(𝑠)=12𝜎ln²𝜎²𝜎−2𝑟𝑠,0<𝑠<².2𝑟 (19)

Set 𝑊(𝑡)=𝐵(𝑣(𝑡)) for 𝑡≥0, 𝑖=1,2. By the time change of Brownian motion, we have that 𝑊_𝑖 is local standard Brownian motion with 𝑊_𝑖(0)=0 running for an amount of time 𝜎²/2𝑟.

Let 𝜙(𝑢,𝑠)=𝑒^𝑟𝑠(𝑢+𝑐₁∫^𝑠₀𝑒^−𝑟𝑡𝑑𝑡) and 𝜙₁(𝑢,𝑠)=𝑒^𝑟𝑠(𝑢+𝑐₂∫^∞₀𝑒^−𝑟𝑡𝑑𝑡). Then, we have the following results.

Theorem 2. When 0<𝑢<𝑏, 12𝜎²𝜕²𝑀₁𝑢,𝑦;𝑏𝜕𝑢²+𝑟𝑢−𝑐₁𝜕𝑀₁𝑢,𝑦;𝑏𝜕𝑢=𝑦𝛿𝜕𝑀₁𝑢,𝑦;𝑏𝜕𝑦+𝜆𝑀₁𝑢,𝑦;𝑏−𝜆₀𝑏−𝑢𝑀₁𝑢+𝑥,𝑦;𝑏𝑓(𝑥)𝑑𝑥+𝜆^∞_𝑏−𝑢𝑀₂𝑢+𝑥,𝑦;𝑏𝑓(𝑥)𝑑𝑥, (20) and when 𝑏≤𝑢<∞, 12𝜎²𝜕²𝑀₂𝑢,𝑦;𝑏𝜕𝑢²+𝑟𝑢−𝑐₂𝜕𝑀₂𝑢,𝑦;𝑏𝜕𝑢=𝑦𝛿𝜕𝑀₂𝑢,𝑦;𝑏+𝑐𝜕𝑦𝜆−₂−𝑐₁𝑀₂𝑢,𝑦;𝑏+𝜆^∞₀𝑀₂𝑢+𝑥,𝑦;𝑏𝑓(𝑥)𝑑𝑥, (21) with the following boundary conditions: 𝑀₁𝑀0,𝑦;𝑏=1,₁𝑏−,𝑦;𝑏=𝑀₂𝑏,𝑦;𝑏,lim_𝑢→∞𝑀₂𝑢,𝑦;𝑏=𝑒^{((𝑐2−𝑐1)/𝑦)𝛿},𝑟𝑏−𝑐₁𝜕𝑀₁𝑢,𝑦;𝑏∣𝜕𝑢^𝑢=𝑏−=𝑟𝑏−𝑐₂𝜕𝑀₂𝑢,𝑦;𝑏∣𝜕𝑢^𝑢=𝑏. (22)(23)(24)(25)

Proof. For 0<𝑢<𝑏, assume that 𝜖,𝑡>0 such that 𝜖<𝑢<𝑏. Define 𝑇₀=𝑇^𝜖_𝑡∧𝑇₁ and 𝑇^𝜖_𝑡=inf{𝑠>0∶𝜙(𝑢+𝑊(𝑠),𝑠)∉(𝜖,𝑏)}∧𝑡. We have 𝑃(𝑇₀<∞)=1 for all 𝑠∈(0,𝑇₀). Then, by the strong Markov property, we get 𝑀₁𝑢,𝑦;𝑏=𝐸^𝑢𝑒^{𝑦𝐷𝑈,𝑏}=𝐸^𝑢𝑀₁𝑈_𝑏𝑇₀,𝑦𝑒^−𝛿𝑇0𝐼𝑇;𝑏=𝐸₁>𝑡𝑀₁𝑈_𝑏𝑇^ℰ_𝑡,𝑦𝑒^{−𝛿𝑇ℰ𝑡}𝐼𝑇;𝑏+𝐸₁≤𝑡𝑀₁𝑈_𝑏𝑇^ℰ_𝑇1,𝑦𝑒^{−𝛿𝑇ℰ𝑇1}(𝑡)(𝑡).;𝑏=𝐈+𝐈𝐈 (26) By the assumption of independence, one gets 𝐈(𝑡)=𝑒^−𝜆𝑡𝐸𝑀₁𝑈_𝑏𝑇^ℰ_𝑡,𝑦𝑒^{−𝛿𝑇ℰ𝑡},(𝑡)=;𝑏𝐈𝐈^𝑡₀𝜆𝑒^−𝜆𝑠𝐸𝐼𝑇^ℰ_𝑡=𝑠𝑀₁𝑈_𝑏(𝑠),𝑦𝑒^{−𝛿𝑇ℰ𝑡};𝑏𝑑𝑠+^𝑡₀𝜆𝑒^−𝜆𝑠𝐸𝐼𝑇^ℰ_𝑡<𝑠𝑀₁𝑈_𝑏(𝑠),𝑦𝑒^{−𝛿𝑇ℰ𝑡}=;𝑏𝑑𝑠_𝑡0𝜆𝑒^−𝜆𝑠𝐸𝐼𝑇^ℰ_𝑡=𝑠𝑔_𝑠𝜙((𝑠)+𝑢+𝑊,𝑠)𝑑𝑠_𝑡0𝜆𝑒^−𝜆𝑠𝐸𝐼𝑇^ℰ_𝑡<𝑠𝑀₁𝑈_𝑏(𝑠),𝑦𝑒^{−𝛿𝑇ℰ𝑡};𝑏𝑑𝑠, (27) where 𝑔_𝑠𝜙(𝑢+𝑊(𝑠),𝑠)=₀^{𝑏−𝜙(𝑢+𝑊(𝑠),𝑠)}𝑀₁𝜙(𝑢+𝑊(𝑠),𝑠)+𝑥,𝑦𝑒^−𝛿𝑠+;𝑏𝑑𝑠_{∞((𝑠))𝑏−𝜙𝑢+𝑊,𝑠}𝑀₂𝜙(𝑢+𝑊(𝑠),𝑠)+𝑥,𝑦𝑒^−𝛿𝑠;𝑏𝑑𝑠. (28) Noting that lim_𝑡→0(𝑇^𝜖_𝑡=𝑡)=1, lim_𝑡→0(𝑇^𝜖_𝑡<𝑡)=0 a.s. and using Ito’s formula, we have lim_𝑡→0𝐈(𝑡)−𝑀₁𝑢,𝑦;𝑏𝑡=12𝜎²𝜕²𝑀⁰₁𝑢,𝑦;𝑏𝜕𝑢²+𝑟𝑢−𝑐₁𝜕𝑀⁰₁𝑢,𝑦;𝑏𝜕𝑢−𝑦𝛿𝜕𝑀⁰₁𝑢,𝑦;𝑏𝜕𝑦−𝜆𝑀⁰₁𝑢,𝑦;𝑏,lim_𝑡→0𝐈𝐈(𝑡)𝑡=𝜆₀^𝑏−𝑢𝑀⁰₁𝑢+𝑥,𝑦;𝑏𝑓(𝑥)𝑑𝑥+𝜆^∞_𝑏−𝑢𝑀⁰₂𝑢+𝑥,𝑦;𝑏𝑓(𝑥)𝑑𝑥. (29) From (26) and (29), we can get (20) for all 𝑢∈(𝜖,). Hence, (20) holds in (0,𝑏).
Similarly, using the above method, we get (21). Conditions (22) and (23) are obvious. By (14) and (15) and by the weak convergence method used in [17], it is easy to check that the boundary conditions (24) and (25) hold. This ends the proof of Theorem 2.

2.2. Integrodifferential Equations for 𝑉(𝑢,𝑏)

Theorem 3. When 0<𝑢<𝑏, 12𝜎²𝑉₁^′′(𝑢)+𝑟𝑢−𝑐₁𝑉^′₁(𝑢)−(𝜆+𝛿)𝑉₁(𝑢)+𝜆₀^𝑏−𝑢𝑉₁(𝑢+𝑥;𝑏)𝑓(𝑥)𝑑𝑥+𝜆^∞_𝑏−𝑢𝑉₂(𝑢+𝑥;𝑏)𝑓(𝑥)𝑑𝑥=0, (30) and when 𝑏≤𝑢<∞, 12𝜎²𝑉₂^′′(𝑢)+𝑟𝑢−𝑐₂𝑉^′₂(𝑢)−(𝜆+𝛿)𝑉₂(𝑢)+𝜆^∞₀𝑉₂(𝑢+𝑥;𝑏)𝑓(𝑥)𝑑𝑥+𝑐₂−𝑐₁=0, (31) with the following boundary conditions: 𝑉₁()𝑉0;𝑏=0,₁(𝑏−;𝑏)=𝑉₂(𝑏;𝑏),lim_𝑢→∞𝑉₂()=𝑐𝑢;𝑏²−𝑐₁𝛿,𝑟𝑏−𝑐₁𝜕𝑉₁(𝑢;𝑏)𝜕𝑢_𝑢=𝑏−=𝑟𝑏−𝑐₂𝜕𝑉₂(𝑢;𝑏)𝜕𝑢_𝑢=𝑏. (32)

3. Numerical Analysis

The second order system of integrodifferential equations such as in Theorem 3 does not have an explicit solution except in some special cases. In this section, we propose the approximate solution of this equation via the use of a collocation based on sinc methods. The sinc method is highly efficient numerical method developed by Frank Stenger. It is widely used in various fields of numerical analysis such as interpolation, quadrature, approximation of transforms, and solution of integral, ordinary differential, and partial differential equations. An introduction to the sinc approximation theory can be found in the Appendix.

3.1. Numerical Solution of the Expected Present Value of Total Dividends 𝑉(𝑢;𝑏)

To construct an approximation on the interval (0,∞), we consider the following conformal map: 𝜙(𝑧)=log𝑧. (33)

The function 𝜙 also provides a one-to-one transformation of (0,∞) onto the real line ℝ. The sinc grid points 𝑧_𝑘 are defined for ℎ>0 and 𝑘=0, ±1,±2,… by 𝑧_𝑘=𝜙⁻¹(𝑘ℎ)=𝑒^𝑘ℎ. (34)

In order to adopt the sinc method procedure, we arrange the systems in Theorem 3 into the following integrodifferential equation: 12𝜎²𝑉^′′(𝑢)+𝑟𝑢−𝑐₂𝐼(𝑢≥𝑏)+𝑟𝑢−𝑐₁𝐼(0<𝑢<𝑏)×𝑉^′(𝑢)−(𝜆+𝛿)𝑉(𝑢)+𝜆^∞₀𝑉()(𝑥)𝑐𝑢+𝑥;𝑏×𝑓𝑑𝑥+₂−𝑐₁𝐼(𝑢>𝑏)=0. (35)

Using the properties of convolution, the above equation can be further written as 12𝜎²𝑉^′′(𝑢)+𝑟𝑢−𝑐₂𝐼()+𝑢≥𝑏𝑟𝑢−𝑐₁𝐼(0<𝑢<𝑏)𝑉^′(𝑢)−(𝜆+𝛿)𝑉(𝑢)+𝜆^∞_𝑢𝑉()𝑓()+𝑐𝑥;𝑏𝑥−𝑢𝑑𝑥₂−𝑐₁𝐼(𝑢>𝑏)=0, (36) with the following boundary conditions: 𝑉(0;𝑏)=0,lim_𝑢→∞𝑉()=𝑐𝑢;𝑏²−𝑐₁𝛿. (37)

Set 𝑊(𝑢)=𝑉(𝑢;𝑏)−𝐿𝑉(𝑢;𝑏)=𝑉(𝑢;𝑏)−(𝑢/(1+𝑢))((𝑐₂−𝑐₁)/𝛿); then 𝛼𝑊(𝑢)∈𝐿_,𝛽_,𝑑(𝜙) and satisfies 𝑊^′′(𝑢)+𝜇(𝑢)𝑊^′(𝑢)−2𝜎²()𝑊(𝑢)+2𝜆+𝛿𝜎²𝜆^∞_𝑢𝑊(𝑥)𝑓()(𝑢)𝑢−𝑥𝑑𝑥+𝑅=0, (38) with the following boundary conditions: 𝑊(0)=0,lim_𝑢→∞𝑊(𝑢)=0, (39) where 𝜇(𝑢)=2𝜎²𝑟𝑢−𝑐₂𝐼()+𝑢≥𝑏𝑟𝑢−𝑐₁𝐼(,𝑅(𝑢)=20<𝑢<𝑏)𝜎²𝑐₂−𝑐₁𝐼()−2𝑢≥𝑏(1+𝑢)³𝑐₂−𝑐₁𝛿(𝑢)1+𝜇(1+𝑢)²𝑐₂−𝑐₁𝛿−2𝜎²()𝑢𝜆+𝛿𝑐1+𝑢²−𝑐₁𝛿2+𝜆𝜎²𝑐₂−𝑐₁𝛿^∞_𝑢𝑥𝑓()1+𝑥𝑥−𝑢𝑑𝑥. (40)

Let 𝜙 be defined by (33). The sinc grid points are defined for ℎ>0 by 𝑢_𝑘=𝑒^𝑘ℎ,𝑘=0,±1,±2,…. (41)

Then, by using Theorems A.3 and A.6 in the Appendix, we have ^∞_𝑢𝑓(𝑥−𝑢)𝑊(𝑥)𝑑𝑥≈^𝑁𝑁_{𝑙=−𝑀𝑖=−𝑀}𝑤^∗_𝑖𝐴_𝑖𝑙𝑊_𝑙,𝑊(𝑥)≈𝑊(𝑥)=^𝑁_𝑙=−𝑀𝑤_𝑙𝑊_𝑙, (42)(43) where 𝐴=𝑌𝐹(𝑆)𝑌⁻¹, 𝐴_𝑖𝑙 is the (𝑖,𝑙)th element of matrix 𝐴. 𝑊_𝑙 denotes an approximate value of 𝑊(𝑢_𝑙), and 𝜙(𝑥)=log𝑥, 𝜙(0)=−∞, 𝜙(+∞)=+∞, and 𝜙^′(𝑥)=1/𝑥. 𝑤_𝑗 and 𝑤^∗_𝑗 for 𝑗=−𝑀,…,𝑁 are defined as Definition A.2 in the Appendix.

Since 𝛼𝑊(𝑥)∈𝐿_,𝛽_,𝑑(𝜙), by Theorem A.3 in the Appendix, then it is convenient to take 𝑤_𝑗=𝛾_𝑗=𝑆(𝑗,ℎ)∘𝜙, 𝑗=−𝑀,…,𝑁. Moreover, by some simple calculations, we have 𝑤^∗_𝑗(𝑢_𝑘)=𝑆(𝑗,ℎ)∘𝜙(𝑢_𝑘), 𝑗,𝑘=−𝑀,…,𝑁 with 𝜙 and 𝑢_𝑘 being defined in (33) and (41), respectively.

Having replaced the integral term on the left-hand side of system (38) with the right-hand side of (42) and having substituted 𝑢=𝑢_𝑘 for 𝑘=−𝑀,…,𝑁, where 𝑢_𝑘 are sinc grid points and also by replacing 𝑊(𝑢) by 𝑊(𝑢) as in (43), we obtain the following collocation result: 𝑊^′′𝑢_𝑘+𝜇𝑢_𝑘𝑊^′𝑢_𝑘−2𝜎²(𝜆+𝛿)𝑊𝑢_𝑘+2𝜎²𝜆^𝑁𝑁_{𝑙=−𝑀𝑖=−𝑀}𝑆𝑗,ℎ∘𝜙𝑢_𝑘𝐴_𝑖𝑙𝑊_𝑙+𝑅𝑢_𝑘=0, (44) where 𝑊(𝑢_𝑘)=∑^𝑁_𝑙=−𝑀𝑊_𝑙𝑆(𝑙,ℎ)∘𝜙(𝑢_𝑘).

Then, by Theorem A.7, we have 𝑊𝑢_𝑘=^𝑁_𝑙=−𝑀𝑊_𝑙𝛿⁽⁰⁾_𝑙𝑘,𝑊^′𝑢_𝑘=^𝑁_𝑙=−𝑀𝑊_𝑙𝜙^′𝑢_𝑘ℎ⁻¹𝛿⁽¹⁾_𝑙𝑘,𝑊^′′𝑢_𝑘=^𝑁_𝑙=−𝑀𝑊_𝑙𝜙^′′𝑢_𝑘ℎ⁻¹𝛿⁽¹⁾_𝑙𝑘+𝜙^′𝑢_𝑘²ℎ⁻²𝛿⁽²⁾_𝑙𝑘. (45)

Replacing (45) in (44), we rewrite (44) as ^𝑁_𝑙=−𝑀𝑊_𝑙𝜙^′′𝑢_𝑘ℎ⁻¹𝛿⁽¹⁾_𝑙𝑘+𝜙^′𝑢_𝑘²ℎ⁻²𝛿⁽²⁾_𝑙𝑘+𝜇𝑢_𝑘𝜙^′𝑢_𝑘ℎ⁻¹𝛿⁽¹⁾_𝑙𝑘−2𝜎²(𝜆+𝛿)𝛿⁽⁰⁾_𝑙𝑘+2𝜎²𝜆×^𝑁_𝑖=−𝑀𝑆𝑗,ℎ∘𝜙𝑢_𝑘𝐴_𝑖𝑙=−𝑅𝑢_𝑘. (46)

Having multiplied the resulting equations by ℎ²/(𝜙^′(𝑢_𝑘))², we have ^𝑁_𝑙=−𝑀𝑊_𝑙𝛿⁽²⁾_𝑙𝑘𝜙+ℎ^′′𝑢_𝑘𝜙^′𝑢_𝑘²+𝜇𝑢_𝑘𝜙^′𝑢_𝑘𝛿⁽¹⁾_𝑙𝑘−2ℎ²(𝜆+𝛿)𝜎²𝜙^′𝑢_𝑘²𝛿⁽⁰⁾_𝑙𝑘+2𝜆ℎ²𝜎²𝜙^′𝑢_𝑘^{2𝑁𝑖=−𝑀}𝑆𝑗,ℎ∘𝜙𝑢_𝑘𝐴_𝑖𝑙ℎ=−²𝑅𝑢_𝑘𝜙^′𝑢_𝑘². (47)

Now, since 𝛿⁽⁰⁾_𝑘𝑙=𝛿⁽⁰⁾_𝑙𝑘, 𝛿⁽¹⁾_𝑘𝑙=−𝛿⁽¹⁾_𝑙𝑘, 𝛿⁽²⁾_𝑘𝑙=𝛿⁽²⁾_𝑙𝑘, and 𝜙^′′(𝑢_𝑘)/(𝜙^′(𝑢_𝑘))²=−(1/𝜙^′(𝑢_𝑘))^′, we obtain the collocation result as ^𝑁_𝑙=−𝑀𝑊_𝑙𝛿⁽²⁾_𝑘𝑙1+ℎ𝜙^′𝑢_𝑘^′−𝜇𝑢_𝑘𝜙^′𝑢_𝑘𝛿⁽¹⁾_𝑘𝑙−2ℎ²(𝜆+𝛿)𝜎²𝜙^′𝑢_𝑘²𝛿⁽⁰⁾_𝑘𝑙+2𝜆ℎ²𝜎²𝜙^′𝑢_𝑘²×^𝑁_𝑖=−𝑀𝑆𝑗,ℎ∘𝜙𝑢_𝑘𝐴_𝑖𝑙ℎ=−²𝑅𝑢_𝑘𝜙^′𝑢_𝑘²,𝑘=−𝑀,…,𝑁. (48)

We set 𝐼^(𝑚)=[𝛿^(𝑚)_𝑘𝑙] and 𝑚=−1,0,1,2, where 𝛿^(𝑚)_𝑘𝑙 denotes the (𝑘,𝑙)th element of the matrix 𝐼^(𝑚). Also, we denote Δ=[𝑆(𝑖,ℎ)∘𝜙(𝑢_𝑘)] and 𝐼^(𝑚), 𝑚=−1,0,1,2 are square matrices of order (𝑁+𝑀+1)×(𝑁+𝑀+1). So, the system (48) can be given in matrix form as 𝐵𝑊=𝐹, (49) where 𝑊=[𝑊_−𝑀,…,𝑊_𝑁]^𝑇,  𝐹=[−(ℎ²𝐻(𝑢_−𝑀)/(𝜙^′(𝑢_−𝑀))²),…,−(ℎ²𝐻(𝑢_𝑁)/(𝜙^′(𝑢_𝑁))²)]^𝑇, and 𝐵=𝐼⁽²⁾+ℎ𝐷_𝑚1𝜙^′′−𝜇𝜙^′𝐼⁽¹⁾−2ℎ²𝜎²(𝜆+𝛿)𝐷_𝑚1𝜙^′2𝐼⁽⁰⁾+2ℎ²𝜆𝜎²𝐷_𝑚1𝜙^′2Δ𝐴. (50)

The above linear system contains (𝑁+𝑀+1) equations with (𝑁+𝑀+1) unknown coefficient 𝑊_𝑘, 𝑘=−𝑀,…,𝑁. Solving this linear system, we obtain the approximate solution of the system as follows: 𝑊(𝑢)≈𝑊(𝑢)=^𝑁_𝑙=−𝑀𝑊_𝑙𝑆()(𝑢),𝑉()(𝑢)+𝑢𝑙,ℎ∘𝜙𝑢;𝑏=𝑊𝑐1+𝑢²−𝑐₁𝛿. (51)

4. Numerical Example

In this section, we consider some numerical samples to illustrate the performance of sinc method and investigate how much the values of 𝑉(𝑢;𝑏) are affected by interest 𝑟. An example is solved under the assumption that the claim size density is given by 𝑓(𝑦)=𝜂𝑒^−𝜂𝑦𝐼(𝑦>0).

Example 1. Let 𝑐₁=0.1, 𝑐₂=0.2, 𝜇=5, 𝛼=𝛽=𝜋/4, 𝜆=3, 𝑏=2, 𝜎=0.05, 𝛿=0.006, 𝛼=1/2, 𝛽=1/4, 𝑑=1/4000, and 𝑟=0.05, the result is shown in Figure 1.
From Figure 1, we see that 𝑉(𝑢;2) is an increasing function with respect to 𝑢.

Figure 1 

The expected present value of total dividends 𝑉(𝑢;2).

Appendix

The sinc method is a highly efficient numerical method developed by Stenger, the pioneer of this field, people in his school, and others [18–22]. It is widely used in various fields of numerical analysis such as interpolation, quadrature, approximation of transforms, and solution of integral, ordinary differential, and partial differential equations.

sinc methods are based on the use of the cardinal function, 𝐶(𝑓,ℎ), which is sinc expansion of function 𝑓, defined by 𝐶𝑓,ℎ(𝑥)=_𝑘∈ℕ𝑓()𝑥𝑘ℎsincℎ−𝑘,−∞<𝑥<∞, (A.1) where the step size ℎ>0, and the function sinc is defined on the whole real line, −∞<𝑥<∞, by (𝑧)=()sincsin𝜋𝑧𝜋𝑧,if𝑧≠0,1,if𝑧=0. (A.2)

For any ℎ>0, the translated sinc functions with evenly spaced nodes are given as 𝑆(𝑧)𝑧𝑗,ℎ=sincℎ−𝑗,𝑗=0,±1,±2,…. (A.3)

The sinc functions are cardinal for the interpolating points 𝑘ℎ in the sense that 𝑆𝑗,ℎ(𝑘ℎ)=𝛿⁽⁰⁾_𝑗𝑘=0,if𝑘≠𝑗,1,if𝑘=𝑗. (A.4)

Definition A.1 (see [19, P73, De. 1.5.2]). Let 𝜙 denote a smooth one-to-one transformation of an arc Γ with end-point 𝑡₁ and 𝑡₂ onto ℝ, such that 𝜙(𝑡₁)=−∞ and 𝜙(𝑡₂)=∞. Let 𝜓=𝜙⁻¹ denote the inverse map, so that Γ=𝑧∈ℂ∶𝑧=𝜓(𝑢),𝑢∈ℝ. (A.5)
Given 𝜙, 𝜓, and a positive number ℎ, define the sinc points 𝑧_𝑘 by 𝑧_𝑘=𝑧_𝑘(ℎ)=𝜓(𝑘ℎ),𝑘=0,±1,±2,…, (A.6) and a function 𝜌 by 𝜌(𝑧)=𝑒^𝜙(𝑧). (A.7)
Observe that 𝜌(𝑧) increases from 0 to ∞ as 𝑧 traverses Γ from 𝑡₁ to 𝑡₂.
Corresponding to positive numbers 𝛼, 𝛽, and 𝑑, let 𝐿𝛼_,𝛽_,𝑑(𝜙) denote the family of all functions 𝑓 defined on Γ for which 𝑓(𝑧)=𝒪𝜌(𝑧)𝛼,if𝑧⟶𝑡₁,𝒪𝜌(𝑧)⁻𝛽,if𝑧⟶𝑡₂, (A.8) such that the Fourier transform {𝑓∘𝜙⁻¹}^∼ satisfies the following relation: 𝑓∘𝜙^−1∼(𝜁)=𝒪𝑒^{−𝑑|𝜁|}, (A.9) for all 𝜁∈ℝ.
Another important family of functions is 𝑀𝛼_,𝛽_,𝑑(𝜙), with 0<𝛼,𝛽≤1, and 0<𝑑<𝜋. It consists of all those functions 𝑓 defined on Γ such that 𝛼𝑔=𝑓−𝐿𝑓∈𝐿_,𝛽_,𝑑(𝜙) and where 𝐿𝑓 is defined by 𝐿𝑓=(𝑓(𝑡₁)+𝜌(𝑥)𝑓(𝑡₂))/(1+𝜌(𝑥)).
Let 𝑁 and 𝑀 denote two positive integers, such that 𝑀=[𝛽𝑁/𝛼], or 𝑁=[𝛼𝑀/𝛽] where [⋅] denotes the greatest integer function. Moreover, we denote 𝑚=𝑀+𝑁+1.

Definition A.2. Given three positive integers 𝑁, 𝑀, and 𝑚 as above, let 𝐷_𝑚 and 𝑉_𝑚 denote linear operators acting on functions 𝑓 defined on Γ given by 𝐷_𝑚𝑓=diag𝑓𝑧_−𝑀,…,𝑓𝑧_𝑁,𝑉_𝑚𝑓=𝑓𝑧_−𝑀,…,𝑓𝑧_𝑁^𝑇, (A.10) where 𝑧_𝑗=𝜓(𝑗ℎ) denote the sinc points. Set ℎ=𝜋𝑑𝛽𝑁^1/2,𝛾_𝑗𝑤=𝑆𝑗,ℎ∘𝜙,𝑗=−𝑀,…,𝑁,_𝑗=𝑤^∗_𝑗=𝛾_𝑗𝑤,𝑗=−𝑀+1,…,𝑁−1,_−𝑀=1−1+𝜌^𝑁_{𝑗=−𝑀+1}𝛾_𝑗1+𝑒^𝑗ℎ,𝑤_𝑁=𝜌−1+𝜌^𝑁−1_𝑗=−𝑀𝑒^𝑗ℎ𝛾_𝑗1+𝑒^𝑗ℎ,𝑤^∗_−𝑀=1+𝑒^−𝑀ℎ𝑤_−𝑀,𝑤^∗_𝑁=1+𝑒^−𝑁ℎ𝑤_𝑁,ℰ_𝑁=𝑁^1/2𝑒^−(𝜋𝑑𝛽_𝑁)^1/2,Δ_𝑚=𝑤_−𝑀,…,𝑤_𝑁,Δ^∗_𝑚=𝑤^∗_−𝑀,…,𝑤^∗_𝑁,𝛿⁽⁻¹⁾_𝑖𝑗=12+₀^𝑖−𝑗sinc(𝑡)𝑑𝑡, (A.11) and we then define an 𝑚×𝑚 matrix 𝐼⁽⁻¹⁾ by 𝐼⁽⁻¹⁾=[𝛿⁽⁻¹⁾_𝑖𝑗], with 𝛿⁽⁻¹⁾_𝑖𝑗 denoting the (𝑖,𝑗)th element of 𝐼⁽⁻¹⁾. Finally, let ‖⋅‖ denote the uniform norm on Γ that is, 𝑓=sup_𝑥∈Γ𝑓(𝑥). (A.12)

Theorem A.3 (see [19, P85, Th. 1.5.13]). Let 𝛼𝑓∈𝑀_,𝛽_,𝑑(𝜙); then, as 𝑁→∞, 𝑓−Δ_𝑚𝑉_𝑚𝑓=𝒪ℰ_𝑁. (A.13) If 𝛼𝑓∈𝐿_,𝛽_,𝑑(𝜙), then it is convenient to take 𝑤_𝑗=𝛾_𝑗=𝑆(𝑗,ℎ)∘𝜙, 𝑗=−𝑀,…,𝑁.

Theorem A.4 (see [19, P87, Th. 1.5.14]). Let 𝛼𝑓∈𝑀_,𝛽_,𝑑(𝜙); then, ℎ𝜙^′𝑘𝑓^(𝑘)−Δ_𝑚^(𝑘)𝑉_𝑚𝑓=𝒪ℰ_𝑁,𝑁>1. (A.14)

Theorem A.5 (see [19, P93, Th. 1.5.19]). Let 𝑓/𝜙^′𝛼∈𝐿_,𝛽^,𝑑(𝜙); then, for all 𝑁>1, 𝒥⁺𝑓−𝒥⁺_𝑚𝑓=𝒪ℰ_𝑁,𝒥⁻𝑓−𝒥⁻_𝑚𝑓=𝒪ℰ_𝑁, (A.15) where 𝒥⁺𝑓(𝑥)=^𝑥_𝑡1𝑓(𝑡)𝑑𝑡,𝒥⁻𝑓(𝑥)=^𝑡2_𝑥𝑓(𝑡)𝒥𝑑𝑡,_+𝑚𝑓(𝑥)=Δ_𝑚(𝑥)𝐴⁺𝑉_𝑚𝑓,𝐴⁺=ℎ𝐼⁽⁻¹⁾𝐷_𝑚1𝜙^′,𝒥⁻_𝑚𝑓(𝑥)=Δ_𝑚(𝑥)𝐴⁻𝑉_𝑚𝑓,𝐴⁻=ℎ𝐼^(−1)𝑇𝐷_𝑚1𝜙^′, (A.16) with (𝐼⁽⁻¹⁾)^𝑇 denoting the transpose of 𝐼⁽⁻¹⁾.

Theorem A.6 (see [19, P95-96, Th. 1.5.20]). Let 𝑔/𝜙^′𝛼∈𝐿_,𝛽^,𝑑(𝜙); suppose that 𝐴^± can be diagonalized with 𝐴⁺=𝑋𝑆𝑋⁻¹ and with 𝐴⁻=𝑌𝑆𝑌⁻¹ and that 𝑆 is a diagonal matrix. If for some positive 𝑐^′ independent of 𝑁, one has 𝐹^′(𝑠)≤𝑐^′ for all ℜ𝑠≥0, then there exists a constant 𝐶 independent of 𝑁 such that 𝑝−𝐹𝒥⁺_𝑚𝑔≤𝐶ℰ_𝑁,𝑞−𝐹𝒥⁻_𝑚𝑔≤𝐶ℰ_𝑁, (A.17) where 𝐹(𝑠)=^𝑐′′₀𝑒^−𝑡/𝑠𝑓(𝑡)𝑝(𝑥)=𝑑𝑡,_𝑥𝑡1𝑓()𝑔(𝑡)𝑞(𝑥)=𝑥−𝑡𝑑𝑡,_𝑡2𝑥𝑓()𝑔(𝑡)𝐹𝒥𝑡−𝑥𝑑𝑡,_+𝑚𝑔(𝑥)=Δ^∗_𝑚(𝑥)𝑋𝐹(𝑆)𝑋⁻¹𝑉_𝑚𝐹𝒥𝑔,_−𝑚𝑔(𝑥)=Δ^∗_𝑚(𝑥)𝑌𝐹(𝑆)𝑌⁻¹𝑉_𝑚𝑔, (A.18)(A.19) with 𝑐^′′≥𝑡₂−𝑡₁.

Theorem A.7 (see [18, P106, Th. 4.1]). Let 𝜙 be a conformal one-to-one transformation of an arc Γ. Then, 𝛿⁽⁰⁾_𝑗𝑘≡𝑆𝑗,ℎ∘𝜙(𝑧)_{𝑧=𝑧𝑘}=𝛿0,𝑖𝑓𝑘≠𝑗,1,𝑖𝑓𝑘=𝑗,_(1)𝑗𝑘𝑑≡ℎ𝑆(𝑧𝑑𝜙𝑗,ℎ∘𝜙)_{𝑧=𝑧𝑘}=(−1)^𝑘−𝑗𝛿𝑘−𝑗,𝑖𝑓𝑘≠𝑗,0,𝑖𝑓𝑘=𝑗,_(2)𝑗𝑘≡ℎ²𝑑²𝑑𝜙²𝑆𝑗,ℎ∘𝜙(𝑧)|_{𝑧=𝑧𝑘}=−2(−1)^𝑘−𝑗𝑘−𝑗²−𝜋,𝑖𝑓𝑘≠𝑗,₂3,𝑖𝑓𝑘=𝑗. (A.20)

Acknowledgments

This work was supported in part by the National Natural Science Foundation of China under Grant nos. 61370096 and 61173012 and the Key Project of Natural Science Foundation of Hunan Province under Grant no. 12JJA005.