- Email: [email protected]
An operator splitting harmonic differential quadrature approach to solve Young’s model for life insurance risk
Insurance: Mathematics and Economics 51 (2012) 442–448
Contents lists available at SciVerse ScienceDirect
Insurance: Mathematics and Economics journal homepage: www.elsevier.com/locate/ime
An operator splitting harmonic differential quadrature approach to solve Young’s model for life insurance risk Luca Vincenzo Ballestra a , Massimiliano Ottaviani b , Graziella Pacelli b,∗ a
Dipartimento di Strategie Aziendali e Metodologie Quantitative, Seconda Università di Napoli, Corso Gran Priorato di Malta, 81043 Capua, Italy
b
Dipartimento di Management, Università Politecnica delle Marche, Piazza Martelli 8, 60121 Ancona, Italy
article
info
Article history: Received May 2012 Received in revised form June 2012 Accepted 30 June 2012 Keywords: Young’s model Life insurance Harmonic differential quadrature
abstract This paper is concerned with the numerical approximation of a mathematical model for life insurance risk that has been presented quite recently by Young (2007, 2008). In particular, such a model, which consists of a system of several non-linear partial differential equations, is solved using a new numerical method that combines an operator splitting procedure with the differential quadrature (DQ) finite difference scheme. This approach allows one to reduce the partial differential problems to systems of linear equations of very small dimension, so that pricing portfolios of many life insurances becomes a relatively easily task. Numerical experiments are presented showing that the method proposed is very accurate and fast. In addition, the limit behavior of portfolios of life insurances as the number of contracts tends to infinity is investigated. This analysis reveals that the prices of portfolios comprising more than five thousand policies can be accurately approximated by solving a linear partial differential equation derived in Young (2007, 2008). © 2012 Elsevier B.V. All rights reserved.
1. Introduction We propose a numerical method for pricing portfolios of life insurances based on a mathematical model that has been presented quite recently by Young (2007, 2008). This model is interesting for several reasons. First of all, the complex stochastic evolution of the mortality rate is taken into account; second, insurers are treated as rational agents, who minimize the risk of their portfolios; finally, the mortality risk premium required by insurance companies is suitably modeled. Note that Young’s model is actually an extension to life insurance risk of a framework for pricing pure endowments that has been previously developed by Milevsky et al. (2005). Furthermore, models similar to Young’s model have also been proposed for pricing portfolios of life insurances and pure endowments (Bayraktar and Young, 2007), and for pricing life annuities (Bayraktar et al., 2009). In Young (2007, 2008) Young considers a portfolio of homogeneous life insurances (i.e., life insurances with same maturity and same face value), and finds that the price of such a portfolio satisfies a system of non-linear partial differential equations of parabolic type. Precisely, we have as many partial differential equations as is the number of policies considered. In addition, in Young (2007, 2008) it is also proven that, as the number of policies
∗
Corresponding author. Tel.: +39 071 2207050; fax: +39 071 2207150. E-mail address: [email protected] (G. Pacelli).
0167-6687/$ – see front matter © 2012 Elsevier B.V. All rights reserved. doi:10.1016/j.insmatheco.2012.06.012
considered tends to infinity, the average price of the life insurances tends to a function that satisfies a linear partial differential equation of parabolic type. Neither the non-linear system of partial differential equations nor the linear partial differential equation have exact closed-form solutions, so they must be approached by numerical techniques. Nevertheless, to the best of our knowledge, no numerical methods to solve the system of non-linear partial differential equations have yet been proposed. Instead, in Young (2008), a finite difference scheme is employed for approximating the linear partial differential equation that holds in the limit as the number of life insurances tends to infinity. Now, the solution of such a limit equation can be used as an estimation of the true portfolio price if the number of policies considered is sufficiently large. However, it would also be interesting to solve the system of non-linear partial differential equations that holds in the finite case, at least for two reasons. First of all, there are practical applications in which the number of life insurances to be priced is not extremely large (say smaller than some thousands). To figure this out, let us think, for instance, to the case of an insurer who wants to sell a basket of new contracts, for example to the employees of a firm; in addition, we shall also consider that Young’s model is to be applied not to all the policies traded by an insurer, but only to policies that can be reasonably assumed to be homogeneous (for maturity and face value). Second, solving the non-linear system of partial differential equations is useful to investigate how well the average life insurance price is approximated by the solution of the linear partial
L.V. Ballestra et al. / Insurance: Mathematics and Economics 51 (2012) 442–448
differential equation. In this respect it should be observed that in Young (2007, 2008) Young does not give any quantitative measure of the distance between the solution of the linear limit equation and the solution of the system of non-linear partial differential equations. In the present paper, the system of non-linear partial differential equations arising in Young’s model is solved using a new method based on the harmonic differential quadrature (HDQ) finite difference scheme. The HDQ approach, originally introduced by Striz et al. (1995), and further developed and analyzed by Chen et al. (1996) and Shu and Xue (1997), has been applied to various problems in science and engineering (see, for example, Janghorban (2011), Malekzadeh and Karami (2005), Shu (2000), Shu and Richard (1992) and Striz et al. (1997)). Its main advantage is that very accurate approximations are achieved using a computational mesh with a very small number of nodes. In the present manuscript, the HDQ method is employed in conjunction with a suitable operator splitting technique, which allows us to decouple the non-linear partial differential problem characteristic of Young’s model into independent smaller problems. By combining the HDQ method with such operator splitting procedure, we obtain systems of linear (algebraic) equations of very small dimension, so that portfolio prices can be calculated in reasonable times also when the number of life insurances is relatively large (say equal to some thousands). Numerical experiments are presented showing that the method proposed is significantly accurate and fast. In fact, for example, portfolios of five thousand policies can be priced with an error of order 10−3 in a time equal to 69.2 s (on a computer with a Pentium P6000 1.87 GHz 4 GB RAM). Furthermore, in this paper the convergence behavior of the portfolio price as the number of policies increases is investigated. In particular, it is shown that when there are five thousand life insurances the difference between the solution of the linear partial differential equation and the true average life insurance price starts to become of order 10−3 . Thus, the prices of portfolios that comprise at least five thousand policies can be accurately approximated by solving the linear partial differential equation, in place of the system of (many) non-linear partial differential equations. Finally, we point out that the numerical method proposed in the present manuscript could also be applied to several other models in actuarial mathematics, which, from the analytical standpoint, are substantially analogous to Young’s model. Some of these models are briefly described in Appendix. The remainder of the paper is organized as follows: in Section 2 the basic facts about Young’s model are briefly recalled; in Section 3 the HDQ-operator splitting finite difference approach is developed; in Section 4 the numerical results obtained are presented and discussed; in Section 5 some conclusions are drawn.
According to Young’s model (Young, 2007, 2008), the mortality rate λ, which measures the probability of an individual dying in an infinitesimal time, is described by the following stochastic differential equation: dλ = µ(λ, t )(λ − λ)dt + σ (t )(λ − λ)dW λ ,
(1)
where µ and σ are suitable functions of their arguments (they will be specified later), λ is a positive constant parameter, and W λ is a Wiener standard process. Furthermore, in Young (2008) the interest rate r is modeled according to the following stochastic differential equation: dr = b(r , t )dt + d(r , t )dW r ,
where b and d are suitable functions of their arguments (they will be specified later), and W r is a Wiener standard process uncorrelated with W λ . Note that in Young (2007, 2008) the vector variable (r , λ) is assumed to vary in the following set:
Ω = [0, +∞) × [λ, +∞).
(2)
(3)
Let Ω denote the interior set of Ω , that is Ω = (0, +∞) × (λ, +∞). Moreover, let A(N ) (r , λ, t ) denote the price of a portfolio of N homogeneous life insurances with face value 1 and maturity T . As shown in Young (2007, 2008), the function A(N ) must satisfy the following system of non-linear partial differential equations:
∂ A(n) ∂ A(n) 1 ∂ 2 A(n) + b(r , t ) + d2 (r , t ) ∂t ∂r 2 ∂r2 (n) ∂A ∂ 2 A(n) 1 + µ(λ, t )(λ − λ) + σ 2 (t )(λ − λ)2 − rA(n) ∂λ 2 ∂λ2 − nλ(A(n) − A(n−1) − 1) (n) 2 2 ∂A = −α σ 2 (t )(λ − λ)2 + nλ A(n) − A(n−1) − 1 , ∂λ (r , λ) ∈ Ω , t < T , n = 1, 2, . . . , N ,
(4)
where α is the so-called Sharpe ratio, which can be thought of as a mortality risk premium (see Young (2008)). The system of partial differential equations (4) must be solved recursively for n = 1, 2, . . . , N, and the function A(0) needed to obtain A(1) is given by A(0) (r , λ, t ) = 0,
(r , λ) ∈ Ω , t ≤ T .
(5)
The differential equations (4) must be equipped with final condition: A(n) (r , λ, T ) = 0, Let p p(N ) =
(N )
A
N
(r , λ) ∈ Ω , n = 1, 2, . . . , N .
(6)
denote the average price of the N life insurances:
(N )
.
(7)
In Young (2007, 2008) it is also shown that as N tends to infinity p(N ) tends to a function p that satisfies the following linear partial differential equation:
∂p ∂p 1 ∂ 2p + b(r , t ) + d2 (r , t ) 2 ∂t ∂r 2 ∂r
∂p 1 ∂ 2p + σ 2 (t )(λ − λ)2 2 ∂λ 2 ∂λ (r , λ) ∈ Ω , t < T ,
+ (µ(λ, t ) − ασ (t )) (λ − λ) − rp − λ(p − 1) = 0,
(8)
with the final condition: p(r , λ, T ) = 0,
2. The mathematical model
443
(r , λ) ∈ Ω .
(9)
In Young (2007, 2008) the problem of which boundary conditions to apply to (4)–(6) and to (8)–(9) is not addressed. Therefore, given that, for economic reasons, the function A(n) is monotone in the λ and the r variables (see also Fig. 1), and is bounded from below and above (being positive and smaller than the total face value of n life insurances), we prescribe:
∂ A(n) (r , λ, t ) = 0, r →+∞ ∂r λ ∈ [λ, +∞), t < T , n = 1, 2, . . . , N , ∂ A(n) (r , λ, t ) lim = 0, λ→+∞ ∂λ r ∈ [0, +∞), t < T , n = 1, 2, . . . , N . lim
(10)
(11)
444
L.V. Ballestra et al. / Insurance: Mathematics and Economics 51 (2012) 442–448
Analogously, we set:
∂ p(r , λ, t ) = 0, λ ∈ [λ, +∞), t < T , ∂r ∂ p(r , λ, t ) lim = 0, r ∈ [0, +∞), t < T . λ→+∞ ∂λ lim
r →+∞
(12) (13)
Instead, it is not straightforward to determine which boundary conditions (if any) should be applied at r = 0 and λ = λ. Therefore, in this paper, following an approach that is very common in quantitative finance (see for example Chiarella et al. (2009)), we assume that the partial differential equations (4) and (8) hold also for r = 0 and λ = λ, that is we let the partial differential equations (4) and (8) themselves impose their boundary conditions at r = 0 and λ = λ. In other words, we assume that Eqs. (4) and (8) hold for (r , λ) ∈ Ω . 3. The numerical method For the sake of brevity, we describe the numerical method directly by considering the partial differential problems (4), (6), (10), (11). The discretization of problem (8), (9), (12), (13) is substantially analogous, and thus is left to the reader. Furthermore, the approach we use is identical for every integer n ≥ 1. Therefore, from now on, for the sake of simplicity, focusing our attention on a single (generic) value of n, the notation n = 1, 2, . . . , N will be omitted. First of all, let us discretize problems (4), (6), (10), (11) in time. In the interval [0, T ] let us consider Mt + 1 equally spaced time levels tk = (k − 1)∆t, k = 1, 2, . . . , Mt + 1, where ∆t = MT . t
Let A(n),k (r , λ) denote an approximate value of A(n) (r , λ, t k ), k = 1, 2, . . . , Mt + 1. Eqs. (4), which, as already said, are assumed to hold for (r , λ) ∈ Ω , and the boundary conditions (10), (11), are approximated in time using the following operator splitting scheme: A(n),k − V (n),k
∆t
+ b(r , tk )
∂ V (n),k ∂r
∂ 2 V (n),k + d (r , t k ) − rV (n),k = 0, 2 ∂r2 1
lim
∂V
2
(n),k
(r , λ) = 0, ∂r
r →+∞
V (n),k − A(n),k−1
∆t
(r , λ) ∈ Ω ,
(14) (15)
∂ A(n),k−1 ∂λ
∂ 2 A(n),k−1 + σ (tk−1 )(λ − λ) 2 ∂λ2 (n),k−1 (n−1),k−1 − nλ(A −A − 1) = (n),k 2 2 ∂A −α σ 2 (tk−1 )(λ − λ)2 + nλ A(n),k − A(n−1),k − 1 , ∂λ 1
2
2
(r , λ) ∈ Ω ,
(16)
∂ A(n),k (r , λ) = 0, λ→+∞ ∂λ lim
r ∈ [0, +∞),
(17)
(r , λ) ∈ Ω .
3.1. Harmonic differential quadrature (HDQ) In the following, for the sake of brevity, we describe the harmonic differential quadrature method directly by considering problem (14)–(15). The numerical approximation of problem (16)–(17) is substantially analogous, and therefore is omitted. Furthermore, we will focus our attention on the single (generic) k-th time-step, and the fact that the index k goes back from Mt + 1 to 2 will be understood. First of all, the semi-infinite domain (3) is replaced with a bounded one:
(18)
That is, starting from the knowledge of A(n),k (r , λ), first of all we compute V (n),k (r , λ) by solving problem (14)–(15), and then we compute A(n),k−1 (r , λ) by solving problem (16)–(17), k = Mt + 1, Mt , . . . , 2. Note that Eq. (14) has been obtained by applying the implicit Euler method to Eq. (4) and by keeping only the time
(19)
where rmax and λmax will be chosen large enough such that the probability of r being larger than rmax and the probability of λ being greater than λmax are negligible (see the next section). Analogously, the boundary condition (15) is replaced by the following one:
∂ V (n),k (r , λ) = 0, ∂r r =rmax
λ ∈ [λ, λmax ].
(20)
In the interval [0, rmax ] let us consider Mr + 1 Chebyshev– Gauss–Lobatto points (Shu and Richard, 1992): ri =
rmax 2
i−1
1 − cos
Mr
π
,
i = 1, 2, . . . , Mr + 1,
(21)
and in the interval [λ, λmax ] let us consider Mλ + 1 Chebyshev– Gauss–Lobatto points:
λj = λ +
for k = Mt + 1, Mt , . . . , 2, with final condition: A(n),Mt +1 (r , λ) = 0,
Remark 1. The time discretization scheme (14)–(17), being based on the implicit/explicit Euler time stepping, is first-order accurate in time, only. Therefore, to gain accuracy, we employ the Richardson extrapolation procedure with halved time-step. More precisely, a second-order accurate solution is obtained by linear extrapolation of the solutions obtained using time-steps ∆t and ∆t (see for example Chang et al. (2007)). 2
Ωb = [0, rmax ] × [λ, λmax ],
λ ∈ [λ, +∞),
+ µ(λ, tk−1 )(λ − λ)
derivative, the terms that contain the derivatives with respect to r, and the term rA(n) . Instead, Eq. (16) has been obtained by applying a mixed implicit/explicit time-discretization scheme to Eq. (4) and by keeping only the time derivative, the terms that contain the derivatives with respect to λ, the term −nλ(A(n) − A(n−1) − 1), and the term at the right hand side. In particular, the terms that contain the derivatives with respect to λ and the term −nλ(A(n) − A(n−1) − 1) are treated according to the implicit Euler approach, whereas the term at the right hand side of (4), being non-linear, is treated according to the explicit Euler approach. The operator splitting scheme (14)–(17) allows us to solve Eq. (4) componentwise. Precisely, Eq. (14) is to be discretized in the r direction only, and Eq. (16) is to be discretized in the λ direction only, which will contribute to obtain systems of linear equations of very small dimension (see the next section).
λmax − λ 2
1 − cos
j−1 Mλ
π
,
j = 1, 2, . . . , Mλ + 1.
(22)
(n),k (n),k (n),k Let Ai,j , Vi,j , ∂ V∂ r , i,j
of A(n),k (ri , λj ), V (n),k (ri , λj ),
∂ 2 V (n),k denote approximate values ∂r2 i,j
∂ V (n),k (r ,λ) ∂r r =r ,λ=λ i
j
,
∂ 2 V (n),k (r ,λ) ∂r2 r =r ,λ=λ
respectively, i = 1, 2, . . . , Mr + 1, j = 1, 2, . . . , Mλ + 1.
i
j
,
L.V. Ballestra et al. / Insurance: Mathematics and Economics 51 (2012) 442–448
The key point of the HDQ method is that the first- and secondorder space derivatives of V (n),k are computed as follows (see Shu and Xue (1997) and Shu et al. (2002)):
Table 1 Model parameters and data.
γ θ κ µ λ σ α
M r +1 ∂ V (n),k = βi,l Vl(,nj ),k , ∂ r i,j l =1 i = 1, 2, . . . , Mr + 1, j = 1, 2, . . . , Mλ + 1,
(23)
(n),k
∂ V ∂r2 2
M r +1 = γi,l Vl(,nj ),k , i,j
l=1
(24)
where
π P (ri ) , r −r 2rmax P (rl ) sin 2ri l π max
i = 1, 2, . . . , Mr + 1, l = 1, 2, . . . , Mr + 1, i ̸= l,
(25)
Mr + 1
βi,i = −
βi,l ,
i = 1, 2, . . . , Mr + 1,
(26)
l=1,l̸=i
ri − rl π cot π , γi,l = βi,l 2βi,l − rmax
(27)
and Mr + 1
sin
ri − rl 2
l=1,l̸=i
π ,
i = 1, 2, . . . , Mr + 1.
(28)
Let us define:
hi,l (tk−1 ) = − b(ri , tk−1 )βi,l +
1 2
d2 (ri , tk−1 )γi,l
∆t ,
i = 1, 2, . . . , Mr , l = 1, 2, . . . , Mr + 1, i ̸= l, hi,i (tk−1 ) = 1 −
b(ri , tk−1 )βi,i +
1 2
d2 (ri , tk−1 )γi,i − ri
.
.. H k−1 = h (t ) Mr ,1 k−1
βMr +1,1
(n),k
Vj
4. Numerical results
∆t , (30)
h1,2 (tk−1 ) h2,2 (tk−1 )
··· ···
.. . hMr ,2 (tk−1 ) β M r + 1 ,2
h1,Mr +1 (tk−1 ) h2,Mr +1 (tk−1 )
.. .
··· ··· T
= V1(,nj),k−1 , V2(,nj),k−1 , . . . , VM(nr),+k1−,j1
, (31) hMr ,Mr +1 (tk−1 ) βMr +1,Mr +1 ,
j = 1, 2, . . . , Mλ + 1,
(32)
T
φj(n),k = A(1n,j),k , A(2n,j),k , . . . , A(Mnr),,jk , 0
,
j = 1, 2, . . . , Mλ + 1.
(33)
By evaluating Eq. (14) and the boundary condition (20) at r = ri and λ = λj , i = 1, 2, . . . , Mr , j = 1, 2, . . . , Mλ + 1, and by using (32)–(33) we obtain: (n),k
Hjk−1 Vj
= φj(n),k ,
j = 1, 2, . . . , Mλ + 1.
Remark 2. We do not perform a theoretical analysis on the stability of the numerical method developed in this paper, as it would be a very complicated task due to the presence of the nonlinear term at the right-hand side of (4), and to the different kinds of approximation involved (the operator-splitting approach, the Richardson extrapolation, the HDQ method). However, the numerical algorithm proposed has been directly tested by numerical simulations. In particular, we have performed several numerical experiments, and no kind of numerical instability has ever been observed. Some of these numerical experiments are reported in the next section.
(29)
i = 1, 2, . . . , Mr , h1,1 (tk−1 ) h2,1 (tk−1 )
The HDQ scheme employed allows us to obtain high levels of spatial accuracy using very small values of Mr and Mλ , so that the (n),k vectors Vj , j = 1, 2, . . . , Mλ + 1, can be computed extremely quickly. As already mentioned, problem (16)–(17) is solved using an analogous approach, and thus portfolios of many life insurances can be priced with a relatively small computational effort. The discretization of problem (8), (9), (12), (13) is substantially identical to that described above, and therefore is omitted. Here we simply observe that by applying the HDQ scheme to problem (8), (9), (12), (13) we obtain a set of Mr + 1 squared systems of Mλ + 1 linear equations analogous to (34) (which again are solved by Gaussian elimination with partial pivoting).
2rmax
i = 1, 2, . . . , Mr + 1,
P (ri ) =
0.1 year−1 0.03 year−1 0.1 year−1 0.04 year−1 0.02 year−1 0.1 year−1/2 0.1 10 year
T
i = 1, 2, . . . , Mr + 1, j = 1, 2, . . . , Mλ + 1,
βi,l =
445
(34)
Relations (34) constitute a set of systems of linear equations (n),k that allow us to determine the unknown vectors Vj , j = 1, 2, . . . , Mλ + 1. Note that we have Mλ + 1 squared linear systems each one of dimension Mr + 1. In this paper the linear systems (34) are solved by Gaussian elimination with partial pivoting (see Quarteroni et al. (2000)).
In Cox et al. (1985), the functions b(r , t ) and d(r , t ) modeling the dynamics of the interest rate are left generic. Therefore, using a very common approach, in the numerical simulations that follow, b(r , t ) and d(r , t ) are chosen according to the Cox–Ingersoll–Ross (CIR) model (see Cox et al. (1985)): b(r , t ) = −γ (r − θ ),
d(r , t ) = κ
√
r,
(35)
where γ , θ and κ are positive constants. Moreover, following (Young, 2008), the parameters µ, λ and σ are assumed to be constant:
µ(λ, t ) = µ,
σ (t ) = σ .
(36)
The values of model parameters and data used in the numerical simulations are reported in Table 1. Note that the values of µ, σ , α , T are the same considered by Young (2008) for solving the linear equation that holds in the limit case. Moreover, we set λmax = 0.4 and rmax = 0.4, since, as is reasonable to expect, both the interest rate and the mortality rate take values larger than 0.4 with extremely small probabilities. The relative error on the current (t = 0) price of the portfolio of life insurances obtained using the HDQ method is measured using the following maximum norm:
(N ),1 (N ) A i,j − A (ri , λj , 0) ErrRel = max , i ,j A(N ) (ri , λj , 0) i = 1, 2, . . . , Mr + 1, j = 1, 2, . . . , Mλ + 1,
(37)
446
L.V. Ballestra et al. / Insurance: Mathematics and Economics 51 (2012) 442–448
Table 2 Results obtained using the HDQ method, Mr = 10, Mλ = 10, Mt = 200. N
ErrRel
10 50 100 500 1000 2000 5000
2.17 × 10 1.57 × 10−3 1.41 × 10−3 1.18 × 10−3 1.14 × 10−3 1.66 × 10−3 3.31 × 10−3 −3
Table 4 Convergence of A(N ) to p.
CPUTime (s)
N
εN
0.109 0.951 1.50 7.50 14.6 26.9 69.2
10 50 100 500 1000 2000 5000
1.76 × 10−1 8.80 × 10−2 6.39 × 10−2 2.97 × 10−2 2.11 × 10−2 1.50 × 10−2 9.43 × 10−3
Table 3 Results obtained using the HDQ method, Mr = 20, Mλ = 20, Mt = 400. N
ErrRel
CPUTime (s)
10 50 100 500 1000 2000 5000
1.62 × 10−4 2.17 × 10−4 2.37 × 10−4 2.66 × 10−4 2.71 × 10−4 2.73 × 10−4 3.73 × 10−4
1.18 6.59 12.9 62.2 117 229 584
(Mt = 400), the portfolio price is computed with a relative error of order 10−4 (see Table 3). Furthermore, the method proposed is computationally fast, as the computer time required to price portfolios comprising up to 5000 policies (with an error of order 10−3 ) is smaller or equal to 69.2 s. In Fig. 1 the numerical solution obtained for a portfolio of five hundred policies (N = 500) is shown. We may see that the HDQ approach proposed allows us to obtain a sharp approximation of the portfolio price also at the boundaries of the computational domain, even though ‘‘artificial’’ boundary conditions have been prescribed (see Section 2). Incidentally, we observe that the portfolio price is monotonically increasing in the λ variable and monotonically decreasing in the r variable. This reflects the fact that the price of life insurance contracts increases as the mortality risk increases (as the probability of dying increases), and decreases as the interest rate increases (as it becomes more convenient to put money in a bank account). Convergence to the limit solution. We want to investigate the convergence of the average price of the N life insurance p(N ) (see (7)) to the function p solution of problem (8), (9), (12), (13). To this aim let us consider the discrete maximum norm of the (relative) difference between the (current) value of p(N ) and the (current) value of p:
ε
(N )
(N ) p (ri , λj , 0) − p(ri , λj , 0) , = max i ,j p(r , λ , 0) i
j
i = 1, 2, . . . , Mr + 1, j = 1, 2, . . . , Mλ + 1. Fig. 1. Solution obtained using the HDQ method, Mr = 10, Mλ = 10, Mt = 200.
(N ),1
where, according to the notation introduced in Section 3, Ai,j
is
the approximate value of A(N ) (ri , λj , 0) calculated using the HDQ method, i = 1, 2, . . . , Mr + 1, j = 1, 2, . . . , Mλ + 1. Note that the true solution A(N ) is not available, so an accurate estimation of it is computed by applying the HDQ scheme with (relatively) large values of Mt , Mr , Mλ (precisely we set Mt = 2000, Mr = 40, Mλ = 40). Moreover, the computer time required to obtain the HDQ numerical solution is denoted with CPUTime (the simulations are carried out on a computer with processor Intel Pentium P6000 1.87 GHz 4 GB RAM, and the software programs are written in FORTRAN). Performances of the HDQ method. In Tables 2 and 3 are reported the errors and computer times experienced by varying the number N of insurance policies and by using different values of Mt , Mr , Mλ . We may note that the HDQ scheme allows one to obtain a very accurate approximation of the price of the portfolio of life insurances. In fact using only 11 nodes in the r-direction (Mr = 10), only 11 nodes in the λ-direction (Mλ = 10), and 200 time-steps (Mt = 200), the portfolio price is computed with a relative error (measured in the maximum norm!) of order 10−3 (see Table 2). Instead of using 21 nodes in the r-direction (Mr = 20), 21 nodes in the λ-direction (Mλ = 20), and 400 time-steps
(38)
The values of ε (N ) experienced in correspondence of different values of N are reported in Table 4. We may note that, in accordance with the theory developed in Young (2008), ε (N ) tends monotonically to zero as N increases. In particular, ε (N ) starts to become of order 10−3 when N is approximately equal to 5000. Therefore, the prices of portfolios comprising more than five thousand life insurances can be accurately approximated by solving the linear partial differential problem (8), (9), (12), (13). This allows one to avoid to pass through the more complex system of non-linear partial differential equations (4), (6), (10), (11). On the contrary, for portfolios with less than five thousand policies, the use of the non-linear partial differential equations (to be solved by an appropriate numerical method such as the one proposed in this paper) seems to be preferable. 5. Conclusions Young’s model for risk insurance is solved using a new numerical method that combines an operator splitting procedure with the differential quadrature finite difference scheme. This approach allows us to reduce the systems of non-linear partial differential problems to systems of linear equations of very small dimension, so that pricing portfolios of many life insurances becomes a relatively easily task. Numerical experiments are presented showing that the method proposed is very accurate and fast. In addition, the limit behavior of portfolios of life insurances
L.V. Ballestra et al. / Insurance: Mathematics and Economics 51 (2012) 442–448
as the number of contracts tends to infinity is investigated. This analysis reveals that portfolios comprising more than five thousand policies can be priced by solving a linear partial differential equation derived in Young (2007, 2008), in place of the more complex system of non-linear partial differential problems. Finally, we point out that the numerical method developed in the present paper could also be applied to several other models in actuarial mathematics, which, from the mathematical standpoint, are substantially analogous to Young’s model. For the reader’s convenience, some of these models are described in Appendix. Appendix The numerical method developed in the present paper could also be used to solve a variety of models in actuarial mathematics, which, from the analytical standpoint, are substantially analogous to Young’s model. Below we briefly describe some of these models, namely the model by Bayraktar et al. (2009), the model by Ludkovski and Young (2008), and the model by Dahl (2004). For the reader’s convenience, we keep as much as possible the same notations used in the papers in which these models have originally been proposed. Therefore, in the following it may happen that the same symbol is used in different models and refers to different quantities, or that the same quantity is used in different models and is referred to by different symbols.
• The model by Bayraktar, Milevsky, Promislow and Young for life annuities In (Bayraktar et al., 2009), Bayraktar et al. consider a portfolio of N (homogeneous) life annuities that pay at a continuous rate of 1 per unit of time until maturity T or until the insured individual dies. Let a(N ) (r , λ, t ) denote the price of this portfolio, as a function of the interest rate r, of the mortality rate λ, and of the time t. According to the model developed in Bayraktar et al. (2009), a(N ) must be obtained solving the following system of partial differential equations:
∂ a(n) 1 ∂ 2 a(n) ∂ a(n) + b(r , t ) + d2 (r , t ) ∂t ∂r 2 ∂r2 ∂ a(n) 1 2 ∂ 2 a(n) + µ(λ, t )(λ − λ) + σ (t )(λ − λ)2 − ra(n) ∂λ 2 ∂λ2 − nλ(a(n) − a(n−1) ) + n (n) 2 2 ∂a + nλ a(n) − a(n−1) , = −α σ 2 (t )(λ − λ)2 ∂λ (r , λ) ∈ [0, +∞) × [λ, +∞), t < T , n = 1, 2, . . . , N ,
∂ H (n) 1 ∂ 2 H (n) ∂ H (n) + (a(t )r + b(t )) + (c (t )r + d(t )) ∂t ∂r 2 ∂r2 2 (n) ∂ H (n) ∂ H 1 2 + µ(λ, t )λ + σ (t )λ2 ∂λ 2 ∂λ2 (n) 2 γ ∂H 1 − rH (n) = − σ 2 (t )λ2 2 F (r , t ) ∂λ F (r , t ) − F (γr ,t ) H (n) −H (n−1) + nλ , 1−e γ (r , λ) ∈ [0, +∞) × [0, +∞), t < T , n = 1, 2, . . . , N , with final condition: H (n) (r , λ, T ) = n,
(r , λ) ∈ [0, +∞) × [0, +∞), n = 1, 2, . . . , N . The function H
(0)
needed in (A.1) to start the recursion is:
a(0) (r , λ, t ) = 0,
(r , λ) ∈ [0, +∞) × [λ, +∞), t ≤ T .
(A.5)
needed in (A.4) to start the recursion is
H (0) (r , λ, t ) = 0,
(r , λ) ∈ Ω , t ≤ T .
(A.6) (N ),a
Instead, the price of the N life annuities H must be obtained solving the following system of partial differential equations:
∂ H (n),a 1 ∂ 2 H (n),a ∂ H (n),a + (a(t )r + b(t )) + (c (t )r + d(t )) ∂t ∂r 2 ∂r2 (n),a 2 (n),a ∂H 1 ∂ H + µ(λ, t )λ + σ 2 (t )λ2 ∂λ 2 ∂λ2 (n),a 2 1 γ ∂H − rH (n),a = − σ 2 (t )λ2 2 F (r , t ) ∂λ γ ( n ), a ( n − F (r , t ) − H −H 1),a 1 − e F (r ,t ) − n, + nλ γ (A.7)
H (n),a (r , λ, T ) = 0, (A.1)
(r , λ) ∈ [0, +∞) × [0, +∞), n = 1, 2, . . . , N .
(A.2)
(A.3)
In (A.1)–(A.3), the quantities b, d, µ, σ , λ and α have exactly the same meaning as in Section 2. From the mathematical standpoint, the system of partial differential equations (A.1)–(A.3) is substantially identical to the system of partial differential equations (4)–(6), and thus could be solved using the HDQ-operator splitting approach developed in the present manuscript. 2
• The model by Ludkovski and Young for pure endowments and life annuities Let r and λ denote the interest rate and the mortality rate, respectively. Moreover, let H (N ) (r , λ, t ) denote the price of N
(A.8)
The function H (0),a needed in (A.7) to start the recursion is H (0) (r , λ, t ) = 0,
(0)
(A.4)
with final condition:
a(n) (r , λ, T ) = 0, The function a
(homogeneous) pure endowments with face value 1 and maturity T , and let H (N ),a (r , λ, t ) denote the value of N (homogeneous) life annuities that pay at a continuous rate of 1 per unit of time until maturity T or until the insured individuals die. In (Ludkovski and Young, 2008), Ludkovski and Young develop a model to compute H (N ) and H (N ),a based on indifference pricing. In particular, H (N ) must be obtained solving the following system of partial differential equations:
(r , λ) ∈ [0, +∞) × [0, +∞), t < T , n = 1, 2, . . . , N ,
with final condition:
(r , λ) ∈ [0, +∞) × [λ, +∞), n = 1, 2, . . . , N .
447
(r , λ) ∈ Ω , t ≤ T .
(A.9)
In (A.4), (A.7), a, b, c, d, σ , µ and F are suitable (known) deterministic functions of their arguments. In particular a, b, c and d model the dynamics of the (stochastic) interest rate, σ and µ model the dynamics of the (stochastic) mortality rate, and F is the market value of a risk-free bond. Finally, in (A.4), (A.7), γ is a constant parameter that models the utility of the insurer. From the mathematical standpoint, the systems of partial differential problems (A.4)–(A.9) are analogous to the system of partial differential problems (4)–(6), and thus could be solved using the HDQ-operator splitting approach developed in the present manuscript.
• The model by Dahl for market reserves
Let J denote a positive integer and let the set J = {1, 2, . . . , J } describe the possible states of a policyholder. As a simple case, it could be J = {1, 2}, with state 1 corresponding to the policyholder being alive, and state 2 corresponding to the policyholder being
448
L.V. Ballestra et al. / Insurance: Mathematics and Economics 51 (2012) 442–448
dead (however other states could be possible as well, for example a state telling if the policyholder is retired). In (Dahl, 2004), Dahl considers a mortality-linked insurance contract to which is associated a generic flow of payments that may depend on the state of the policyholder and on the price of a traded stock. Let s denote the price of the stock, and let µ denote the mortality rate. Finally, let V j (s, µ, t ) denote the market reserve for the insurance contract, which depends also on the state j of the policyholder, j ∈ J . According to the model developed in Dahl (2004), V 1 , V 2 , . . . , V J must satisfy the following system of partial differential equations:
∂ 2V j ∂V j 1 ∂V j + α µ,Q (µ, t ) + (σ µ (µ, t ))2 ∂t ∂µ 2 ∂µ2 2 j j ∂ V ∂V 1 2 + r (s, t )s + σ s (s, t ) s2 2 − r (s, t )V j ∂s 2 ∂s J = −aj (s, t ) − ajk (s, t ) + V k − V j k=1,k̸=j
× 1 + g (µ, t ) Rjk (µ, t ),
jk
(s, µ) ∈ [0, +∞) × [0, +∞), t < T , j = 1, 2, . . . , J ,
(A.10)
with final condition: V j (s, µ, T ) = ∆Aj (s, T ),
(s, µ) ∈ [0, +∞) × [0, +∞), j = 1, 2, . . . , J . µ,Q
µ
(A.11)
In (A.10)–(A.11), α , σ , r, σ , a , a , ∆A , g and Rjk are suitable (known) deterministic functions of their arguments, j = 1, 2, . . . , J, k = 1, 2, . . . , J. In particular α µ,Q and σ µ model the dynamics of the (stochastic) mortality rate, r (s, t ) and σ s model the dynamics of the (stochastic) stock price, aj , ajk , and ∆Aj models the flows of payments associated to the insurance contract, and g jk and Rjk model the transition from state j to state k, j = 1, 2, . . . , J, k = 1, 2, . . . , J. From the mathematical standpoint, the system of partial differential problems (A.10)–(A.11) is very similar to the system of partial differential problems (4)–(6), and thus could be solved using the HDQ-operator splitting approach proposed in this manuscript. Finally, we briefly recall that in (Dahl and Møller, 2006) Dahl and Møller propose a model to compute the market reserve value of a portfolio of pure endowments and derive a (single) partial differential equation analogous to the partial differential equations (A.10) (see Dahl and Møller (2006, p. 210)). Then, the numerical method developed in the present paper could be used to solve such a partial differential equation as well. s
j
jk
j
jk
References Bayraktar, E., Milevsky, M.A., Promislow, S.D., Young, V.R., 2009. Valuation of mortality risk via the instantaneous Sharpe ratio: applications to life annuities. Journal of Economic Dynamics and Control 33 (676–691). Bayraktar, E., Young, V.R., 2007. Hedging life insurance with pure endowments. Insurance: Mathematics and Economics 40, 435–444. Chang, C.-C., Chung, S.-L., Stapleton, R.C., 2007. Richardson extrapolation techniques for pricing American-style options. Journal of Futures Markets 27, 791–817. Chen, W., Wang, X., Zhong, T., 1996. The structure of weighting coefficient matrices of harmonic differential quadrature and its application. Communications in Numerical Methods in Engineering 12, 455–460. Chiarella, C., Kang, B., Meyer, G.H., Ziogas, A., 2009. The evaluation of American option prices under stochastic volatility and jump-diffusion dynamics using the method of lines. International Journal of Theoretical and Applied Finance 12, 393–425. Cox, J.C., Ingersoll, J.E., Ross, S.A., 1985. A theory of the term structure of interest rates. Econometrica 53, 385–408. Dahl, M., 2004. Stochastic mortality in life insurance: market reserves and mortality-linked insurance contracts. Insurance: Mathematics and Economics 35 (113–136). Dahl, M., Møller, T., 2006. Valuation and hedging of life insurance liabilities with systematic mortality risk. Insurance: Mathematics and Economics 39 (193–217). Janghorban, M., 2011. Two different types of differential quadrature methods for static analysis of microbeams based on nonlocal thermal elasticity theory in thermal environment. Archive of Applied Mechanics 1–7. Ludkovski, M., Young, V.R., 2008. Indifference pricing of pure endowments and life annuities under stochastic hazard and interest rates. Insurance: Mathematics and Economics 42 (14–30). Malekzadeh, P., Karami, G., 2005. Polynomial and harmonic differential quadrature methods for free vibration of variable thickness thick skew plates. Journal of Engineering Structures 27, 1563–1574. Milevsky, M.A., Promislow, S.D., Young, V.R., 2005. Financial valuation of mortality risk via the instantaneous Sharpe ratio: applications to pricing pure endowments, Technical Report. Available at http://arxiv.org/abs/0705.1302 (accessed May 2012). Quarteroni, A., Sacco, R., Saleri, F., 2000. Numerical Mathematics. Springer-Verlag, New-York. Shu, C., 2000. Differential Quadrature and Its Application in Engineering. SpringerVerlag, London. Shu, C., Richard, B.E., 1992. Application of generalized differential quadrature to solve two-dimensional incompressible Navier–Stokes stress equations. International Journal for Numerical Methods in Fluids 15, 791–798. Shu, C., Xue, H., 1997. Explicit computation of weighting coefficients in the harmonic differential quadrature. Journal of Sound and Vibration 204, 549–555. Shu, C., Yao, Q., Yeom, K.S., Zhu, Y.D., 2002. Numerical analysis of flow and thermal fields in arbitrary eccentric annulus by differential quadrature method. Heat and Mass Transfer 38, 597–608. Striz, A.G., Chen, W.L., Bert, C.W., 1997. Free vibration of plates by the high accuracy quadrature element method. Journal of Sound and Vibration 202, 689–702. Striz, A.G., Wang, X., Bert, C.W., 1995. Harmonic differential quadrature method and applications to analysis of structural components. ACTA Mechanica 111, 85–94. Young, V.R., 2007. Pricing life insurance under stochastic mortality via the instantaneous Sharpe ratio: theorems and proofs, Technical Report. Available at http://arxiv.org/abs/0705.1297 (accessed May 2012). Young, V.R., 2008. Pricing life insurance under stochastic mortality via the instantaneous Sharpe ratio. Insurance: Mathematics and Economics 42 (691–703).
Contents lists available at SciVerse ScienceDirect
Insurance: Mathematics and Economics journal homepage: www.elsevier.com/locate/ime
An operator splitting harmonic differential quadrature approach to solve Young’s model for life insurance risk Luca Vincenzo Ballestra a , Massimiliano Ottaviani b , Graziella Pacelli b,∗ a
Dipartimento di Strategie Aziendali e Metodologie Quantitative, Seconda Università di Napoli, Corso Gran Priorato di Malta, 81043 Capua, Italy
b
Dipartimento di Management, Università Politecnica delle Marche, Piazza Martelli 8, 60121 Ancona, Italy
article
info
Article history: Received May 2012 Received in revised form June 2012 Accepted 30 June 2012 Keywords: Young’s model Life insurance Harmonic differential quadrature
abstract This paper is concerned with the numerical approximation of a mathematical model for life insurance risk that has been presented quite recently by Young (2007, 2008). In particular, such a model, which consists of a system of several non-linear partial differential equations, is solved using a new numerical method that combines an operator splitting procedure with the differential quadrature (DQ) finite difference scheme. This approach allows one to reduce the partial differential problems to systems of linear equations of very small dimension, so that pricing portfolios of many life insurances becomes a relatively easily task. Numerical experiments are presented showing that the method proposed is very accurate and fast. In addition, the limit behavior of portfolios of life insurances as the number of contracts tends to infinity is investigated. This analysis reveals that the prices of portfolios comprising more than five thousand policies can be accurately approximated by solving a linear partial differential equation derived in Young (2007, 2008). © 2012 Elsevier B.V. All rights reserved.
1. Introduction We propose a numerical method for pricing portfolios of life insurances based on a mathematical model that has been presented quite recently by Young (2007, 2008). This model is interesting for several reasons. First of all, the complex stochastic evolution of the mortality rate is taken into account; second, insurers are treated as rational agents, who minimize the risk of their portfolios; finally, the mortality risk premium required by insurance companies is suitably modeled. Note that Young’s model is actually an extension to life insurance risk of a framework for pricing pure endowments that has been previously developed by Milevsky et al. (2005). Furthermore, models similar to Young’s model have also been proposed for pricing portfolios of life insurances and pure endowments (Bayraktar and Young, 2007), and for pricing life annuities (Bayraktar et al., 2009). In Young (2007, 2008) Young considers a portfolio of homogeneous life insurances (i.e., life insurances with same maturity and same face value), and finds that the price of such a portfolio satisfies a system of non-linear partial differential equations of parabolic type. Precisely, we have as many partial differential equations as is the number of policies considered. In addition, in Young (2007, 2008) it is also proven that, as the number of policies
∗
Corresponding author. Tel.: +39 071 2207050; fax: +39 071 2207150. E-mail address: [email protected] (G. Pacelli).
0167-6687/$ – see front matter © 2012 Elsevier B.V. All rights reserved. doi:10.1016/j.insmatheco.2012.06.012
considered tends to infinity, the average price of the life insurances tends to a function that satisfies a linear partial differential equation of parabolic type. Neither the non-linear system of partial differential equations nor the linear partial differential equation have exact closed-form solutions, so they must be approached by numerical techniques. Nevertheless, to the best of our knowledge, no numerical methods to solve the system of non-linear partial differential equations have yet been proposed. Instead, in Young (2008), a finite difference scheme is employed for approximating the linear partial differential equation that holds in the limit as the number of life insurances tends to infinity. Now, the solution of such a limit equation can be used as an estimation of the true portfolio price if the number of policies considered is sufficiently large. However, it would also be interesting to solve the system of non-linear partial differential equations that holds in the finite case, at least for two reasons. First of all, there are practical applications in which the number of life insurances to be priced is not extremely large (say smaller than some thousands). To figure this out, let us think, for instance, to the case of an insurer who wants to sell a basket of new contracts, for example to the employees of a firm; in addition, we shall also consider that Young’s model is to be applied not to all the policies traded by an insurer, but only to policies that can be reasonably assumed to be homogeneous (for maturity and face value). Second, solving the non-linear system of partial differential equations is useful to investigate how well the average life insurance price is approximated by the solution of the linear partial
L.V. Ballestra et al. / Insurance: Mathematics and Economics 51 (2012) 442–448
differential equation. In this respect it should be observed that in Young (2007, 2008) Young does not give any quantitative measure of the distance between the solution of the linear limit equation and the solution of the system of non-linear partial differential equations. In the present paper, the system of non-linear partial differential equations arising in Young’s model is solved using a new method based on the harmonic differential quadrature (HDQ) finite difference scheme. The HDQ approach, originally introduced by Striz et al. (1995), and further developed and analyzed by Chen et al. (1996) and Shu and Xue (1997), has been applied to various problems in science and engineering (see, for example, Janghorban (2011), Malekzadeh and Karami (2005), Shu (2000), Shu and Richard (1992) and Striz et al. (1997)). Its main advantage is that very accurate approximations are achieved using a computational mesh with a very small number of nodes. In the present manuscript, the HDQ method is employed in conjunction with a suitable operator splitting technique, which allows us to decouple the non-linear partial differential problem characteristic of Young’s model into independent smaller problems. By combining the HDQ method with such operator splitting procedure, we obtain systems of linear (algebraic) equations of very small dimension, so that portfolio prices can be calculated in reasonable times also when the number of life insurances is relatively large (say equal to some thousands). Numerical experiments are presented showing that the method proposed is significantly accurate and fast. In fact, for example, portfolios of five thousand policies can be priced with an error of order 10−3 in a time equal to 69.2 s (on a computer with a Pentium P6000 1.87 GHz 4 GB RAM). Furthermore, in this paper the convergence behavior of the portfolio price as the number of policies increases is investigated. In particular, it is shown that when there are five thousand life insurances the difference between the solution of the linear partial differential equation and the true average life insurance price starts to become of order 10−3 . Thus, the prices of portfolios that comprise at least five thousand policies can be accurately approximated by solving the linear partial differential equation, in place of the system of (many) non-linear partial differential equations. Finally, we point out that the numerical method proposed in the present manuscript could also be applied to several other models in actuarial mathematics, which, from the analytical standpoint, are substantially analogous to Young’s model. Some of these models are briefly described in Appendix. The remainder of the paper is organized as follows: in Section 2 the basic facts about Young’s model are briefly recalled; in Section 3 the HDQ-operator splitting finite difference approach is developed; in Section 4 the numerical results obtained are presented and discussed; in Section 5 some conclusions are drawn.
According to Young’s model (Young, 2007, 2008), the mortality rate λ, which measures the probability of an individual dying in an infinitesimal time, is described by the following stochastic differential equation: dλ = µ(λ, t )(λ − λ)dt + σ (t )(λ − λ)dW λ ,
(1)
where µ and σ are suitable functions of their arguments (they will be specified later), λ is a positive constant parameter, and W λ is a Wiener standard process. Furthermore, in Young (2008) the interest rate r is modeled according to the following stochastic differential equation: dr = b(r , t )dt + d(r , t )dW r ,
where b and d are suitable functions of their arguments (they will be specified later), and W r is a Wiener standard process uncorrelated with W λ . Note that in Young (2007, 2008) the vector variable (r , λ) is assumed to vary in the following set:
Ω = [0, +∞) × [λ, +∞).
(2)
(3)
Let Ω denote the interior set of Ω , that is Ω = (0, +∞) × (λ, +∞). Moreover, let A(N ) (r , λ, t ) denote the price of a portfolio of N homogeneous life insurances with face value 1 and maturity T . As shown in Young (2007, 2008), the function A(N ) must satisfy the following system of non-linear partial differential equations:
∂ A(n) ∂ A(n) 1 ∂ 2 A(n) + b(r , t ) + d2 (r , t ) ∂t ∂r 2 ∂r2 (n) ∂A ∂ 2 A(n) 1 + µ(λ, t )(λ − λ) + σ 2 (t )(λ − λ)2 − rA(n) ∂λ 2 ∂λ2 − nλ(A(n) − A(n−1) − 1) (n) 2 2 ∂A = −α σ 2 (t )(λ − λ)2 + nλ A(n) − A(n−1) − 1 , ∂λ (r , λ) ∈ Ω , t < T , n = 1, 2, . . . , N ,
(4)
where α is the so-called Sharpe ratio, which can be thought of as a mortality risk premium (see Young (2008)). The system of partial differential equations (4) must be solved recursively for n = 1, 2, . . . , N, and the function A(0) needed to obtain A(1) is given by A(0) (r , λ, t ) = 0,
(r , λ) ∈ Ω , t ≤ T .
(5)
The differential equations (4) must be equipped with final condition: A(n) (r , λ, T ) = 0, Let p p(N ) =
(N )
A
N
(r , λ) ∈ Ω , n = 1, 2, . . . , N .
(6)
denote the average price of the N life insurances:
(N )
.
(7)
In Young (2007, 2008) it is also shown that as N tends to infinity p(N ) tends to a function p that satisfies the following linear partial differential equation:
∂p ∂p 1 ∂ 2p + b(r , t ) + d2 (r , t ) 2 ∂t ∂r 2 ∂r
∂p 1 ∂ 2p + σ 2 (t )(λ − λ)2 2 ∂λ 2 ∂λ (r , λ) ∈ Ω , t < T ,
+ (µ(λ, t ) − ασ (t )) (λ − λ) − rp − λ(p − 1) = 0,
(8)
with the final condition: p(r , λ, T ) = 0,
2. The mathematical model
443
(r , λ) ∈ Ω .
(9)
In Young (2007, 2008) the problem of which boundary conditions to apply to (4)–(6) and to (8)–(9) is not addressed. Therefore, given that, for economic reasons, the function A(n) is monotone in the λ and the r variables (see also Fig. 1), and is bounded from below and above (being positive and smaller than the total face value of n life insurances), we prescribe:
∂ A(n) (r , λ, t ) = 0, r →+∞ ∂r λ ∈ [λ, +∞), t < T , n = 1, 2, . . . , N , ∂ A(n) (r , λ, t ) lim = 0, λ→+∞ ∂λ r ∈ [0, +∞), t < T , n = 1, 2, . . . , N . lim
(10)
(11)
444
L.V. Ballestra et al. / Insurance: Mathematics and Economics 51 (2012) 442–448
Analogously, we set:
∂ p(r , λ, t ) = 0, λ ∈ [λ, +∞), t < T , ∂r ∂ p(r , λ, t ) lim = 0, r ∈ [0, +∞), t < T . λ→+∞ ∂λ lim
r →+∞
(12) (13)
Instead, it is not straightforward to determine which boundary conditions (if any) should be applied at r = 0 and λ = λ. Therefore, in this paper, following an approach that is very common in quantitative finance (see for example Chiarella et al. (2009)), we assume that the partial differential equations (4) and (8) hold also for r = 0 and λ = λ, that is we let the partial differential equations (4) and (8) themselves impose their boundary conditions at r = 0 and λ = λ. In other words, we assume that Eqs. (4) and (8) hold for (r , λ) ∈ Ω . 3. The numerical method For the sake of brevity, we describe the numerical method directly by considering the partial differential problems (4), (6), (10), (11). The discretization of problem (8), (9), (12), (13) is substantially analogous, and thus is left to the reader. Furthermore, the approach we use is identical for every integer n ≥ 1. Therefore, from now on, for the sake of simplicity, focusing our attention on a single (generic) value of n, the notation n = 1, 2, . . . , N will be omitted. First of all, let us discretize problems (4), (6), (10), (11) in time. In the interval [0, T ] let us consider Mt + 1 equally spaced time levels tk = (k − 1)∆t, k = 1, 2, . . . , Mt + 1, where ∆t = MT . t
Let A(n),k (r , λ) denote an approximate value of A(n) (r , λ, t k ), k = 1, 2, . . . , Mt + 1. Eqs. (4), which, as already said, are assumed to hold for (r , λ) ∈ Ω , and the boundary conditions (10), (11), are approximated in time using the following operator splitting scheme: A(n),k − V (n),k
∆t
+ b(r , tk )
∂ V (n),k ∂r
∂ 2 V (n),k + d (r , t k ) − rV (n),k = 0, 2 ∂r2 1
lim
∂V
2
(n),k
(r , λ) = 0, ∂r
r →+∞
V (n),k − A(n),k−1
∆t
(r , λ) ∈ Ω ,
(14) (15)
∂ A(n),k−1 ∂λ
∂ 2 A(n),k−1 + σ (tk−1 )(λ − λ) 2 ∂λ2 (n),k−1 (n−1),k−1 − nλ(A −A − 1) = (n),k 2 2 ∂A −α σ 2 (tk−1 )(λ − λ)2 + nλ A(n),k − A(n−1),k − 1 , ∂λ 1
2
2
(r , λ) ∈ Ω ,
(16)
∂ A(n),k (r , λ) = 0, λ→+∞ ∂λ lim
r ∈ [0, +∞),
(17)
(r , λ) ∈ Ω .
3.1. Harmonic differential quadrature (HDQ) In the following, for the sake of brevity, we describe the harmonic differential quadrature method directly by considering problem (14)–(15). The numerical approximation of problem (16)–(17) is substantially analogous, and therefore is omitted. Furthermore, we will focus our attention on the single (generic) k-th time-step, and the fact that the index k goes back from Mt + 1 to 2 will be understood. First of all, the semi-infinite domain (3) is replaced with a bounded one:
(18)
That is, starting from the knowledge of A(n),k (r , λ), first of all we compute V (n),k (r , λ) by solving problem (14)–(15), and then we compute A(n),k−1 (r , λ) by solving problem (16)–(17), k = Mt + 1, Mt , . . . , 2. Note that Eq. (14) has been obtained by applying the implicit Euler method to Eq. (4) and by keeping only the time
(19)
where rmax and λmax will be chosen large enough such that the probability of r being larger than rmax and the probability of λ being greater than λmax are negligible (see the next section). Analogously, the boundary condition (15) is replaced by the following one:
∂ V (n),k (r , λ) = 0, ∂r r =rmax
λ ∈ [λ, λmax ].
(20)
In the interval [0, rmax ] let us consider Mr + 1 Chebyshev– Gauss–Lobatto points (Shu and Richard, 1992): ri =
rmax 2
i−1
1 − cos
Mr
π
,
i = 1, 2, . . . , Mr + 1,
(21)
and in the interval [λ, λmax ] let us consider Mλ + 1 Chebyshev– Gauss–Lobatto points:
λj = λ +
for k = Mt + 1, Mt , . . . , 2, with final condition: A(n),Mt +1 (r , λ) = 0,
Remark 1. The time discretization scheme (14)–(17), being based on the implicit/explicit Euler time stepping, is first-order accurate in time, only. Therefore, to gain accuracy, we employ the Richardson extrapolation procedure with halved time-step. More precisely, a second-order accurate solution is obtained by linear extrapolation of the solutions obtained using time-steps ∆t and ∆t (see for example Chang et al. (2007)). 2
Ωb = [0, rmax ] × [λ, λmax ],
λ ∈ [λ, +∞),
+ µ(λ, tk−1 )(λ − λ)
derivative, the terms that contain the derivatives with respect to r, and the term rA(n) . Instead, Eq. (16) has been obtained by applying a mixed implicit/explicit time-discretization scheme to Eq. (4) and by keeping only the time derivative, the terms that contain the derivatives with respect to λ, the term −nλ(A(n) − A(n−1) − 1), and the term at the right hand side. In particular, the terms that contain the derivatives with respect to λ and the term −nλ(A(n) − A(n−1) − 1) are treated according to the implicit Euler approach, whereas the term at the right hand side of (4), being non-linear, is treated according to the explicit Euler approach. The operator splitting scheme (14)–(17) allows us to solve Eq. (4) componentwise. Precisely, Eq. (14) is to be discretized in the r direction only, and Eq. (16) is to be discretized in the λ direction only, which will contribute to obtain systems of linear equations of very small dimension (see the next section).
λmax − λ 2
1 − cos
j−1 Mλ
π
,
j = 1, 2, . . . , Mλ + 1.
(22)
(n),k (n),k (n),k Let Ai,j , Vi,j , ∂ V∂ r , i,j
of A(n),k (ri , λj ), V (n),k (ri , λj ),
∂ 2 V (n),k denote approximate values ∂r2 i,j
∂ V (n),k (r ,λ) ∂r r =r ,λ=λ i
j
,
∂ 2 V (n),k (r ,λ) ∂r2 r =r ,λ=λ
respectively, i = 1, 2, . . . , Mr + 1, j = 1, 2, . . . , Mλ + 1.
i
j
,
L.V. Ballestra et al. / Insurance: Mathematics and Economics 51 (2012) 442–448
The key point of the HDQ method is that the first- and secondorder space derivatives of V (n),k are computed as follows (see Shu and Xue (1997) and Shu et al. (2002)):
Table 1 Model parameters and data.
γ θ κ µ λ σ α
M r +1 ∂ V (n),k = βi,l Vl(,nj ),k , ∂ r i,j l =1 i = 1, 2, . . . , Mr + 1, j = 1, 2, . . . , Mλ + 1,
(23)
(n),k
∂ V ∂r2 2
M r +1 = γi,l Vl(,nj ),k , i,j
l=1
(24)
where
π P (ri ) , r −r 2rmax P (rl ) sin 2ri l π max
i = 1, 2, . . . , Mr + 1, l = 1, 2, . . . , Mr + 1, i ̸= l,
(25)
Mr + 1
βi,i = −
βi,l ,
i = 1, 2, . . . , Mr + 1,
(26)
l=1,l̸=i
ri − rl π cot π , γi,l = βi,l 2βi,l − rmax
(27)
and Mr + 1
sin
ri − rl 2
l=1,l̸=i
π ,
i = 1, 2, . . . , Mr + 1.
(28)
Let us define:
hi,l (tk−1 ) = − b(ri , tk−1 )βi,l +
1 2
d2 (ri , tk−1 )γi,l
∆t ,
i = 1, 2, . . . , Mr , l = 1, 2, . . . , Mr + 1, i ̸= l, hi,i (tk−1 ) = 1 −
b(ri , tk−1 )βi,i +
1 2
d2 (ri , tk−1 )γi,i − ri
.
.. H k−1 = h (t ) Mr ,1 k−1
βMr +1,1
(n),k
Vj
4. Numerical results
∆t , (30)
h1,2 (tk−1 ) h2,2 (tk−1 )
··· ···
.. . hMr ,2 (tk−1 ) β M r + 1 ,2
h1,Mr +1 (tk−1 ) h2,Mr +1 (tk−1 )
.. .
··· ··· T
= V1(,nj),k−1 , V2(,nj),k−1 , . . . , VM(nr),+k1−,j1
, (31) hMr ,Mr +1 (tk−1 ) βMr +1,Mr +1 ,
j = 1, 2, . . . , Mλ + 1,
(32)
T
φj(n),k = A(1n,j),k , A(2n,j),k , . . . , A(Mnr),,jk , 0
,
j = 1, 2, . . . , Mλ + 1.
(33)
By evaluating Eq. (14) and the boundary condition (20) at r = ri and λ = λj , i = 1, 2, . . . , Mr , j = 1, 2, . . . , Mλ + 1, and by using (32)–(33) we obtain: (n),k
Hjk−1 Vj
= φj(n),k ,
j = 1, 2, . . . , Mλ + 1.
Remark 2. We do not perform a theoretical analysis on the stability of the numerical method developed in this paper, as it would be a very complicated task due to the presence of the nonlinear term at the right-hand side of (4), and to the different kinds of approximation involved (the operator-splitting approach, the Richardson extrapolation, the HDQ method). However, the numerical algorithm proposed has been directly tested by numerical simulations. In particular, we have performed several numerical experiments, and no kind of numerical instability has ever been observed. Some of these numerical experiments are reported in the next section.
(29)
i = 1, 2, . . . , Mr , h1,1 (tk−1 ) h2,1 (tk−1 )
The HDQ scheme employed allows us to obtain high levels of spatial accuracy using very small values of Mr and Mλ , so that the (n),k vectors Vj , j = 1, 2, . . . , Mλ + 1, can be computed extremely quickly. As already mentioned, problem (16)–(17) is solved using an analogous approach, and thus portfolios of many life insurances can be priced with a relatively small computational effort. The discretization of problem (8), (9), (12), (13) is substantially identical to that described above, and therefore is omitted. Here we simply observe that by applying the HDQ scheme to problem (8), (9), (12), (13) we obtain a set of Mr + 1 squared systems of Mλ + 1 linear equations analogous to (34) (which again are solved by Gaussian elimination with partial pivoting).
2rmax
i = 1, 2, . . . , Mr + 1,
P (ri ) =
0.1 year−1 0.03 year−1 0.1 year−1 0.04 year−1 0.02 year−1 0.1 year−1/2 0.1 10 year
T
i = 1, 2, . . . , Mr + 1, j = 1, 2, . . . , Mλ + 1,
βi,l =
445
(34)
Relations (34) constitute a set of systems of linear equations (n),k that allow us to determine the unknown vectors Vj , j = 1, 2, . . . , Mλ + 1. Note that we have Mλ + 1 squared linear systems each one of dimension Mr + 1. In this paper the linear systems (34) are solved by Gaussian elimination with partial pivoting (see Quarteroni et al. (2000)).
In Cox et al. (1985), the functions b(r , t ) and d(r , t ) modeling the dynamics of the interest rate are left generic. Therefore, using a very common approach, in the numerical simulations that follow, b(r , t ) and d(r , t ) are chosen according to the Cox–Ingersoll–Ross (CIR) model (see Cox et al. (1985)): b(r , t ) = −γ (r − θ ),
d(r , t ) = κ
√
r,
(35)
where γ , θ and κ are positive constants. Moreover, following (Young, 2008), the parameters µ, λ and σ are assumed to be constant:
µ(λ, t ) = µ,
σ (t ) = σ .
(36)
The values of model parameters and data used in the numerical simulations are reported in Table 1. Note that the values of µ, σ , α , T are the same considered by Young (2008) for solving the linear equation that holds in the limit case. Moreover, we set λmax = 0.4 and rmax = 0.4, since, as is reasonable to expect, both the interest rate and the mortality rate take values larger than 0.4 with extremely small probabilities. The relative error on the current (t = 0) price of the portfolio of life insurances obtained using the HDQ method is measured using the following maximum norm:
(N ),1 (N ) A i,j − A (ri , λj , 0) ErrRel = max , i ,j A(N ) (ri , λj , 0) i = 1, 2, . . . , Mr + 1, j = 1, 2, . . . , Mλ + 1,
(37)
446
L.V. Ballestra et al. / Insurance: Mathematics and Economics 51 (2012) 442–448
Table 2 Results obtained using the HDQ method, Mr = 10, Mλ = 10, Mt = 200. N
ErrRel
10 50 100 500 1000 2000 5000
2.17 × 10 1.57 × 10−3 1.41 × 10−3 1.18 × 10−3 1.14 × 10−3 1.66 × 10−3 3.31 × 10−3 −3
Table 4 Convergence of A(N ) to p.
CPUTime (s)
N
εN
0.109 0.951 1.50 7.50 14.6 26.9 69.2
10 50 100 500 1000 2000 5000
1.76 × 10−1 8.80 × 10−2 6.39 × 10−2 2.97 × 10−2 2.11 × 10−2 1.50 × 10−2 9.43 × 10−3
Table 3 Results obtained using the HDQ method, Mr = 20, Mλ = 20, Mt = 400. N
ErrRel
CPUTime (s)
10 50 100 500 1000 2000 5000
1.62 × 10−4 2.17 × 10−4 2.37 × 10−4 2.66 × 10−4 2.71 × 10−4 2.73 × 10−4 3.73 × 10−4
1.18 6.59 12.9 62.2 117 229 584
(Mt = 400), the portfolio price is computed with a relative error of order 10−4 (see Table 3). Furthermore, the method proposed is computationally fast, as the computer time required to price portfolios comprising up to 5000 policies (with an error of order 10−3 ) is smaller or equal to 69.2 s. In Fig. 1 the numerical solution obtained for a portfolio of five hundred policies (N = 500) is shown. We may see that the HDQ approach proposed allows us to obtain a sharp approximation of the portfolio price also at the boundaries of the computational domain, even though ‘‘artificial’’ boundary conditions have been prescribed (see Section 2). Incidentally, we observe that the portfolio price is monotonically increasing in the λ variable and monotonically decreasing in the r variable. This reflects the fact that the price of life insurance contracts increases as the mortality risk increases (as the probability of dying increases), and decreases as the interest rate increases (as it becomes more convenient to put money in a bank account). Convergence to the limit solution. We want to investigate the convergence of the average price of the N life insurance p(N ) (see (7)) to the function p solution of problem (8), (9), (12), (13). To this aim let us consider the discrete maximum norm of the (relative) difference between the (current) value of p(N ) and the (current) value of p:
ε
(N )
(N ) p (ri , λj , 0) − p(ri , λj , 0) , = max i ,j p(r , λ , 0) i
j
i = 1, 2, . . . , Mr + 1, j = 1, 2, . . . , Mλ + 1. Fig. 1. Solution obtained using the HDQ method, Mr = 10, Mλ = 10, Mt = 200.
(N ),1
where, according to the notation introduced in Section 3, Ai,j
is
the approximate value of A(N ) (ri , λj , 0) calculated using the HDQ method, i = 1, 2, . . . , Mr + 1, j = 1, 2, . . . , Mλ + 1. Note that the true solution A(N ) is not available, so an accurate estimation of it is computed by applying the HDQ scheme with (relatively) large values of Mt , Mr , Mλ (precisely we set Mt = 2000, Mr = 40, Mλ = 40). Moreover, the computer time required to obtain the HDQ numerical solution is denoted with CPUTime (the simulations are carried out on a computer with processor Intel Pentium P6000 1.87 GHz 4 GB RAM, and the software programs are written in FORTRAN). Performances of the HDQ method. In Tables 2 and 3 are reported the errors and computer times experienced by varying the number N of insurance policies and by using different values of Mt , Mr , Mλ . We may note that the HDQ scheme allows one to obtain a very accurate approximation of the price of the portfolio of life insurances. In fact using only 11 nodes in the r-direction (Mr = 10), only 11 nodes in the λ-direction (Mλ = 10), and 200 time-steps (Mt = 200), the portfolio price is computed with a relative error (measured in the maximum norm!) of order 10−3 (see Table 2). Instead of using 21 nodes in the r-direction (Mr = 20), 21 nodes in the λ-direction (Mλ = 20), and 400 time-steps
(38)
The values of ε (N ) experienced in correspondence of different values of N are reported in Table 4. We may note that, in accordance with the theory developed in Young (2008), ε (N ) tends monotonically to zero as N increases. In particular, ε (N ) starts to become of order 10−3 when N is approximately equal to 5000. Therefore, the prices of portfolios comprising more than five thousand life insurances can be accurately approximated by solving the linear partial differential problem (8), (9), (12), (13). This allows one to avoid to pass through the more complex system of non-linear partial differential equations (4), (6), (10), (11). On the contrary, for portfolios with less than five thousand policies, the use of the non-linear partial differential equations (to be solved by an appropriate numerical method such as the one proposed in this paper) seems to be preferable. 5. Conclusions Young’s model for risk insurance is solved using a new numerical method that combines an operator splitting procedure with the differential quadrature finite difference scheme. This approach allows us to reduce the systems of non-linear partial differential problems to systems of linear equations of very small dimension, so that pricing portfolios of many life insurances becomes a relatively easily task. Numerical experiments are presented showing that the method proposed is very accurate and fast. In addition, the limit behavior of portfolios of life insurances
L.V. Ballestra et al. / Insurance: Mathematics and Economics 51 (2012) 442–448
as the number of contracts tends to infinity is investigated. This analysis reveals that portfolios comprising more than five thousand policies can be priced by solving a linear partial differential equation derived in Young (2007, 2008), in place of the more complex system of non-linear partial differential problems. Finally, we point out that the numerical method developed in the present paper could also be applied to several other models in actuarial mathematics, which, from the mathematical standpoint, are substantially analogous to Young’s model. For the reader’s convenience, some of these models are described in Appendix. Appendix The numerical method developed in the present paper could also be used to solve a variety of models in actuarial mathematics, which, from the analytical standpoint, are substantially analogous to Young’s model. Below we briefly describe some of these models, namely the model by Bayraktar et al. (2009), the model by Ludkovski and Young (2008), and the model by Dahl (2004). For the reader’s convenience, we keep as much as possible the same notations used in the papers in which these models have originally been proposed. Therefore, in the following it may happen that the same symbol is used in different models and refers to different quantities, or that the same quantity is used in different models and is referred to by different symbols.
• The model by Bayraktar, Milevsky, Promislow and Young for life annuities In (Bayraktar et al., 2009), Bayraktar et al. consider a portfolio of N (homogeneous) life annuities that pay at a continuous rate of 1 per unit of time until maturity T or until the insured individual dies. Let a(N ) (r , λ, t ) denote the price of this portfolio, as a function of the interest rate r, of the mortality rate λ, and of the time t. According to the model developed in Bayraktar et al. (2009), a(N ) must be obtained solving the following system of partial differential equations:
∂ a(n) 1 ∂ 2 a(n) ∂ a(n) + b(r , t ) + d2 (r , t ) ∂t ∂r 2 ∂r2 ∂ a(n) 1 2 ∂ 2 a(n) + µ(λ, t )(λ − λ) + σ (t )(λ − λ)2 − ra(n) ∂λ 2 ∂λ2 − nλ(a(n) − a(n−1) ) + n (n) 2 2 ∂a + nλ a(n) − a(n−1) , = −α σ 2 (t )(λ − λ)2 ∂λ (r , λ) ∈ [0, +∞) × [λ, +∞), t < T , n = 1, 2, . . . , N ,
∂ H (n) 1 ∂ 2 H (n) ∂ H (n) + (a(t )r + b(t )) + (c (t )r + d(t )) ∂t ∂r 2 ∂r2 2 (n) ∂ H (n) ∂ H 1 2 + µ(λ, t )λ + σ (t )λ2 ∂λ 2 ∂λ2 (n) 2 γ ∂H 1 − rH (n) = − σ 2 (t )λ2 2 F (r , t ) ∂λ F (r , t ) − F (γr ,t ) H (n) −H (n−1) + nλ , 1−e γ (r , λ) ∈ [0, +∞) × [0, +∞), t < T , n = 1, 2, . . . , N , with final condition: H (n) (r , λ, T ) = n,
(r , λ) ∈ [0, +∞) × [0, +∞), n = 1, 2, . . . , N . The function H
(0)
needed in (A.1) to start the recursion is:
a(0) (r , λ, t ) = 0,
(r , λ) ∈ [0, +∞) × [λ, +∞), t ≤ T .
(A.5)
needed in (A.4) to start the recursion is
H (0) (r , λ, t ) = 0,
(r , λ) ∈ Ω , t ≤ T .
(A.6) (N ),a
Instead, the price of the N life annuities H must be obtained solving the following system of partial differential equations:
∂ H (n),a 1 ∂ 2 H (n),a ∂ H (n),a + (a(t )r + b(t )) + (c (t )r + d(t )) ∂t ∂r 2 ∂r2 (n),a 2 (n),a ∂H 1 ∂ H + µ(λ, t )λ + σ 2 (t )λ2 ∂λ 2 ∂λ2 (n),a 2 1 γ ∂H − rH (n),a = − σ 2 (t )λ2 2 F (r , t ) ∂λ γ ( n ), a ( n − F (r , t ) − H −H 1),a 1 − e F (r ,t ) − n, + nλ γ (A.7)
H (n),a (r , λ, T ) = 0, (A.1)
(r , λ) ∈ [0, +∞) × [0, +∞), n = 1, 2, . . . , N .
(A.2)
(A.3)
In (A.1)–(A.3), the quantities b, d, µ, σ , λ and α have exactly the same meaning as in Section 2. From the mathematical standpoint, the system of partial differential equations (A.1)–(A.3) is substantially identical to the system of partial differential equations (4)–(6), and thus could be solved using the HDQ-operator splitting approach developed in the present manuscript. 2
• The model by Ludkovski and Young for pure endowments and life annuities Let r and λ denote the interest rate and the mortality rate, respectively. Moreover, let H (N ) (r , λ, t ) denote the price of N
(A.8)
The function H (0),a needed in (A.7) to start the recursion is H (0) (r , λ, t ) = 0,
(0)
(A.4)
with final condition:
a(n) (r , λ, T ) = 0, The function a
(homogeneous) pure endowments with face value 1 and maturity T , and let H (N ),a (r , λ, t ) denote the value of N (homogeneous) life annuities that pay at a continuous rate of 1 per unit of time until maturity T or until the insured individuals die. In (Ludkovski and Young, 2008), Ludkovski and Young develop a model to compute H (N ) and H (N ),a based on indifference pricing. In particular, H (N ) must be obtained solving the following system of partial differential equations:
(r , λ) ∈ [0, +∞) × [0, +∞), t < T , n = 1, 2, . . . , N ,
with final condition:
(r , λ) ∈ [0, +∞) × [λ, +∞), n = 1, 2, . . . , N .
447
(r , λ) ∈ Ω , t ≤ T .
(A.9)
In (A.4), (A.7), a, b, c, d, σ , µ and F are suitable (known) deterministic functions of their arguments. In particular a, b, c and d model the dynamics of the (stochastic) interest rate, σ and µ model the dynamics of the (stochastic) mortality rate, and F is the market value of a risk-free bond. Finally, in (A.4), (A.7), γ is a constant parameter that models the utility of the insurer. From the mathematical standpoint, the systems of partial differential problems (A.4)–(A.9) are analogous to the system of partial differential problems (4)–(6), and thus could be solved using the HDQ-operator splitting approach developed in the present manuscript.
• The model by Dahl for market reserves
Let J denote a positive integer and let the set J = {1, 2, . . . , J } describe the possible states of a policyholder. As a simple case, it could be J = {1, 2}, with state 1 corresponding to the policyholder being alive, and state 2 corresponding to the policyholder being
448
L.V. Ballestra et al. / Insurance: Mathematics and Economics 51 (2012) 442–448
dead (however other states could be possible as well, for example a state telling if the policyholder is retired). In (Dahl, 2004), Dahl considers a mortality-linked insurance contract to which is associated a generic flow of payments that may depend on the state of the policyholder and on the price of a traded stock. Let s denote the price of the stock, and let µ denote the mortality rate. Finally, let V j (s, µ, t ) denote the market reserve for the insurance contract, which depends also on the state j of the policyholder, j ∈ J . According to the model developed in Dahl (2004), V 1 , V 2 , . . . , V J must satisfy the following system of partial differential equations:
∂ 2V j ∂V j 1 ∂V j + α µ,Q (µ, t ) + (σ µ (µ, t ))2 ∂t ∂µ 2 ∂µ2 2 j j ∂ V ∂V 1 2 + r (s, t )s + σ s (s, t ) s2 2 − r (s, t )V j ∂s 2 ∂s J = −aj (s, t ) − ajk (s, t ) + V k − V j k=1,k̸=j
× 1 + g (µ, t ) Rjk (µ, t ),
jk
(s, µ) ∈ [0, +∞) × [0, +∞), t < T , j = 1, 2, . . . , J ,
(A.10)
with final condition: V j (s, µ, T ) = ∆Aj (s, T ),
(s, µ) ∈ [0, +∞) × [0, +∞), j = 1, 2, . . . , J . µ,Q
µ
(A.11)
In (A.10)–(A.11), α , σ , r, σ , a , a , ∆A , g and Rjk are suitable (known) deterministic functions of their arguments, j = 1, 2, . . . , J, k = 1, 2, . . . , J. In particular α µ,Q and σ µ model the dynamics of the (stochastic) mortality rate, r (s, t ) and σ s model the dynamics of the (stochastic) stock price, aj , ajk , and ∆Aj models the flows of payments associated to the insurance contract, and g jk and Rjk model the transition from state j to state k, j = 1, 2, . . . , J, k = 1, 2, . . . , J. From the mathematical standpoint, the system of partial differential problems (A.10)–(A.11) is very similar to the system of partial differential problems (4)–(6), and thus could be solved using the HDQ-operator splitting approach proposed in this manuscript. Finally, we briefly recall that in (Dahl and Møller, 2006) Dahl and Møller propose a model to compute the market reserve value of a portfolio of pure endowments and derive a (single) partial differential equation analogous to the partial differential equations (A.10) (see Dahl and Møller (2006, p. 210)). Then, the numerical method developed in the present paper could be used to solve such a partial differential equation as well. s
j
jk
j
jk
References Bayraktar, E., Milevsky, M.A., Promislow, S.D., Young, V.R., 2009. Valuation of mortality risk via the instantaneous Sharpe ratio: applications to life annuities. Journal of Economic Dynamics and Control 33 (676–691). Bayraktar, E., Young, V.R., 2007. Hedging life insurance with pure endowments. Insurance: Mathematics and Economics 40, 435–444. Chang, C.-C., Chung, S.-L., Stapleton, R.C., 2007. Richardson extrapolation techniques for pricing American-style options. Journal of Futures Markets 27, 791–817. Chen, W., Wang, X., Zhong, T., 1996. The structure of weighting coefficient matrices of harmonic differential quadrature and its application. Communications in Numerical Methods in Engineering 12, 455–460. Chiarella, C., Kang, B., Meyer, G.H., Ziogas, A., 2009. The evaluation of American option prices under stochastic volatility and jump-diffusion dynamics using the method of lines. International Journal of Theoretical and Applied Finance 12, 393–425. Cox, J.C., Ingersoll, J.E., Ross, S.A., 1985. A theory of the term structure of interest rates. Econometrica 53, 385–408. Dahl, M., 2004. Stochastic mortality in life insurance: market reserves and mortality-linked insurance contracts. Insurance: Mathematics and Economics 35 (113–136). Dahl, M., Møller, T., 2006. Valuation and hedging of life insurance liabilities with systematic mortality risk. Insurance: Mathematics and Economics 39 (193–217). Janghorban, M., 2011. Two different types of differential quadrature methods for static analysis of microbeams based on nonlocal thermal elasticity theory in thermal environment. Archive of Applied Mechanics 1–7. Ludkovski, M., Young, V.R., 2008. Indifference pricing of pure endowments and life annuities under stochastic hazard and interest rates. Insurance: Mathematics and Economics 42 (14–30). Malekzadeh, P., Karami, G., 2005. Polynomial and harmonic differential quadrature methods for free vibration of variable thickness thick skew plates. Journal of Engineering Structures 27, 1563–1574. Milevsky, M.A., Promislow, S.D., Young, V.R., 2005. Financial valuation of mortality risk via the instantaneous Sharpe ratio: applications to pricing pure endowments, Technical Report. Available at http://arxiv.org/abs/0705.1302 (accessed May 2012). Quarteroni, A., Sacco, R., Saleri, F., 2000. Numerical Mathematics. Springer-Verlag, New-York. Shu, C., 2000. Differential Quadrature and Its Application in Engineering. SpringerVerlag, London. Shu, C., Richard, B.E., 1992. Application of generalized differential quadrature to solve two-dimensional incompressible Navier–Stokes stress equations. International Journal for Numerical Methods in Fluids 15, 791–798. Shu, C., Xue, H., 1997. Explicit computation of weighting coefficients in the harmonic differential quadrature. Journal of Sound and Vibration 204, 549–555. Shu, C., Yao, Q., Yeom, K.S., Zhu, Y.D., 2002. Numerical analysis of flow and thermal fields in arbitrary eccentric annulus by differential quadrature method. Heat and Mass Transfer 38, 597–608. Striz, A.G., Chen, W.L., Bert, C.W., 1997. Free vibration of plates by the high accuracy quadrature element method. Journal of Sound and Vibration 202, 689–702. Striz, A.G., Wang, X., Bert, C.W., 1995. Harmonic differential quadrature method and applications to analysis of structural components. ACTA Mechanica 111, 85–94. Young, V.R., 2007. Pricing life insurance under stochastic mortality via the instantaneous Sharpe ratio: theorems and proofs, Technical Report. Available at http://arxiv.org/abs/0705.1297 (accessed May 2012). Young, V.R., 2008. Pricing life insurance under stochastic mortality via the instantaneous Sharpe ratio. Insurance: Mathematics and Economics 42 (691–703).