On a stochastic nonlocal system with discrete diffusion modeling life tables

On a stochastic nonlocal system with discrete diffusion modeling life tables Tom´as Caraballo1, Francisco Morillas2, Jos´e Valero3 1Universidad de Sevilla, Departamento de Ecuaciones Diferenciales y An´alisis Num´erico Apdo. de Correos 1160, 41080-Sevilla Spain E.mail: [email protected] 2Universitat de Val`encia, Departament d’Economia Aplicada, Facultat d’Economia, Campus dels Tarongers s/n, 46022-Val`encia, Spain. E.mail: F[email protected] 3Universidad Miguel Hern´andez de Elche, Centro de Investigaci´on Operativa Avda. Universidad s/n, Elche (Alicante), 03202, Spain. E.mail: jv[email protected] Dedicated to Prof. Bjorn Schmalfuss on the occasion of his 65-th birthday Abstract In this paper we study an stochastic system of differential equations with nonlocal discrete diffusion. For two type of noises we study the existence of either positive or probability solutions. Also, we analyse the asymptotic behaviour of solutions in the long term, showing that under suitable assumptions they tend to a neighborhood of the unique deterministic fixed point. Finally, we perform numerical simulations and discuss the application of the results to life tables for mortality in Spain. AMS Subject Classification (2020): 34A12,34C60, 34D05, 34D20, 34D23, 60J65, 62P05, 65C99 Key words: stochastic equations, ordinary differential equations, nonlocal problems, discrete diffusion, life tables 1 Introduction Diffusion processes can be sometimes modeled by nonlocal operators, which take into account the influence of the value of the variable in the whole domain (or in part of them) and not only in the current point. For example, we can consider 1 the following nonlocal parabolic problem ut(x, t) = ZRN J(x−y) (u(y, t)−u(x, t)) dy,x∈Ω, t > 0, u(x, t) = g(x, t) , x /∈Ω, t > 0, u(x, 0) = u0(x) , x∈Ω, where Ω ⊆Ris some interval (see [6], [13], [32]). Also, nonlinear nonlocal problems modeling phase transitions have been studied for example in [4], [5]. In [29] we considered the following nonlocal problem with discrete diffusion: d dtui(t) = X r∈D ji−rur(t)−ui(t) + X r∈Z\D ji−rgr(t) , i∈D,t > 0,(1) ui(0) = u0 i,i∈D, where D={m1, . . . , m2},m1< m2,mi∈Z, and j:Z→R+. The term Pr∈Dji−rur(t)−ui(t) can be seen as a space discretization of the integral RRNJ(x−y) (u(y, t)−u(x, t)) dy (see e.g. [32]). This problem was shown to be appropriate to model dynamical life tables by performing numerical simulations with data of mortality from Spanish population. In actuarial or demographic sciences life tables are used to study some biometrics functions, which are related, for example, to the probability of survival or death. They are important in order to calculate the insurance premium or to analyse the sustainability of the social welfare system The predictions for mortality in Spain obtained from model (1) were shown to be adequate for a period of about three years but no longer. This is why in [30] we considered the following model with delay: d dtui(t) = X r∈DZ0 −h ji−rur(t+s)αi(s)dµ(s)−ui(t) (2) +X r∈Z\DZ0 −h ji−rgr(t+s)αi(s)dµ(s), i∈D, t > 0, ui(τ+s)≡φi(s) , i∈D,s∈[−h, 0], where αi: [−s, 0] →R+and dµ(s) = ξ(s)ds being ξ(·) a probability density. In [30] it was shown that with this model the prediction horizon can be extended up to 8 years. In addition, it gives coherent values, in magnitude, when comparing it with other classical techniques such as the Lee-Carter model, up to 18 years. However, in the real world there is always some kind of noise which can be intrinsic to the process under study or that can be related to errors in the observed data. Hence, we will consider now the following stochastic system of 2 differential equations: d dtui(t) = X r∈D ji−rur(t)−ui(t) + X r∈Z\D ji−rgr(t) + bσi(ui(t))dwi dt ,i∈D,t > 0, (3) ui(0) = u0 i,i∈D, where wi(t) are independent Brownian motions and b > 0 is the intensity of the white noise. We will assume the following assumptions on the kernel {ji}and the functions {gi(·)}: (H1) ji≥0 for all i∈Z. (H2) Pi∈Zji= 1. (H3) g∈C([0,+∞), l∞ 2), where l∞ 2={(ui)i∈Z\D: supi∈Z\D|ui|}. We will consider two specific type of noises. Namely: 1) σi(v) = v(linear case); 2) σi(v) = v(1−v). The choice of the noise in the second case is motivated by the fact that we are interested in studying variables like the probability of death, which take values in the interval [0,1]. Noises of this type have been used for example for the logistic equation [22]. Our main aim in this paper is to present some theoretical results concerning the properties of solutions of the stochastic system. First, we prove in the linear case the existence of a unique positive solution whenever the initial datum is positive. Second, we establish in the second case that if the initial condition takes values in the interval (0,1), then the solution remains in this interval for any future time. Third, we analyse the asymptotic behaviour of solutions, showing under certain assumptions that for large times they remain in a neighborhood of the unique fixed point of the deterministic system. Finally, we perform numerical simulations of solutions and analyse the efficiency of model (3) for prediction of mortality using data from Spain. 2 Properties of solutions In this section, under certain conditions on the functions of the problem, we will establish first that for any positive initial condition there exists a unique globally defined solution which is almost sure positive. Then we will prove that if the initial condition lies in the interval (0,1), then the solution remains there for all positive times. This is important for the application to life tables. Throughout this paper we will assume for the kernel {ji}and the functions {gi(·)}the conditions (H1) −(H3) given above. 3 System (3) can be rewritten in the matricial form du dt =Bu (t) + h(t) + bσ(u(t))dw dt ,(4) u(0) = u0∈Rm, where u(t) = (um1(t), ..., um2(t))T, w (t) = (wm1(t), ..., wm2(t))T,h(t) = (hm1(t), ..., hm2(t))T,m=m2−m1+ 1, σ:Rm→Rm×mand (σ(v))ij = 0 if i6=j, (σ(v))ii =σi(vi) for i=m1, ..., m2, hi(t) = X r∈Z\D ji−rgr(t), B=          j0j−1j−2· · · · · · jm1−m2 j1j0j−1j−2· · · jm1−m2+1 j2j1j0j1· · · jm1−m2+2 . . .............. . . jm2−m1−1· · · j2j1j0j1 jm2−m1· · · · · · j2j1j0          . 2.1 A linear noise First, we will consider a standard linear noise, that is, we study the system d dtui(t) = X r∈D ji−rur(t)−ui(t) + X r∈Z\D ji−rgr(t) + bui(t)dwi dt ,i∈D,t > 0, (5) ui(0) = u0 i,i∈D. Denote Rm +={v∈Rm:vj>0 for all j}.We will prove the existence of global positive solutions. For this aim we follow a standard procedure in stochastic differential equations (see e.g. [10], [11], [26]). Lemma 1 Assume that gr(t)≥0for all rand t. Then for any u0∈Rm +there exists a unique globally defined solution u(·)such that u(t)∈Rm +almost sure for t≥0. Proof. The existence of a unique locally defined solution for problem (4) follows from standard results as the functions involved in the systems are locally Lipschitz (see e.g. [15, Chapter 2, p.45]). If the solution is defined in the maximal interval [0, τe), then we need to establish that τe= +∞and that u(t)∈Rm +for all t≥0 a.s. 4 Let k0>0 be large enough so that u0 i∈[1/k0, k0] for all i∈D. As usual, we define the stopping time τk= inf{t∈[0,∞) : ui(t)6∈ (1 k, k) for some i∈D}. This sequence is increasing as k%+∞. If τ∞= limk→+∞τk= +∞a.s., then τe= +∞and u(t)∈Rm +for t≥0 almost sure, proving the assertion. By contradiction, assume the existence of T, ε > 0 such that P(τ∞≤T)> ε. In such a case there exists k1≥k0for which P(τk≤T)≥εfor k≥k1. Further, we consider the C2-function V:Rm +→R1 +given by V(u) = X i∈D (ui−1−log(ui)). For u(t)∈Rm +by Itˆo’s formula we have dV (u(t)) = X i∈D1−1 ui(t) X r∈D ji−rur(t)−ui(t) + X r∈Z\D ji−rgr(t) dt (6) +X i∈D 1 2b2dt +X i∈D1−1 ui(t)bui(t)dwi(t) =I(t)dt +X i∈D Ii(t)dwi(t), where Ii(t) = 1−1 ui(t)bui(t). By using (H1) −(H3) and gr(t)≥0 the first term is estimated as follows for u(t)∈Rm +: I(t)≤X i∈DX r∈D ji−rur(t)−X i∈D ui(t) + m+X i∈DX r∈Z\D ji−rgr(t) + m 2b2 −X i∈D 1 ui(t) X r∈D ji−rur(t) + X r∈Z\D ji−rgr(t)  ≤KT, where we have used that X i∈D 1 ui(t) X r∈D ji−rur(t) + X r∈Z\D ji−rgr(t) ≥0, 5 X i∈DX r∈D ji−rur(t)−X i∈D ui(t) = X r∈DX i∈D ji−rur(t)−X i∈D ui(t) ≤X r∈D ur(t)−X i∈D ui(t) = 0, X i∈DX r∈Z\D ji−rgr(t)≤CTX i∈DX r∈Z\D ji−r≤mCT. Integrating in (6) over (0, τk∧T) and taking expectations we obtain that 0≤EV(u(τk∧T)) ≤V(u0) + EZτk∧T 0 KTdt +EZτk∧T 0X i∈D1−1 ui(t)bui(t)dwi(t) =V(u0) + KTE(τk∧T)≤V(u0) + KTT. Let Ωk={ω:τk≤T}, which satisfies P(Ωk)≥εfor k≥k1. For any ω∈Ωkthere is i∈Dsuch that either ui(τk, ω) = kor ui(τk, ω) = 1/k, which implies that V(u(τk∧T, ω)) ≥(k−1−log(k)) ∧(1 k−1 + log(k)). Hence, V(u0)+KTT≥E(1ΩkV(u(τk∧T))) ≥ε((k−1−log(k))∧(1 k−1+log(k))) = εR(k), where 1Astands for the indicator function of the set A. Passing to the limit as k→+∞we get a contradiction as R(k)→+∞. Corollary 2 Let u0, v0∈Rmbe two initial conditions satisfying u0 i> v0 ifor any i∈D. Also, gu, gv∈C([0,+∞), l∞ 2)are such that gu i(t)≥gv i(t)for all i∈Dand t≥0. Then, ui(t)> vi(t), for all i∈Dand t≥0, where u(·),v(·) are the unique solutions to problem (3) corresponding to {u0, gu}and {v0, gv}, respectively. Proof. It follows immediately by defining the function w(·) = u(·)−v(·) and applying Lemma 1 to it. 2.2 A nonlinear noise Second, we are going to consider the system d dtui(t) = X r∈D ji−rur(t)−ui(t) + X r∈Z\D ji−rgr(t) + bui(t)(1 −ui(t))dwi dt , (7) ui(0) = u0 i, 6 where i∈D,t > 0. We are interesting in proving that if the components of the initial condition lie within the interval (0,1) and the functions gr(t) take values in [0,1], then the components of the solution u(t) also lie in (0,1) for any t≥0 a.s. This is important for applying the results to life tables, as in such situation the variable is a probability. We start with the deterministic case. Lemma 3 Let b= 0. Assume that gr(t)∈[0,1] for all rand t. Then for any u0∈Rmsuch that u0 i∈[0,1], for any i∈D, the unique solution u(·)to (7) satisfies that ui(t)∈[0,1] for all i∈Dand t≥0. Proof. The fact that ui(t)≥0 was proved in [29], so it remains to check that ui(t)≤1. Let z+= max{z, 0}for a real number z. Multiplying the equation by (u−1)+and using the equality [1, Lemma 2.2] dy dt , y+Rm =1 2 d dt  y+  2 Rm we have 1 2 d dt   (u(t)−1)+   2 Rm =X i∈DX r∈D ji−rur(t) (ui(t)−1)+ +X r∈Z\D ji−rgr(t) (ui(t)−1)+−X i∈D ui(t) (ui(t)−1)+ =X i∈DX r∈D ji−r(ur(t)−1)(ui(t)−1)++X i∈DX r∈D ji−r(ui(t)−1)+ +X i∈DX r∈Z\D ji−rgr(t) (ui(t)−1)+−X i∈D (ui(t)−1) (ui(t)−1)+−X i∈D (ui(t)−1)+. Using gr(t)≤1 and (H1) −(H2) we obtain that X i∈DX r∈D ji−r(ui(t)−1)++X i∈DX r∈Z\D ji−rgr(t) (ui(t)−1)+−X i∈D (ui(t)−1)+ ≤X i∈D X r∈D ji−r+X r∈Z\D ji−r (ui(t)−1)+−X i∈D (ui(t)−1)+= 0. 7 Hence, by (H2) we deduce that 1 2 d dt   (u(t)−1)+   2 Rm ≤X i∈DX r∈D ji−r(ur(t)−1)(ui(t)−1)+ ≤X i∈DX r∈D ji−r(ur(t)−1)+(ui(t)−1)+ ≤1 2X i∈DX r∈D ji−r(ur(t)−1)+2+(ui(t)−1)+2 ≤  (u(t)−1)+   2 Rm. Thus,   (u(t)−1)+   2 Rm≤  u0−1+   2 Rme2t, so u(t)≤1 for all t≥0. Further, we study system (7). Lemma 4 Assume that gr(t)∈[0,1] for all rand t. Then for any u0∈Rm such that u0 i∈(0,1), for any i∈D, the unique solution u(·)to (7) satisfies almost sure that ui(t)∈(0,1) for all i∈Dand t≥0. Proof. Let k0>0 be large enough so that u0 i∈[1/k0,1−1/k0] for all i∈D. We define now the stopping time τk= inf{t∈[0,∞) : ui(t)6∈ (1 k,1−1 k) for some i∈D}. This sequence is increasing as k%+∞. If τ∞= limk→+∞τk= +∞a.s., then 0< ui(t)<1 almost sure for t≥0, proving the assertion. By contradiction, assume the existence of T, ε > 0 such that P(τ∞≤T)> ε. In such a case there exists k1≥k0for which P(τk≤T)≥εfor k≥k1. Let K0={u= (um1, ..., um2)∈Rm +: 0 < ui<1}. We define the C2-function V:K0→R1 +given by V(u) = −X i∈D (log(1 −ui) + log(ui)). 8 For u(t)∈K0by Itˆo’s formula we have dV (u(t)) = X i∈D1 1−ui(t)−1 ui X r∈D ji−rur(t)−ui(t) + X r∈Z\D ji−rgr(t) dt +X i∈D 1 2u2 i(t) + (1 −ui(t))2b2dt +X i∈D b(2ui(t)−1)dwi(t) =I(t)dt +X i∈D Ii(t)dwi(t), where Ii(t) = b(2ui(t)−1). We note that X i∈D−1 ui X r∈D ji−rur(t)−ui(t) + X r∈Z\D ji−rgr(t) ≤X i∈D 1. Then the term I(t) is estimated as follows when t≤τk: I(t)≤X i∈D1 1−ui(t) X r∈D ji−rur(t)−ui(t) + X r∈Z\D ji−r +X i∈D1 + b2 ≤X i∈D1 1−ui(t) X r∈D ji−r(ur(t)−1) + X r∈D ji−r+X r∈Z\D ji−r−ui(t) +m1 + b2 ≤X i∈D1 1−ui(t)(1 −ui(t)) + m1 + b2=m2 + b2. Integrating over (0, τk∧T) and taking expectations we infer that 0≤EV(u(τk∧T)) ≤V(u0) + EZτk∧T 0 m2 + b2dt +EZτk∧T 0X i∈D b(2ui(t)−1)dwi(t) =V(u0) + m2 + b2E(τk∧T)≤V(u0) + m2 + b2T. Let Ωk={ω:τk≤T}, which satisfies P(Ωk)≥εfor k≥k1. For any ω∈Ωkthere exists i∈Dsuch that either ui(τk, ω) = 1 kor ui(τk, ω)=1−1 k, so V(u(τk∧T, ω)) ≥log(k)−log(1 −1 k). Thus, V(u0) + m2 + b2T≥E(1ΩkV(u(τk∧T))) ≥εlog(k)−log(1 −1 k). Passing to the limit as k→+∞we arrive at a contradiction. 9 model does not take into account the improvements in mortality, which, as it was already mentioned, is the essential nature of dynamic mortality models. It is easy to obtain an initial correction of the estimates using the average annual improvement rate, IR, which is calculated as the arithmetic mean of the improvement rates by age: IRx=qt x qt+1 x . We have obtained that IR = 0.03.Figure 2 shows the new graphic after this correction. Now, the set of corrected predictions lies between the observed values in 2010 and 2019. Figure 2 shows how a simple correction provides reasonable estimates in magnitude for a short period of time. However, it can be seen that these estimates are not equally adequate for all ages. This is due to the fact that the improvement in mortality does not affect all ages in the same way and, therefore, it is necessary to introduce appropriate correction rates for each age, for example incorporating the history of mortality through delay terms as has benn carried out, in the deterministic system, in the previous work [30]. In this sense, the study of the stochastic model with delay would be an improvement of the current work. This is outside the scope of this study, which is a first approximation to the stochastic model, and will be considered in future works. Figure 2: Ensemble of estimates of mortality Acknowledgments This work has been partially supported by the projects PGC2018-096540- B-I00 (Spanish Ministry of Science, Innovation and Universities), PID2019- 108654GB-I00 (Spanish Ministry of Science and Innovation), P18-FR-4509 (Junta de Andaluc´ıa) and by the Generalitat Valenciana, project PROMETEO/2021/063.. 16 References [1] J.M. Amig´o, A. Gim´enez, F. Morillas and J. Valero, Attractors for a lattice dynamical system generated by non-newtonian fluids modelling suspensions, Internat. J. Bifur. Chaos. 20 (2010) 2681-2700. [2] M. Ayuso, H. Corrales, M. Guill´en, A.M. Perez-Mar´ın and J.L. Rojo, Estad´ıstica Actuarial Vida, (Barcelona, Universitat de Barcelona Edicions, 2007). [3] E. Barbi, F. Lagona, M. Marsili, J.W. Vaupel and W.W. Kenneth, The plateau of human mortality: demography of longevity pioneers, Science 360 (2018) 1459-1461. [4] P. Bates, P. Fife, X. Ren and X. Wang, Travelling waves in a convolution model for phase transitions, Arch. Ratio. Mech. Anal. 138 (1997), 105–136. [5] P. Bates and A. Chmaj, An integrodifferential model for phase transitions: stationary solutions in higher dimensions, J. Stat. Phys. 95 (1999) 1119- 1139. [6] M. Bogoya and C.A. G´omez S., Discrete model of a nonlocal diffusion equation, Rev. Colombiana Mat. 47 (2013) 83-94. [7] A.J.G. Cairns, D. Blake and K. Dowd, Modelling and management of mortality risk: a review. Scandinavian Actuarial Journal 2008 (2008) 79-113. [8] A.J.G. Cairns, D. Blake, K. Dowd, G.D. Coughlan, D. Epstein and M. Khalaf-Allah, Mortality density forecasts: an analysis of six stochastic mortality models, Insurance: Mathematics and Economics 48 (2011) 355-367. [9] A.J.G. Cairns, D. Blake, K. Dowd, G.D. Coughlan, D. Epstein, A. Ong and I.A. Balevich, Quantitative comparison of stochastic mortality models using data from England and Wales and the United States. North American Actuarial Journal 13 (2009) 1-35. [10] T. Caraballo, M. El Fatini, R. Pettersson and R. Taki, A stochastic SIRI epidemic model with relapse and media coverage, Discrete Contin. Dyn. Syst. Ser. B.23 (2018) 3483-3501. [11] T. Caraballo, M. El Fatini, M. El Khalifi, R. Gerlach and R. Pettersson, Analysis of a stochastic distributed delay epidemic model with relapse and Gamma distribution kernel, Chaos Solit. Fractals 133 (2020) 109643. [12] J.B. Copas and S. Haberman Non-parametric graduation using kernel methods, J. Inst. Actuaries 110 (1983) 135-156. [13] C. Cortazar, M. Elgueta and J.D. Rossi, Nonlocal diffusion problems that approximate the heat equations with Dirichlet boundary conditions, Israel J. Math. 170 (2009) 53-60. 17 [14] E. Dolgin, Longevity data hint at no natural limit on lifespan, Nature 559 (2018) 14-15. [15] I.I. Gihman and A.V. Skorohod, Stochastic Differential Equations, (Berlin, Springer-Verlag, 1972). [16] J. Gavin, S. Haberman and R. Verrall, Graduation by Kernel and adaptive kernel methods with a boundary correction, Transactions of the Society of Actuaries, XLVII (1995) 1-37. [17] B. Gompertz, On the nature of the function of the law of human mortality and on a new mode of determining the value of life contingencies, Transactions of The Royal Society 115 (1825) 513-585. [18] P. Hall and T. Wehrly, A geometrical method for removing edge effects from kernel-type nonparametric regression estimators, J. Am. Stat. Assoc. 86 (1991) 665-672. [19] L.M.A. Helligman and J.H. Pollard, The age pattern of mortality, J. Inst. Actuaries,107 (1980) 49-82. [20] P.E. Kloeden, E. Platen and H. Schurz, Numerical Solution of SDE through Computer Experiments, (Berlin, Springer, 1994). [21] Life tables: Methodolgy. Instituto Nacional de Estad´ıstica (Statistics National Institute of Spain). [22] M. Liu and K. Wang, A note on stability of stochastic logistic equation, Appl. Math. Letters. 26 (2013) 601-606. [23] R. Lee and L. Carter, Modelling and forecasting US mortality, J. Am. Stat. Assoc. 87 (1992) 659-671. [24] W. Makeham, On the law of mortality and the construction of annuity tables, Journal of the Institute of Actuaries 8(1869) 301-310. [25] X. Mao, Stochastic Differential Equations and Applications, (Oxford, Woodhead Publishing, 2007). [26] X. Mao, G. Marion and E. Renshaw, Environmental Brownian noise suppresses explosions in population dynamics, Stoch. Process Their Appl. 97 (2002) 95-110. [27] T. Missov, Gamma-Gompertz life expectancy at birth, Demographic Research 28 (2013) 259-270. [28] F.G. Morillas and I.B. Sampere, Using wavelet techniques to approximate the subjacent risk of death, In Modern Mathematics and Mechanics, V.A. Sadovnichiy and M.Z. Zgurovsky eds. (Cham, Switzerland, Springer International Publishing, 2019), Chapter 28, pp. 541-557. 18 [29] F. Morillas and J. Valero, On a nonlocal discrete diffusion system modeling life tables, RACSAM 108 (2014) 935-955. [30] F. Morillas and J. Valero, On a retarded nonlocal ordinary differential system with discrete diffusion modeling life tables, Mathematics 9(2021) 220. [31] J.M. Pav´ıa, F.G. Morillas and J.C. Bosch-Rodr´ıguez, Using parametric bootstrap to introduce and manage uncertainty: replicated loaded insurance life tables. North American Actuarial Journal 23 (2019) 434-446. [32] M. P´erez-Llanos and J.D. Rossi, Numerical approximations for a nonlocal evolution equations, SIAM J. Numer. Anal. 49 (2012) 2103-2123. 19