The Thermodynamic Approach to Whole-Life Insurance: A Method for Evaluation of Surrender Risk
aa r X i v : . [ q -f i n . GN ] D ec The Thermodynamic Approach to Whole-LifeInsurance: A Method for Evaluation ofSurrender Risk
Jirˆo Akahori, Yuuki Ida, Maho Nishida, and Shuji TamadaDepartment of Mathematical Sciences, Ritsumeikan UniversityDecember 2020
Abstract
We introduce a collective model for life insurance where the het-erogeneity of each insured, including the health state, is modeled bya diffusion process. This model is influenced by concepts in statisti-cal mechanics. Using the proposed framework, one can describe thetotal pay-off as a functional of the diffusion process, which can beused to derive a level premium that evaluates the risk of lapses due tothe so-called adverse selection. Two numerically tractable models arepresented to exemplify the flexibility of the proposed framework.
Keywords:
Life insurance, Surrender Risk, Collective model, Feynmann-Kac Formula
The risk of lapses, also referred to as surrender risk, is the instabilityassociated with unexpected lapses of insurance contracts, which mayresult in a huge loss for the insurer. As insurance contracts are collec-tive in nature, we can divide the cause of lapses into two categories:homogeneous and heterogeneous (among the policy holders). The for-mer is basically due to macro-economic shocks — typically changes ininterest rates, recessions, inflation, and so forth. The surrender riskarising from interest rate fluctuations, for example, can be evaluatedusing option-pricing type technologies (see e.g. [3]). he latter, i.e. heterogeneous causes, can be further divided intotwo groups: economic and non-economic. A policyholder may surren-der because she is unexpectedly short of money while the economyis, as a whole, healthy. Such a cause is classified as heterogeneous-economic (see e.g. [13]).Among non-economic causes, demographic heterogeneity, whichis variation in the force of mortality, has been a central issue bothacademically and practically, as it has been recognized as a (potential)cause of the so-called moral hazard or adverse selection. The riskfrom heterogeneity has been recognized in insurance (since the 19-th century!), as is pointed out in the seminal paper by G.A. Akerlof[2], which made clear the role of the “asymmetry of information” inadverse selection. M. Rothschild and Nobel laureate J. Stiglitz provedthe non-existence of an economic equilibrium under asymmetry in [16].According to that paper, in a “rational world”, adverse selection leadsto the non-existence of an equilibrium, implying potential instabilityof a life insurance contract.Since then there have been many studies on life insurance, howevermost empirical studies do not seem to observe the predicted negativeeffects of adverse selection (see e.g. [6]). Some recent studies, like[12], suggest the effect of so-called advantageous selection: health-ier people might be more risk averse. Thus many theoretical studieshave been rather interested in modelling the effect of heterogeneity onthe surrender risk without the hypothesis of rationality. One of themain streams can be found in direct modelling of the dependence ofdemographic heterogeneity and surrender inclination. In most cases,however, this is modeled by a simple stochastic model, a discrete-time,discrete-state Markov chain at best; see e.g. [5], [8], [18], and morerecently [1], to name a few. By contrast, the present paper proposesa model using diffusion processes.We will start in section 2.1 with introducing the basic framework ofour model, without mortality or surrender risk, where the heterogene-ity of each insured person is modeled by a multi-dimensional diffusionprocess aiming to describe the various causes discussed above. Theintroduced diffusion processes naturally define a probability measure,which describes the distribution of the “heterogeneity”, and we assumethat it is approximated by its infinite-agent limit, which is analogous toa thermodynamic limit (Assumption 2.2 and Theorem 2.3). The“thermodynamic” procedure is the core of our framework. Then insection 2.2, we introduce the “lifetime” of an insured, modeled bythe killing time of an associated diffusion process (Assumption 2.4).The thermodynamic probability measure is then compensated by thekilling rate (Theorem 2.7). The basic framework with the lifetime is sed in section 2.3 to model the cash-flow of a whole-life insurancewith a level premium. Its continuous-time version is presented in sec-tion 2.4, where, after taking the “thermodynamic limit”, the premiumcan be obtained by calculating the Laplace transform of expectationswith respect to the heterogeneity distribution (see equation (16) inTheorem 2.9). The Laplace transform method is another core of ourframework and enables computational tractability.Then, in section 3, we add to the introduced model the surrendertime, which is again modelled by a killing time (Assumption 3.1). Un-der the continuous-time thermodynamic model, the pay-off, which isthe revenue minus expenditure (of the insurance company), is mod-elled by the Laplace transform of (27), after deriving the thermody-namic limit under the extended setting in Theorem 3.2.In sections 4 and 5, we introduce two specific models where we cancalculate the formula (27) analytically. In the former, we rely on thelocal constancy of the killing rates, which can however approximatefairly general rates. If the approximated rate is highly non-linear, thenwe need to resort to the numerical algorithm proposed in Theorem4.2, which is among the mathematical contributions of the presentpaper. In contrast, the latter model gives us completely an analyticalexpression of (27) using the symmetry of the 2-dimensional Besselprocess.There are some studies analysing the surrender risk in the spiritof quantitative finance, such as [10] and more recently [4], which usecontinuous time processes similar to ours; however, these are not con-cerned with heterogeneity. They model the lapses by “jumps”, that is,exogenous events. Such an approach might be called a “reduced formapproach”. From this point of view, our model can be understoodas a structural version of, for example, the model proposed by O. LeCourtois and H. Nakagawa [10] (see Remark 3.4 in Section 3 below). Acknowledgments
The authors would like to thank Corina Constantinescu (Univ. Liv-erpool) and Gregory Markowsky (Monash Univ.) for review and forvaluable comments and suggestions. Model without Surrender Risk
Let X t be a time-homogeneous diffusion process in R d ,( { X t : t ∈ R + } , { P x : x ∈ R n } ) , P x ( X = x ) = 1 . (1)Here, we consider X to be a quantified personal profile such as health,economic state, or other conditions of a person, which depend on atime parameter t ∈ R + . The infinitesimal generator of X will bedenoted by L . Assumption 2.1.
We assume that the transition probability of X has a smooth density; that is, there exists a smooth function q suchthat P x ( X t ∈ A ) = Z A q ( t, x, y ) dy, ( A ∈ B ( R d )) . (2)We further assume that P x ( X t ∈ A ) is a continuous function in x forany fixed t > A ∈ B ( R d ).Our model for life insurance is as follows. Let N ∈ N and I = { , · · · , N } be the number and the set of the initial insured, respec-tively. The initial personal profile of each insured is x i ∈ R d for i ∈ I .We assume that each of X i ( i = 1 , · · · , N ) is distributed as P x i . Wealso assume that the initial profile of each insured takes values only in N − Z d ≡ { k/N : k ∈ Z d } We define a probability measure µ N on R d whose support is N − Z d by µ N ( A ) := 1 N ♯ { i ∈ I : x i ∈ A } for A ∈ B ( R d ) . (3)We call it the initial distribution measure . Assumption 2.2.
We assume that there exists f ≥ k f k = 1such that lim N →∞ Z R d h (cid:16) xN (cid:17) µ N ( dx ) = Z R d h ( x ) f ( x ) dx (4)for any bounded continuous function h .Let a random counting measure v N be defined by v N ( t, A ) := 1 N ♯ { i ∈ I : X it ∈ A } , ( A ∈ B ( R d )) . This expresses the proportion of the insured whose profile is in A attime t . heorem 2.3. Under
Assumptions 2.1 and , we have that lim N →∞ E [ v N ( t, A )] = Z A E [ f ( X ∗ t ) | X ∗ = x ] dx, where X ∗ is the adjoint process of X , that is, a diffusion process whoseinfinitesimal generator is L ∗ , the adjoint operator of L with respect to L ( R d , Leb) .Proof.
Observe that E [ v N ( t, A )] = E X j ∈I N { i ∈I : X it ∈ A } ( j ) = 1 N X i ∈I P ( X it ∈ A ) = 1 N X i ∈I P x i ( X t ∈ A )= 1 N X k ∈ Z d P kN ( X t ∈ A ) µ N ( { k/N } ) N = X k ∈ Z d P kN ( X t ∈ A ) µ N ( { k/N } ) . Letting N → ∞ , we getlim N →∞ E [ v N ( t, A )] = Z R d P x ( X t ∈ A ) f ( x ) dx by using (4) in Assumption 2.2. Moreover, this formula can be rewrit-ten (where we can use Fubini’s theorem since both P x ( X t ∈ A ) and f ( x ) are positive functions) aslim N →∞ E [ v N ( t, A )] = Z P y ( X t ∈ A ) f ( y ) dy = Z (cid:18)Z q ( t, y, x )1 A ( x ) dx (cid:19) f ( y ) dy = Z A ( x ) dx Z q ( t, y, x ) f ( y ) dy. Therefore, A u ( t, A ) is an absolutely continuous measure with re-spect to Lebesgue measure, and its density is given by u ( t, x ) := Z q ( t, y, x ) f ( y ) dy. Using the adjoint process (with respect to Lebesgue measure) X ∗ , thefunction u can be expressed as u ( t, x ) = E [ f ( X ∗ t ) | X ∗ = x ] =: E x [ f ( X ∗ t )] , (5)where we abuse the notation E x a bit. .2 Mortality Model We introduce a mortality model, where the force of mortality is depen-dent only on the current personal profile. In the model, the insureddo not surrender the policy.Let ˆ X be a killed process obtained from the Markov process (1)defined with a random time ζ as follows,ˆ X t = X t { ζ>t } + ∞ { ζ ≤ t } . (6) Assumption 2.4.
We assume that P ( ζ > t | σ ( X s : s ≤ t )) = e − R t V ( X s ) ds (7)and lim t →∞ P ( ζ > t | σ ( X s : s ≤ t )) = 0 (8)with a positive function V . Remark 2.5.
The mortality in our model is consistent with the ones[17], and [19], [20] in demography.Then, we know that (see e.g. [15, chapter III, Theorem 18.6]) thetransition semigroup ( ˆ T t ) t ≥ of ˆ X is given byˆ T t f ( x ) := E [ f ( ˆ X s ) | X = x ]= E [ f ( X s ) e − R t V ( X s ) ds | X = x ] , ( f ∈ C ( R d )) . By the Feynman-Kac formula (see e.g. [15, chapter III 19.]), theinfinitesimal generator ˆ L of the semigroup ˆ T is given, using the in-finitesimal generator L of T , as followsˆ L f ( x ) = L f ( x ) − V ( x ) f ( x ) , ( f ∈ D ( L )) , where the domain of L is, as usual, the space of such f thatlim t → ˆ T t f − ft exists in C ( R d ). We will sometimes use the notationˆ T t f ( x ) = E x [ f ( ˆ X )] = E x [ f ( X ) e − R t V ( X s ) ds ]for the purpose of clarifying both the starting point and the samplepath which we are looking at, even though this may again be a bit ofan abuse of notation. ssumption 2.6. We assume that ˆ T t has a density, that is, thereexists a smooth function q V such thatˆ T t f ( x ) = Z R d q V ( t, x, y ) f ( y ) dy. Let the model size be N as in the preceding section, and the initialdistribution measure µ N be as (3). Let a random counting measure v N be redefined by v N ( t, A ) = 1 N ♯ { i ∈ I : ˆ X it ∈ A } , ( A ∈ B ( R d )) . Theorem 2.7.
Let A ∈ B ( R d ) . Under Assumptions 2.6 and , wehave that lim N →∞ E [ v N ( t, A )] = Z A E [ f ( X ∗ t ) e − R t V ( X ∗ s ) ds | X ∗ = x ] dx. (9) Proof.
Let us calculate E [ v N ( t, A )] as E [ v N ( t, A )] = E X j ∈I N { i ∈I : ˆ X it ∈ A } ( j ) = 1 N X i ∈I P ( ˆ X it ∈ A ) = 1 N X i ∈I P x i ( ˆ X t ∈ A )= X k ∈ Z d P k/N ( ˆ X t ∈ A ) µ N ( { k/N } ) . Letting N → ∞ , we get the formulalim N →∞ E [ v N ( t, A )] = Z R d P y ( ˆ X t ∈ A )( y ) f ( y ) dy by using (4) in Assumption 2.2. Then, by the same procedure as weused in the proof of Theorem 2.3, we get (9). We consider an insurance where both the premium and the insuranceproceeds are paid discretely at time t = 0 , , · · · . Let p and A be the(level) premium and the sum insured, respectively. Then, the revenueof the insurance company v t ( p ) at each time t is given by v t ( p ) = p · ♯ { i ∈ I : ˆ X it = ∞} ( t = 0 , , , ... ) . n our model, the expenditure of the insurance company c t at eachtime t is given by c = 0 ,c t = A · ♯ { i ∈ I : ˆ X it = ∞ , ˆ X it − = ∞} ( t = 1 , , ... ) . Let the expected return R d of the discrete-time insurance model attime 0 be defined by R d ( N, p ) := ∞ X t =0 e − rt E [ v t ( p ) − c t ] , (10)and p d ( N ) be the solution of R d ( N, p ; N ) = 0. Note that the solutionis unique and strictly positive since (10) is linear in p and we haveclearly ∞ X t =0 e − rt P ( i ∈ I : ˆ X it = ∞ ) > A ∞ X t =0 e − rt E [ ♯ { i ∈ I : ˆ X it = ∞ , ˆ X it − = ∞} ] > Theorem 2.8.
We have that p d ( N ) = A ∞ X t =1 e − rt X k ∈ Z d E kN [ e − R t − V ( X s ) ds − e − R t V ( X s ) ds ] µ N ( { k/N } ) ∞ X t =0 e − rt X k ∈ Z d E kN [ e − R t V ( X s ) ds ] µ N ( { k/N } ) . Moreover, p d ( ∞ ) defined by p d ( ∞ ) := lim N →∞ p d ( N ) is given by p d ( ∞ ) = A ∞ X t =1 e − rt Z R d E y [ e − R t − ˜ V ( X ∗ s ) ds f ( X ∗ t − ) − e − R t ˜ V ( X ∗ s ) ds f ( X ∗ t )] dy ∞ X t =0 e − rt Z R d E y [ e − R t ˜ V ( X ∗ s ) ds f ( X ∗ t )] dy , using the adjoint process. roof. We first obtain the expected return R d ( N, p ) at time 0. Theexpected revenue at time 0 is ∞ X t =0 e − rt E [ v t ( p )] = ∞ X t =0 e − rt E h p · { i ∈I : ˆ X it = ∞} i = p ∞ X t =0 e − rt X i ∈I P ( ˆ X it = ∞ )= p ∞ X t =0 e − rt X k ∈ Z d P kN ( ˆ X t = ∞ ) µ N ( { k/N } ) N. Since { ˆ X t = ∞} = { ζ > t } , the expected revenue is now given by= p ∞ X t =0 e − rt X k ∈ Z d P kN ( ζ > t ) µ N ( { k/N } ) N = p ∞ X t =0 e − rt X k ∈ Z d E kN [ E kN [1 { ζ>t } | σ ( X s : s ≤ t )]] µ N ( { k/N } ) N = p ∞ X t =0 e − rt X k ∈ Z d E kN [ e − R t V ( X s ) ds ] µ N ( { k/N } ) N. (11)Next, we calculate the expected expenditure at time 0 as E [ c t ] = A E h { i ∈I : ˆ X it = ∞ , ˆ X it − = ∞} i = A X i ∈I P ( ˆ X it = ∞ , ˆ X it − = ∞ )= A X k ∈ Z d P kN ( ˆ X t = ∞ , ˆ X t − = ∞ ) µ N ( { k/N } ) N. As above we can use expressions with ζ instead, and the expectedexpenditure is now= A X k ∈ Z d P kN ( t − < ζ ≤ t ) µ N ( { k/N } ) N = A X k ∈ Z d (cid:16) E kN [ E kN [1 { ζ>t − } − { ζ>t } | σ ( X s : s ≤ t )]] (cid:17) µ N ( { k/N } ) N = A X k ∈ Z d E kN [ e − R t − V ( X s ) ds − e − R t V ( X s ) ds ] µ N ( { k/N } ) N. (12)Using the formulas (11) and (12), the expected return R d ( N, p ) is alculated as follows: R d ( N, p )= ∞ X t =0 e − rt E [ v t ( p ) − c t ]= p ∞ X t =0 e − rt X k ∈ Z d E kN [ e − R t V ( X s ) ds ] µ N ( { k/N } ) N − A ∞ X t =1 e − rt X k ∈ Z d E kN [ e − R t − V ( X s ) ds − e − R t V ( X s ) ds ] µ N ( { k/N } ) N. (13)Thus, the premium p d ( N ) is given by p d ( N ) = A ∞ X t =1 e − rt X k ∈ Z d E kN [ e − R t − V ( X s ) ds − e − R t V ( X s ) ds ] µ N ( { k/N } ) ∞ X t =0 e − rt X k ∈ Z d E kN [ e − R t V ( X s ) ds ] µ N ( { k/N } ) . (14)Letting N → ∞ , which is possible since each term in (14) convergesby Assumption 2.2, we get the following formula: p d ( ∞ ) = A ∞ X t =1 e − rt Z R d E x [ e − R t − V ( X s ) ds − e − R t V ( X s ) ds ] f ( x ) dx ∞ X t =0 e − rt Z R d E x [ e − R t V ( X s ) ds ] f ( x ) dx . Moreover, this formula can be rewritten by using the adjoint processand the potential ˜ V of L ∗ as follows: p d ( ∞ ) = A ∞ X t =1 e − rt Z R d E y [ e − R t − ˜ V ( X ∗ s ) ds f ( X ∗ t − ) − e − R t ˜ V ( X ∗ s ) ds f ( X ∗ t )] dy ∞ X t =0 e − rt Z R d E y [ e − R t ˜ V ( X ∗ s ) ds f ( X ∗ t )] dy . .4 Continuous-Time Level-Premium InsuranceModel We will consider a continuous-time payment model in this section.The revenue of the insurance company during [0 , t ] is given by Z t e − rs v s ( p ) ds, where v t ( p ) = p · ♯ { i ∈ I : ˆ X it = ∞} , while the expenditure of the insurance company during [0 , t ] is givenby C t = X j ∈I Ae − rζ j { i ∈I : ζ i Theorem 2.9. We have that p c ( N ) = A Z ∞ e − rt X k ∈ Z d E kN [ V ( X s ) e − R t V ( X s ) ds ] µ N ( { k/N } ) dt Z ∞ e − rt X k ∈ Z d E kN [ e − R t V ( X s ) ds ] µ N ( { k/N } ) dt . (15) Moreover, p c ( ∞ ) defined by p c ( ∞ ) := lim N →∞ p c ( N ) is given by p c ( ∞ ) = A Z ∞ e − rt Z R d E y [ f ( X ∗ t ) ˜ V ( X ∗ t ) e − R t ˜ V ( X ∗ s ) ds ] dydt Z ∞ e − rt Z R d E y [ f ( X ∗ t ) e − R t ˜ V ( X ∗ s ) ds ] dydt , (16) where ˜ V is the potential of L ∗ . roof. We first obtain the expected return R c ( N, p ). Calculating v t ( p )like we did in (11), we get Z ∞ e − rt E [ v t ( p )] dt = p Z ∞ e − rt X k ∈ Z d E kN [ e − R t V ( X s ) ds ] µ N ( { k/N } ) N dt, (17)while the expected expenditure at time 0:, E [ C t ] = E X j ∈I Ae − rζ j { i ∈I : ζ i In this section we consider a model where insurers can surrender theirpolicy by balancing their personal conditions with the premium.In addition to ζ , we introduce a new random time ξ ( p ), which isdependent on a parameter p , satisfying the following assumption. Assumption 3.1. We assume that1. P ( ξ ( p ) > t | σ ( X s : s ≤ t )) = e − R t D ( p,X s ) ds with a positive measurable function D : [0 , ∞ ) × R d ∋ ( p, x ) D ( p, x ) ∈ R > , which is increasing in p and decreasing in x .2. The random times ζ and ξ ( p ) are conditionally independent inthe following sense: P ( ξ ( p ) > t, ζ > t | σ ( X s : s ≤ t ))= P ( ξ ( p ) > t | σ ( X s : s ≤ t )) P ( ζ > t | σ ( X s : s ≤ t ))= e − R t V ( X s ) ds − R t D ( X s ,p ) ds . (21) The revenue of the insurance company during [0 , t ] is given by Z t e − rs v s ( p ) ds, where v t ( p ) = p · ♯ { i ∈ I : ζ i > t, ξ i ( p ) > t } , (22)while the expenditure of the insurance company during [0 , t ] is givenby C t = A X j ∈I e − rζ j { i ∈I : ζ i ≤ t, ξ i ( p ) >t } ( j ) . (23)The expected return at time 0 is the same as the one in the previoussection, namely, R s ( N, p ) = Z ∞ e − rt E [ v t ( p )] dt − E [ C ∞ ] , where the subscript s is put to indicate that it is the one with surrenderrisk. It should be noted that the expected return is no longer a linearfunction in p , and thus we may not have uniqueness of the solution p for R c ( N, p ) = 0. heorem 3.2. We have that R s ( N, p )= Z ∞ e − rt E [ v t ( p )] dt − E [ C ∞ ]= Z ∞ e − rt X k ∈ Z d E kN [( p − AV ( X t )) e − R t V ( X s )+ D ( X s ,p ) ds ] µ N ( { k/N } ) N dt. (24) Proof. The proof is almost the same as for Theorem 2.9. First, Z ∞ e − rt E [ v t ( p )] dt = Z ∞ e − rt E (cid:2) p · ♯ { i ∈ I : ζ i > t, ξ i ( p ) > t } (cid:3) dt = p Z ∞ e − rt X i ∈I P ( ζ i > t, ξ i ( p ) > t ) dt = p Z ∞ e − rt X k ∈ Z d P kN ( ζ > t, ξ ( p ) > t ) µ N ( { k/N } ) N dt = p Z ∞ e − rt X k ∈ Z d E kN [ E [1 { ζ>t,ξ ( p ) >t } | σ ( X s : s ≤ t )]] µ N ( { k/N } ) N dt This formula can be rewritten, by (21), as= p Z ∞ e − rt X k ∈ Z d E kN [ e − R t ( V ( X s )+ D ( X s ,p )) ds ] µ N ( { k/N } ) N dt. (25)Next, the expected expenditure at time 0, E [ C ∞ ] = lim t →∞ E A X j ∈I e − rζ j { i ∈I : ζ i Since R s ( N, < R s is continuous in p , thesolution p ∗ of R S ( N, p ∗ ) = 0 exists if and only if R s ( N, p ) > p . An evident sufficient condition that ensures the exis-tence of the solution is that V and D are bounded, since this implieslim p →∞ R s ( N, p ) = + ∞ . Remark 3.4. In the paper by O. Le Courtois and H. Nakagawa[10], the number of surrenders is modeled by a counting process,and the expected number of the remaining participants at t is ex-pressed by its intensity process, while in our framework it is given by R f ( x ) E [ e − R t D ( X s ,p ) ds ] dx using the “state process” X . In this section, we specifically assume that the personal conditionprocess X t is one dimensional Brownian motion with drift; that is, X t = aW t + bt, here W t is standard Brownian motion, a > , b ∈ R , so that P x ( X t ∈ A ) = Z A √ πa t e − ( y − x − bt )22 a t dy, A ∈ B ( R ) . Clearly, Assumption 2.1 is satisfied.Moreover, we assume that the killing rate functions V ( y ) and D ( y, p ) are step functions; that is, V ( y ) = M X i =1 λ i ( y i − ,y i ] ( y ) ,D ( y, p ) = M X i =1 µ i ( p )1 ( y i − ,y i ] ( y ) , ( M ∈ N , λ i , µ i ( p ) ∈ R , −∞ = y < y < ... < y M = ∞ ) . The model is so designed that we can fit/calibrate any data to a certainextent provided that the personal condition is expressed by a realnumber.In this case, the expected return R s ( N, p ) can be calculated, butonly if we consider the inversion of a large scale matrix to be tractable.Below we first show how the expected return is calculated out of aninversion of a large matrix. Then, we propose a numerical scheme toreduce the burden.We define z V and z as follows: for y ∈ R , z V ( y ) := Z ∞ e − rt E y [ V ( X t ) e − R t ( V ( X s )+ D ( X s ,p )) ds ] dt, and z ( y ) := Z ∞ e − rt E y [ e − R t ( V ( X s )+ D ( X s ,p )) ds ] dt. Then we can express the expected return as R s ( N, p ) = ∞ X k = −∞ N µ N (cid:18) kN (cid:19) (cid:18) pu (cid:18) kN (cid:19) − Au V (cid:18) kN (cid:19)(cid:19) . Lemma 4.1. We have that z V ( y ) = M X i =1 n C i, + e α i, + y + C i, − e α i, − y + γ i o { y i − Using the Feynman-Kac formula, we have that z V and z sat-isfy − V ( y ) + rz V ( y ) = 12 a z ′′ V ( y ) + bz ′ V ( y ) − ( V ( y ) + D ( y, p )) z V ( y ) (30) − rz ( y ) = 12 a z ′′ ( y ) + bz ′ ( y ) − ( V ( y ) + D ( y, p )) z ( y ) (31) or y ∈ R \ { y , · · · , y M − } , respectively. Since z , z V ∈ C (see e.g.[9, Theorem 6.4.1]), z ∗ ( y i +) = z ∗ ( y i − ) , z ′∗ ( y i +) = z ′∗ ( y i − ) , ( i = 1 , , ..., M − ∗ = V, 1. Specifically, C i, + e α i, + y i + C i, − e α i, − y i + γ i = C i +1 , + e α i +1 , + y i + C i +1 , − e α i +1 , − y i + γ i +1 (32)and C i, + α i, + e α i, + y i + C i, − α i, − e α i, − y i = C i +1 , + α i +1 , + e α i +1 , + y i + C i +1 , − α i +1 , − e α i +1 , − y i (33)for i = 1 , , ..., M − 1, which is equivalent to the equation (28). Simi-larly the equations derived for e C i, ± , which is obtained by replacing γ i with e γ i , is equivalent to (29). As we remarked already, the inversion of the matrix − e α + − y + e − α −− y t J e α − + y − e α ++ y J ( − e α + − y + e − α −− y t J ) α + − ( e α − + y − e α ++ y J ) α − + ! might become too heavy if M is large. We then propose a numericalscheme to solve the equations (28) and (29) which might work wellwhen δ := min ≤ i ≤ M − ( y i +1 − y i )is very small.For simplicity, we assume that y , y , · · · , y M − are equally spaced: δ = y i +1 − y i , i = 1 , , · · · , M − , where we can get an exact expression. First let us observe that − e α + − y + e − α −− y t J e α − + y − e α ++ y J ( − e α + − y + e − α −− y t J ) α + − ( e α − + y − e α ++ y J ) α − + ! = I − t J e α + − δ I − J e α − + δ ( I − t J e α + − δ ) α + − I − J e α − + δ α − + ! − e α + − y R M − ⊗ R M − R M − ⊗ R M − e α − + y ! =: L ( δ ) − e α + − y R M − ⊗ R M − R M − ⊗ R M − e α − + y ! . ince we have immediately − e α + − y R M − ⊗ R M − R M − ⊗ R M − e α − + y ! − = − e − α + − y R M − ⊗ R M − R M − ⊗ R M − e − α − + y ! , we can concentrate on the inversion of L ( δ ).Noting that L ( δ ) = (cid:18) I − t J I − J ( I − t J ) α + − ( I − J ) α − + (cid:19) + X k =1 ∞ δ k k ! (cid:18) t J ( α + − ) k J ( α − + ) kt J ( α + − ) k +1 J ( α − k ) k +1 (cid:19) , we propose the following scheme to invert the matrix L ( δ ), where weactually will need the inversion of L (0) = (cid:18) I − t J I − J ( I − t J ) α + − ( I − J ) α − + (cid:19) . For convenience, we denote (cid:18) t J ( α + − ) k J ( α − + ) kt J ( α + − ) k +1 J ( α − k ) k +1 , (cid:19) =: L ( k ) (0) . Theorem 4.2. (i) Let c ∈ R M − be fixed. Define x k , k = 0 , , · · · recursively by x k = ( L (0) − c k = 0 , − L (0) − P k − j =0 1( k − j )! L ( k − j ) (0) x j k ≥ . Then x := P ∞ k =0 δ k x k satisfies L ( δ ) x = c . (ii) In particular, x − P nk =0 x k = O ( δ n +1 ) for each n . (iii) We have explicitly L (0) − = (cid:18) A α − + ( I − J ) − A ( I − J ) − A α + − ( I − t J ) − A ( I − t J ) − (cid:19) where A = (cid:18) ( α M, − − α , + ) − ( ˆ α − − ˆ α + ) − ( ˆ α + I − α M, − I )1 R M − ( ˆ α − − ˆ α + ) − ( α , + − α M, − ) − t R M − . (cid:19) and A = (cid:18) t R M − ( α M, − − α , + ) − ( ˆ α + − ˆ α − ) − ( α M, − − α , + ) − ( ˆ α + − ˆ α − ) − ( ˆ α − I − α , + I )1 R M − (cid:19) , with ˆ α ± := diag( α , ± , · · · , α M − , ± ) ∈ R M − ⊗ R M − . roof. (i), (ii) For each n ∈ N define L n := n X k =0 δ k k ! L ( k ) (0) . Then, L − L n = O ( δ n +1 ) . Since L n n X k =0 x k = O ( δ n +1 ) , by a standard argument we have the assertion (ii), and hence (i).(iii) Put K := ( I − J ) − ( I − t J )= (cid:18) t R M − − I R M − ⊗ R M − R M − (cid:19) . Then we have that L (0) − = (cid:18) I − t J I − J ( I − t J ) α + − ( I − J ) α − + (cid:19) − = (cid:18) ( α − + K − Kα + − ) − α − + ( I − J ) − ( Kα + − − α − + K ) − ( I − J ) − ( α + − K − − K − α − + ) − α + − ( I − t J ) − ( K − α − + − α + − K − ) − ( I − t J ) − (cid:19) . We then see that α + K − Kα − = (cid:18) α , + t R M − R M − ˆ α + I R M − ⊗ R M − (cid:19) (cid:18) t R M − − I R M − ⊗ R M − R M − (cid:19) − (cid:18) t R M − − I R M − ⊗ R M − R M − (cid:19) (cid:18) ˆ α − I R M − ⊗ R M − R M − t R M − α M, − (cid:19) = (cid:18) t R M − α , + − α M, − ˆ α − − ˆ α + ( ˆ α + I − α M, − I )1 R M − (cid:19) = A − and α + − K − − K − α − + = (cid:18) ( ˆ α − I − α , + I )1 R M − ˆ α + − ˆ α − α M, − − α , + t R M − (cid:19) = A − . A Model by the 2-Dimensional SquaredBessel Process In this section, we specifically assume that the personal conditionprocess X t is the 2-dimensional squared Bessel Process dX t = 2 p X t dW t + 2 dt ( a > . (34)We further assume that the initial condition distribution is approx-imated by the exponential distribution whose mean is γ ; the limitdensity f in Assumption 2.2 is given by f ( x ) = γe − γx ( γ > . Moreover, we assume that the killing rate functions are V ( y ) and D ( y, p ) V ( x ) = mx + n,D ( x, p ) = ϕ ( p ) x + ̺ ( p ) , ( m > , ϕ ( p ) < , n, ̺ ( p ) ∈ R ) . We call this the 2-dimensional Squared Bessel model, 2SB model forshort. The 2SB model, by nature, satisfies Assumptions 2.1, 2.2, and3.1. Theorem 5.1. The virtual average expected return in the 2SB modelis explicitly calculated as: VAR s ( p ) = Amλ + (cid:18) γp + Amcλ − γAn (cid:19) γ √ λ + 2 λ × (cid:18) F , c √ λ − , c √ λ − − γ − √ λγ + √ λ ! − c + √ λ (cid:19) , where F is the hypergeometric function, c := r + n + ̺ ( p ) , and λ := m + ϕ ( p ) . Proof. First note thatVAR s ( p ) = Z ∞ e − rt dt Z ∞ γe − γx dx × E [( p − A ( mX t + n )) e − ( m + ϕ ( p )) R t X s ds ] | X = x ] e − ( n + ̺ ( p )) t . ow we see that we need to calculate γ Z ∞ e − ct Z ∞ e − γx E [( aX t + b ) e − λ R t X s ds | X = x ] dx dt for λ = m + ϕ ( p ) ,c = r + n + ϕ ( p ) ,a = − Am, and b = p − An, which can be further reduced to the calculation of I γ ( c, λ ) := Z ∞ e − ct Z ∞ e − γx E [ e − λ R t X s ds | X = x ] dx dt since γ Z ∞ e − ct Z ∞ e − γx E [ X t e − λ R t X s ,ds | X = x ] dx dt = − γλ Z ∞ e − ct Z ∞ e − γx ∂ t E [ e − λ R t X s ds | X = x ] dx dt = − λ − cλ Z ∞ e − ct Z ∞ e − γx E [ e − λ R t X s ds | X = x ] dx dt. That is,VAR s ( p ) = Amm + ϕ ( p )+ (cid:18) γp + Am ( r + n + ̺ ( p )) m + ϕ ( p ) − γAn (cid:19) I γ ( r + n + ̺ ( p ) , m + ϕ ( p )) . It is well-known that (see e.g. Ikeda-Watanabe [7]) that E [ e − λ R t X s ds | X = x ] = e − x √ λ tanh √ λt cosh √ λt , herefore we have I = Z ∞ e − ct ( γ + √ λ tanh √ λt ) cosh √ λt dt = Z ∞ e − ct ( γ cosh √ λt + √ λ sinh √ λt ) dt = Z ∞ γ + √ λ e − ( c + √ λ ) t γ −√ λγ + √ λ e − √ λt dt = Z ∞ e − ( c + √ λ ) t γ + √ λ ∞ X j =1 − γ − √ λγ + √ λ e − t ! j dt = 2 γ + √ λ ∞ X j =1 − γ − √ λγ + √ λ ! j Z ∞ e − ( c +(2 j +1) √ λ ) t dt = 2 γ + √ λ ∞ X j =1 − γ − √ λγ + √ λ ! j c + (2 j + 1) √ λ = 1 γ √ λ + 2 λ ∞ X j =1 − γ − √ λγ + √ λ ! j (cid:16) c √ λ − (cid:17) j (cid:16) c √ λ − (cid:17) j = 1 γ √ λ + 2 λ × (cid:18) F , c √ λ − , c √ λ − − γ − √ λγ + √ λ ! − c + √ λ (cid:19) . Note that here ( · ) n is the Pochhammer symbol, that is,( x ) n = n − Y k =0 ( x + k )for a complex number x . The present paper proposed a totally new framework to evaluate theheterogeneous risks in whole-life insurance. We have employed a large-agent limit which is analogous to the thermodynamic one, and boththe life-time and the surrender-time are modelled by the killing time ofthe diffusion process, while the cash-flow is evaluated by the Laplacetransform. The two specific models have shown the potential of ourframework. ven though we have worked only on the determination of thelevel-premium, the proposed framework can be used for more gen-eral cases, including forward-looking models where a mean-field typeapproximation can work. References [1] Adams, C. J., Donnelly, C. and Macdonald, A. S. The impact of knownbreast cancer polygenes on critical illness insurance. Scandinavian Actu-arial Journal , Volume 2015, Issue 2 141–171 , 2015.[2] Akerlof, G. A. The market for “lemons”: quality uncertainty and themarket mechanism. The Quarterly Journal of Economics , pages 488–500,1970[3] Albizzati, M-O. and Geman, H. Interest Rate Risk Management andValuation of the Surrender Option in Life Insurance Policies The Journalof Risk and Insurance Vol. 61, No. 4 (1994), pp. 616-637.[4] Ballotta, L. Eberlein, E. Schmidt, T. and Zeineddine, R. Variable an-nuities in a L´evy-based hybrid model with surrender risk, QuantitativeFinance , 20:5, (2020), 867-886,[5] Bluhm, W. F. Cumulative antiselection theory, Transactions of Societyof actuaries ,34 (1982).[6] Cawley, J. and Philipson, T. An Empirical Examination of Barriers toTrade in Insurance” American Economic Review , 89 (4) (1999), 827-846.[7] Ikeda, N. and Watanabe, S. Stochastic Differential Equations and Diffu-sion Processes , 2nd edtion, Kodansha-North Holland, 1989.[8] Jones, B. L. A Model for Analyzing the Impact of Selective Lapsation onMortality, North American Actuarial Journal , 2:1 (1998), 79–86.[9] Karatzas, I. and Shreve, S. Brownian Motion and Stochastic Calculus ,2nd edition. Springer, 1991.[10] Le Courtois, O. and Nakagawa, H., On surrender and default risks.Math. Finance 23 (1) (2013), 143–168.[11] Loisel, S. and Milhaud, X. From deterministic to stochastic surrenderrisk models: Impact of correlation crises on economic capital, EuropeanJournal of Operational Research , Volume 214, Issue 2, 16 (2011), Pages348-357.[12] de Meza D, and Webb, D.C. Advantageous Selection in Insurance Mar-kets The RAND Journal of Economics Vol. 32, No. 2 (2001), pp. 249–262.[13] Milhaud, X., S. Loisel, and V. Maume-Deschamps, Surrender Triggersin Life Insurance: Classification and Risk Predictions, Bulletin Franc¸caisd’Actuariat, 11(22), (2011):5-48,[14] Gatzert, N. Hoermann, G. and Schmeiser, H. The Impact of the Sec-ondary Market on Life Insurers’ Surrender Profits The Journal of Riskand Insurance Vol. 76, No. 4 (2009), pp. 887-908. 15] Rogers, L.C.G. and Williams, D. Diffusions, Markov Processes and Mar-tingales Volume 1:Foundations Cambridge The Quar-terly Journal of Economics , Vo. 90, No.4 (1976) pages 629–649.[17] Vaupel, J.W., Manton, K.G. and Stallard, E. The impact of hetero-geneity in individual frailty on the dynamics of mortality. Demography Theoretical Population Biology 11 (1977), 37–48.[20] Yashin, A.I., Manton, K.G., and Vaupel, J.W. Mortality and Aging in aHeterogeneous Population: A Stochastic Process Model with Observed andUnobserved Variables, Theoretical Population Biology 27, (1985) 154–75. A Appendix Lemma A.1. We have that E [ e − rζ i { ζ i ≤ t } ] = Z t e − rs E [ V ( X is ) e − R s V ( X iu ) du ] ds. (35) Proof. Let G ( u ) = P ( ζ i > u ) for u ≥ 0. Then, by taking the expecta-tion of both sides of (7), G ( s ) = E [ e − R s V ( X u ) du ] , s ≥ . Since e − R s V ( X u ) du is differentiable in s almost surely and uniformlybounded by 1, we see that G is also differentiable and G ′ ( s ) = − E [ V ( X is ) e − R s V ( X iu ) du ] , s > . Then, since E h e − rζ i { ζ i ≤ t } i = Z t e − rs P ( ζ i ∈ ds )= − Z t e − rs dG s = − Z t e − rs G ′ ( s ) ds, we get (35).we get (35).