Aggregate Markov models in life insurance: properties and valuation

In multi-state life insurance, an adequate balance between analytic tractability, computational efficiency, and statistical flexibility is of great importance. This might explain the popularity of Markov chain modelling, where matrix analytic methods allow for a comprehensive treatment. Unfortunately, Markov chain modelling is unable to capture duration effects, so this paper presents aggregate Markov models as an alternative. Aggregate Markov models retain most of the analytical tractability of Markov chains, yet are non-Markovian and thus more flexible. Based on an explicit characterization of the fundamental martingales, matrix representations of the expected accumulated cash flows and corresponding prospective reserves are derived for duration-dependent payments with and without incidental policyholder behaviour. Throughout, special attention is given to a semi-Markovian case. Finally, the methods and results are illustrated in a numerical example.


Introduction
In this paper, we propose a new class of multi-state models, the so-called aggregate Markov models, and study the valuation of life insurance liabilities for this class of models.The results established range from a characterization of the fundamental martingales to genuine computational schemes for the quantities of interest, including prospective reserves.In contrast, the companion paper [2] deals with the statistical aspects.

INTRODUCTION
The classic approach to multi-state modelling consists of Markov chain modelling, where the process governing the state of the insured Z is taken to be a time-inhomogeneous Markov chain and the payments between the insured and the insurer are required to consist of deterministic sojourn and transition payments.This approach dates back to at least [17], but has recently been given new life via matrix analytic methods, see [6,1].Although Markov chain modelling is attractive due to its inherent simplicity, it suffers from a number of defects, not at least its inability to properly capture duration effects.In recent years, semi-Markov modelling has therefore gained considerable attention.In semi-Markov modelling, see [18,15,11,10], only the joint process (Z, U ), where U denotes the time spent in the current state, is assumed to be Markovian, and the sojourn and transition payments are allowed to depend on U .Unfortunately, semi-Markov modelling entails less analytical tractability and increased computational load.
In an aggregate Markov model, each observable state is assumed to consist of multiple unobservable sub-states, and only the full model consisting of all sub-states is assumed to be Markovian.Different to the classic Markov chains, which require the Markov property already for the observable states, aggregate Markov models are non-Markovian and thus flexible, yet they retain most of the analytical tractability of Markov chains.In particular, matrix analytic methods related to inhomogeneous phase-type distributions, confer with [3], are applicable and lead to a unifying and transparent treatment.
Phase-type distributions have a long history of extensive use in applied probability.They have been employed in areas such as queueing theory [4,22,23,20], actuarial science [5,7], and telecommunications [4,20], where the phase-type assumption leads to exact and in many cases explicit formulas for properties such as waiting time distributions, queue length, ruin probabilities, and buffer overflows.
Both homogeneous as well as inhomogeneous phase-type distributions are dense in the class of distributions on the positive reals, confer with [7], and therefore able to approximate any non-negative distribution arbitrarily close -in the sense of weak convergence as the number of phases increases to infinity.Hence the class of phase-type distributions has been considered as striking a balance between tractability and generality.Inhomogeneous phase-type distributions may be used instead of homogeneous ones, and this might be particularly relevant if the tail behaviour is known to be different from exponential, see [3].
In the area of queuing theory, so-called quasi-birth-and-death (QBD) processes have been extensively studied, confer with [20] and references therein, as a model for the number of customers in a queue.They constitute the time-homogeneous analogue to the aggregate Markov models considered here.
The first main contribution of the paper is an explicit characterization of the martingales for the associated counting processes N , which reveals that aggregate Markov models may be highly non-Markovian.For many practical purposes, less might suffice.To this end, we provide a sort of reset property under which the aggregate Markov model is actually semi-Markovian.The second main contribution of the paper are matrix representations for the expected accumulated cash flows, and hereby the prospective reserves, for duration-dependent payments with and without incidental policyholder behaviour.Special attention is given to the case where the payments are duration independent; here our results indicate that aggregate Markov modelling may hold a competitive advantage over semi-Markov modelling.
The remainder of the paper is structured as follows.Section 2 provides some background, with Subsection 2.1 devoted to the basics of inhomogeneous phase-type distributions and Subsection 2.2 to the basics of multi-state modelling in life insurance.These subsections might be passed over by readers who are familiar with the subject matter.In Section 3, we introduce the aggregate Markov models and state the aforementioned reset property.The main contributions take place in Section 4 and Section 5.The former is devoted to the distributional properties of Z, including the characterization of the fundamental martingales, while the latter deals with the valuation of life insurance liabilities and contains, in particular, matrix representations of the expected accumulated cash flows.To showcase the practical potential of aggregate Markov modelling, Section 6 contains a numerical example.Finally, Section 7 concludes.Proofs may be found in Appendix A.

Preliminaries
Before introducing the setting of the paper, we provide some background.Subsection 2.1 contains a short review on imhomogeneous phase-type distributions, which play a critical role later.This review is followed by Subsection 2.2, which collects some insights on multi-state modelling in life insurance, in particular in relation to Markov chain and semi-Markov modelling, hereby motivating our aggregate setup.The actual presentation of our setup is postponed to Section 3.
In what follows and throughout the paper, we denote the product integral of a square matrix function A(x) as where I is the identity matrix.Under suitable regularity conditions, it may equivalently be cast as the solution to Kolmogorov's forward and backward differential equations: In other words, the computation of a product integral is equivalent to solving a system of first order, in general inhomogeneous, ordinary differential equations.In particular, a standard Runge-Kutta scheme may be employed.Alternatively, the solution may be approximated by a product of matrix exponentials, using a piece-wise constant square matrix function.Further details may be found in [6, Section 2].For a survey on the probabilistic and statistical aspects of product integration, we refer to [14] and, concerning applications to life insurance, also [21,6].

Inhomogeneous phase-type distributions
In this subsection, we review the notion of inhomogeneous phase-type (IPH) distributions introduced in [3].Consider a smooth and suitably regular time-inhomogeneous Markov jump process X = {X(t)} t≥0 on the finite state space J = {1, . . ., J − 1, J}, where the states {1, . . ., J − 1} are transient while J is absorbing.The transition intensity matrix function where T (t) is a sub-intensity matrix function consisting of transition rates between the transient states and t(t) = −T (t)1 1 1 J is a column vector of transition rates to the absorbing state, the so-called exit rate vector function.Further, assume that P(X(0) = J) = 0 and denote by π = (π 1 , . . ., π J−1 ) the remaining vector of initial probabilities π j = P(X(0) = j).The time until absorption, given by is then said to be an inhomogeneous phase-type distribution with representation (π π π, T T T ), and we write τ ∼ IPH(π π π, T T T ).
The transition probability matrix function P (t, s) = {p ij (t, s)} i,j∈J with elements is given as the product integral of the transition intensity matrix function: The probability density function f (t) and distribution function F (t) of τ may then be obtained through product integrals of the sub-intensity matrix function T (t): From these, one finds the following conditional distribution: which entails that where α α α(s) is given by In other words, for IPH distributions, the overshoot is again IPH-distributed.
Example 2.1.In the case of a single phase, that is J = 2, we have giving the density and distribution functions while the conditional distribution (2.3) takes the form

Multi-state modelling
Insurance contracts may be modelled as a stream of payments B = {B(t)} t≥0 , benefits less premiums, between the insured and the insurer.In life insurance, including health and disability insurance and pensions, the payments depend on the state of the insured, leading to so-called multi-state modelling.In general, the state of the insured Z = {Z(t)} t≥0 is a non-explosive jump process on a typically finite state space J = {1, 2, . . ., J}, J ∈ N, while the payments are typically finite variation processes adapted to the information generated by Z.

Markov chain models
The most classic approach to multi-state modelling is (smooth) Markov chain models, where Z is taken to be a time-inhomogeneous Markov jump process (Markov chain) with suitably regular transition rates ν jk (t), so that Using the standard convention ν jj (t) = − k∈J k̸ =j ν jk (t), the square matrix function with indices ν jk (t) is then the transition intensity matrix function of Z.In addition to the Markov assumption, the payments are assumed to take the form for suitably regular deterministic sojourn payment rates b j (t) and transition payments b jk (t) depending only on time.Here N is the multivariate counting process associated to Z with components N jk = {N jk (t)} t≥0 given by Markov chain modelling dates back to at least [17] and was popularized in [24].
Regarding the valuation of life insurance liabilities, calculating the so-called expected accumulated cash flow A(t, s) is key.For Markov chain models, the expected accumulated cash flow A(t, s) is given by A(t, s) = i∈J 1 (Z(t)=i) A i (t, s), where The transition probabilities, considered as a square matrix function, are given as the product integral of the transition rates, also considered as a square matrix function.In other words, the transition probabilities may be found simply by solving Kolmogorov's forward differential equations.
It is, of course, possible, and of interest, to relax both the Markov assumption as well as the structure of the payments.In doing so, it is critical to strike an adequate balance between analytic tractability, computational efficiency, and statistical flexibility.

Semi-Markov models
A more modern approach is semi-Markov modelling, see for instance [15,11,10].There are two major differences between (smooth) semi-Markov modelling and (smooth) Markov chain modelling.First, the jump process Z describing the state of the insured is no longer required to be Markovian; rather, (Z, U ) is assumed Markovian, where U = {U (t)} t≥0 is the duration since the last transition given by Therefore, the model can no longer be described by transition rates that solely depend on time.Instead, the transition rates are now functions of both time and duration, written ν jk (t, u).
Second, the payments take the more general form for suitably regular deterministic sojourn payment rates b j (t, u) and transition payments b jk (t, u) depending on time and duration.
For semi-Markov models, the expected accumulated cash flow A(t, s) depends on both the current state and current duration.To clarify, it may actually be decomposed according to A(t, s) = i∈J 1 (Z(t)=i) A i,U (t) (t, s), where The transition probabilities may be calculated by solving a system of integro-differential equations, confer with [10, Section 3].Numerical methods for integro-differential equations can, generally speaking, be rather intricate.The implementation of semi-Markov models is, therefore, non-trivial and may carry some operational risk.

Aggregate Markov models
In this paper, we introduce a class of aggregate models that, similar to semi-Markov models, allow for added flexibility, such as duration dependence, but avoid some of the aforementioned numerical challenges posed by semi-Markov modelling.
Denote by (T n ) n∈N 0 the jump times of Z, where we employ the convention T 0 = 0. Returning to the case where Z is Markovian with transition rates ν jk (t), recall that We conclude that In other words, the sojourn times follow one-dimensional IPH distributions that are mostly independent of the past history of the jump process.This paper considers instead jump processes with sojourn times admitting conditional IPH distributions of general dimension.Hereby we shall be able to capture, for instance, duration dependence while avoiding the need for intricate numerical methods.

Setup
In this section, we present the general setup of the paper.Subsection 3.1 introduces the probabilistic model for the state of the insured, while Subsection 3.2 introduces the payments between the insured and the insurer.

Probabilistic model
Similar to Subsection 2.2, let Z be a jump process governing the state of the insured, thus taking values in the finite set of (macrostates) J = {1, 2, ..., J}, J ∈ N.This set consists of biometric or behavioural states that are actually observed, for example active, disabled, free-policy, and dead.To allow for added flexibility, to each macrostate we may introduce additional sub-states (microstates) that are not observable.
To be specific, to each macrostate j, a number d j ≥ 1 of microstates are assigned.The resulting state space is therefore and the total number of microstates is d = j∈J d j .Elements of E are generally denoted by bold letters such as j ∈ E. Now introduce a time-inhomogeneous Markov jump process X X X = {X X X(t)} t≥0 = {(X 1 (t), X 2 (t))} t≥0 on the state space E with transition intensity matrix function M (t).Then X 1 (t) keeps track of the macrostate, that is Z(t) = X 1 (t), while X 2 (t) identifies the current microstate contingent on the state of X 1 (t).
The transition intensity matrix function M (t) can be written on the following block form: where M jj (t) are sub-intensity matrix functions of dimension d j ×d j providing transition rates between the microstates of macrostate j, and M jk (t) are non-negative matrix functions of dimension d j × d k providing transition rates from microstates within macrostate j to microstates within macrostate k.
We denote an element of M (t) by µ jk (t), j, k ∈ E. The off-diagonal elements are nonnegative, providing the jump rates between different states, while the diagonal equals the negative of the row sums of the off-diagonal elements.Consequently, rows all sum to zero, so M (t) is a proper transition intensity matrix function.

The column vector function
contains the exit rate function out of macrostate j.This function is non-negative due to M jj (t) being a sub-intensity matrix function.The last equality follows from the row sums of M (t) being zero.
In this paper, we give special attention to that case where M jk , j, k ∈ J , j ̸ = k, is a matrix of rank one on the form where β β β jk (t) is a d j -dimensional non-negative column vector function and Here β β β jk (t) provides the vector of jump rates from the microstates of macrostate j to the macrostate k, and 3) holds, we say that k has the reset property from state j.If all states k have the reset property from all states j ̸ = k, we simply say the reset property is satisfied.
Another way of writing (3.3) is Remark 3.2.Focusing only on the transition from macrostate j to macrostate k, we look at the following elements of M (t): Here β β β jk (t) is a column vector function of exit rates from states (j, j) in {(j, 1), ..., (j, d j )}.Therefore it seems natural to pair j and j, which explains the seemingly awkward indexation of elements of β β β jk (t).Since a new microstate is picked independently of (j, j) from {(k, 1), ..., (k, d k )}, records of where the process transitioned from are lost upon transition, which explains the term reset property.△ Let F Z = {F Z (t)} t≥0 denote the natural filtration generated by the macrostate process Z. Since only Z is observed, the filtration F Z represents the available information.We may, as previously, associate to Z a multivariate counting process N with components N jk = {N jk (t)} t≥0 given by as well as a marked point process (T n , Y n ) ∞ n=0 with T n the n'th jump time of Z and Y n = Z(T n ); we use the convention T 0 = 0. Disregarding null-sets, the jump process Z, the multivariate counting process N , and the marked point process (T n , Y n ) ∞ n=0 generate the same information.
Although the microstate process X X X is Markovian, this is generally not the case for the macrostate process Z.In this paper, we derive distributional properties of Z by deriving distributional properties of the multivariate counting process N and the marked point process (T n , Y n ) ∞ n=0 .We are especially interested in the special case where the reset property holds.Here it turns out that (Z, U ) becomes Markovian, where U is the duration process defined in (2.4).

Payments
Having specified the probabilistic model, we now turn our attention to the insurance contract itself.Again, we denote by B = {B(t)} t≥0 the payments, benefits less premiums, between the insured and the insurer.We suppose that B takes the form prescribed in (2.5).Furthermore, throughout the paper, we assume a maximal contract time η > 0 such that all sojourn payment rates and transition payments are zero after time η.
In the later stages of the paper, we add another layer of complications by turning to so-called scaled payments that appear in connection with policyholder behaviour such as free-policy conversion and stochastic retirement.To be precise, here we furthermore consider payments where τ is the exercise time of some policyholder option (modelled incidentally) and 0 < ρ(t, j, k) ≤ 1 is a suitable regular deterministic scaling factor.
The remainder of the paper now focuses on deriving distributional properties of the macrostate process Z, establishing computational schemes for relevant expected accumulated cash flows and prospective reserves, and finally relating these findings to existing models and methods in the life insurance literature.

Properties of Z
In this section, we derive some distributional properties of Z.In Subsection 4.1, we consider the general setup and derive the conditional finite-dimensional distributions of the marked point process (T n , Y n ) n∈N 0 associated to Z as well as the predictable compensators of the multivariate counting processes N associated to Z.In Subsection 4.2, we impose the reset property, which we show, by applying the results of Subsection 4.1, leads to (Z, U ) being Markovian.

General results
Since Z is generally not Markovian, we introduce to keep track of the history of Z. Write for a generic realization of S n whenever T n < ∞.Let denote the conditional survival function of T n+1 given S n , and let denote the conditional probability mass function of Y n+1 given (S n , T n+1 ).These quantities determine the distribution of Z.The following result provides a characterization of them within our setup.
Proposition 4.1.The conditional finite-dimensional distributions of the marked point process (T n , Y n ) ∞ n=1 are given by where the d yn -dimensional row vector α α α(s n ) is given by Proof.Please refer to Appendix A.

2). △
The compensators of the multivariate counting process associated to Z, which also determine the distribution of Z, are key quantities in the context of estimation and valuation.
In our setup, they take the following form.
Theorem 4.3.The counting process N jk has (F Z , P)-compensator given by dΛ jk (t) = λ jk (t) dt with Proof.Please refer to Appendix A.

Reset property and semi-Markovianity
Suppose now that the reset property holds, that is M jk (t) satisfies (3.3) for all j ̸ = k.We now make the following observations.From (3.3), we see that for k ̸ = j, is a 1 × 1-dimensional matrix, implying it cancels if appearing in both the numerator and denominator of a fraction.In particular, Combined with the results of Subsection 4.1, this yields the following corollaries.
Corollary 4.4.Assume that the reset property holds.Then the conditional finitedimensional distributions of the marked point process (T n , Y n ) ∞ n=1 are given by where m m m yn (t n+1 ) is obtained from (3.4).
Remark 4.5.The first statement of Corollary 4.4 corresponds to Assume that the reset property holds.Then the counting process N jk has (F Z , P)-compensator given by dΛ jk (t) = λ jk (t) dt with The results show that the general path dependence of Z through α α α(S n ) is significantly reduced whenever the reset property is imposed.We may actually write which shows that the macrostate process Z is a time-inhomogeneous semi-Markov process with transition rates ν jk (t, u), j ̸ = k, which are functions of both time and duration, given by Conversely, the class of aggregate Markov models is quite flexible.In light also of the many prized denseness results for phase type distributions, we therefore suggest the following conjecture: Conjecture 4.7.The class of aggregate Markov models with the reset property is dense in the class of time-inhomogeneous semi-Markovian models with absolutely continuous sojourn time distributions.
In other words, we conjecture that any (smooth) semi-Markovian model can be approximated arbitrarily well by an aggregate Markov model with the reset property simply by letting the number of microstates increase to infinity.The clarification of this conjecture is outside the scope of this paper, but it nevertheless constitutes an interesting research direction.

Valuation
In this section, we consider the valuation of the life insurance liabilities corresponding to the payment process B. In Subsection 5.1, we provide expressions for the expected accumulated cash flows and, hereby, the prospective reserves.In particular, we provide matrix representations that are useful in implementing the models.The expected accumulated cash flows are composed of conditional occupation probabilities, for which we derive formulas in Subsection 5.2.Special emphasize is given to the semi-Markovian case of Subsection 4.2.Finally, in Subsection 5.3 and Subsection 5.4, we investigate the impact of duration-independent payments and the inclusion of policyholder options, respectively.
Throughout the section, the time value of money is described by a deterministic and suitably regular interest rate r(t).As long as financial and insurance risks are assumed independent, the extension to stochastic interest rates is straightforward.

General results
The expected accumulated cash flow A(t, s) valued at time t is given by, confer with Definition 2.2 in [8], so that the prospective reserve reads e − s t r(v) dv A(t, ds). (5.1) If the payments are on the form (2.5), then where the conditional occupation probabilities are given by If the reset property is also satisfied, then Z is a time-inhomogeneous semi-Markovian process and thus A(t, s) = i∈J 1 (Z(t)=i) A i,U (t) (t, s), where Furthermore, it holds that where the transition probabilities are given by For implementation purposes, it may be beneficial to use matrix representations of the expected accumulated cash flow A(t, s) following along the lines of [6], since it allows for more compact and direct computations.
In the general case where the aforementioned reset property is not satisfied, the process Z is non-Markovian, so it is not sensible to form a transition probability matrix function in the usual way.Instead, we form a d-dimensional vector function according to p(t, s, dz) = p j (t, s, dz) j∈E . (5.4) In regards to payments and transition rates, however, the fact that X X X is assumed to be Markovian allows us to follow more closely the approach of [6].In the present setup, we have a set of sojourn payment rates and transition payments that are all identical across microstates (of the same macrostate).Hence, the d-dimensional vector of sojourn payment rates on the micro level is given by where b b b j (t, u) = b j (t, u)1 1 1 d j .The matrices of transition payments must, in a similar fashion, be identical across microstates (of the same macrostate), so that the transition payment matrix function on a micro level is given by where B ij (t, u), i, j ∈ J , j ̸ = i, is a d i ×d j -dimensional matrix with b ij (t, u) in all entries, and B ii (t, u) = 0 is a d i × d i -dimensional matrix of zeroes.Based hereon, we define the reward matrix function as where • denotes the Schur product, that is The expected accumulated cash flow A(t, s) may then be seen to have the following matrix representation: where the original sums over the state space E are reduced to matrix multiplications.
In the case of the reset property, the semi-Markovianity of Z implies that it suffices to consider the transition probabilities of (5.3).Thus, it is sensible to form a J × ddimensional matrix function according to p(t, u, s, dz) = {p ij (t, u, s, dz)} i∈J ,j∈E .
Similarly, we may form the J-dimensional vector of state-wise expected accumulated cash flows, which then can be calculated as follows: where the reward matrices R(t, u) of (5.7) are modified according to (3.3)-(3.4).We can also cast (5.9) as A u (t, ds) = a a a u (t, s) ds, where then is a vector of state-wise expected cash flows.

Conditional occupation and transition probabilities
We now provide calculation schemes for the conditional occupation and transition probabilities.Rather than working directly with these quantities, it turns out to be fruitful to focus instead on In the following, we require the d × d j -dimensional matrices e e e j e e e ′ j , where e e e j is a d j -dimensional column with a one in entry j and otherwise zeroes, and e e e j is a d-dimensional column vector with a one in entry d 1 + • • • + d j−1 + j and otherwise zeroes.Here and in the following, primes denotes matrix transposition.The entries of E j are zero, except in the j'th block row, where they consist of the d j -dimensional identity matrix.Roughly speaking, they allow us to extend a distribution on microstates in a single macrostate to the whole state space E (and vice versa).
Theorem 5.1.It holds that where the auxiliary quantities p j (t, s, z; s n ) are zero for t n ≥ s − z and (5.10) Proof.Please refer to Appendix A.
If the reset property is satisfied, in which case Z is a time-inhomogeneous semi-Markovian process, we can use (4.1) to immediately obtain the following corollary.
Corollary 5.2.Assume that the reset property holds.Then pij (t, u, s, z) is zero for t − u ≥ s − z and for t − u < s − z.
We note that z → pij (t, u, s, z) is continuous on [0, u + s − t) and actually constant on [s − t, u + s − t).If i ̸ = j, then the continuity extends to [0, u + s − t], while The fact that z → pij (t, u, s, z) is constant on [s − t, u + s − t) may be utilized to reduce the computational load when calculating the expected accumulated cash flows.
We conclude this subsection by presenting an algorithm for the computation of expected cash flows in models with the reset property, see the following page.The computational scheme is similar to the algorithm for general semi-Markov models based on Kolmogorov's forward integro-differential equation proposed in [10, Section 3].Both algorithms employ a two-dimensional time and duration grid, and one would therefore expect the computational loads to be comparable.

Duration-independent payments
We now consider the simplifications arising from duration-independent payments, that is, when where According to Theorem 5.1, )e e e j . (5.12) Algorithm 1 Computation of expected cash flows in an aggregate Markov model with the reset property.
Input: Current time t ∈ [0, η), current duration u ∈ [0, t], and a grid T : 1) Calculate initial conditional distributions at time t: 2) Compute transition probabilities for the Markovian state process X X X, by solving Kolmogorov's forward differential equation on T .
ii) Calculate the vector of state-wise expected cash flows for time t ℓ : using numerical integration methods on the grid T for the integral, and where Output: For each ℓ ∈ {1, . . ., n}, a vector of state-wise expected cash flows a a a u (t, t ℓ ).
If the reset property is satisfied, then we are rather interested in A A A u (t, s), which subject to (5.11) reads where This should lead to a significant reduction in computational load since the above simplification allows one to adapt Algorithm 1 to employ only a one-dimensional time grid.For general semi-Markov models, where the computation of transition probabilities relies on Kolmogorov's forward integro-differential equation, such a simplification is not possible.It should be noted, however, that the computation of the term might still be rather burdensome if d is large.To conclude, if the number of microstates per macrostates is not too large, aggregate Markov models might hold a competitive advantage over general semi-Markov models if the payments of interest are durationindependent.
The above discussion relates to duration-independent payments.However, it is also applicable to certain crude duration-dependent payments.This is partly illustrated by the numerical example in Section 6, where we consider a contract stipulating a waiting period.

Policyholder behaviour
We now extend the results of Subsections 5.1-5.2 to include incidental policyholder behaviour such as free-policy conversion and expedited or postponed retirement.The inclusion of policyholder options is quite popular in the life insurance literature, confer with [16,9,10,13], especially for Markov chains.General insights based on change of measure techniques were recently provided in [12].In the following, we adapt the general methods and results of [12] to our setting.
Suppose that the macrostates J can be decomposed as with 1 ∈ J 0 and where the transition intensity matrix function M (t) of the microstate process X X X is composed of block matrix functions satisfying In that case, the macrostate process Z almost surely never hits ∇, and, upon entering the states J 1 , the process never returns to J 0 .Recall that Z(0) ≡ 1, so the process starts in J 0 .We may thus interpret J 0 as the states of the insured prior to exercising their policyholder option and J 1 as the states of the insured after exercising the option.
(The role of the 'dummy' state ∇ will be clear later.)The first hitting time of J 1 , given by then constitutes the exercise time of the option.Since Z almost surely never hits ∇, we may as well take b ∇ (s, u) = 0, b j∇ (s, u) = 0, and b ∇k (s, u) = 0.At exercise, the original contractual payments are scaled with the factor ρ(τ, Z(τ −), Z(τ )), where 0 < ρ(t, j, k) ≤ 1 is some suitably regular deterministic function.The payment process of interest B ρ = {B ρ (t)} t≥0 thus takes the form The scaling factor is typically selected as to maintain actuarial equivalence with respect to a safe-side valuation basis, the so-called technical basis; we just consider it given.The corresponding expected accumulated cash flow A ρ (t, s) valued at time t is The following result is a consequence of [12, Theorem 3.6 and Proposition 3.10].
Proposition 5.3.It holds that where E denotes expectation with respect to another probability measure P. Furthermore, the (F Z , P)-compensators of the counting processes are given by Recall that the compensators determine the distribution of Z. Thus from the expressions for the compensators obtained in Theorem 4.3 and the above proposition, we find that Z under P follows an aggregate Markov model with transition intensity matrix function M (t) composed of block matrix functions (1 − ρ(t, j, k))M jk (t)e e e, j ∈ J 0 , Furthermore, if the reset property is satisfied under P, this is also the case under P.All in all, according to Proposition 5.3 the expected accumulated cash flow and thus also the expected accumulated cash flow A ρ (t, s), can be calculated using the formulas of Subsections 5.1-5.2, but with M (t) replaced by M (t).

Numerical example
We conclude the paper by presenting a numerical example that serves to illustrate the methods presented in Section 5.The probabilistic models, described by transition rates on the micro level, are taken from the numerical example in [2], where aggregate Markov models corresponding to Figure 1 with the reset property are fitted to simulated data on a macro level for different numbers of disability microstates, d 2 .The simulations are based on a (smooth) semi-Markovian disability model employed by a large Danish life insurance company that has been reported to and published by the Danish Financial Supervisory Authority.The only duration effects present in this model concern the rates from the disability state, which also explains why we do not add extra microstates to the active macrostate.The rates from the disability state are, at least after some months, decreasing as functions of duration.
In the following, the analysis of [2] is extended with calculations and comparisons of state-wise expected cash flows and prospective reserves.We focus on coverage which admits duration-dependent payments, namely disability coverage with a waiting period.To be specific, we consider a male of age t = 40 years with a disability annuity of rate 1 per year, starting 3 months after the onset of disability, but only until retirement at age 65.The only non-zero payments function is thus b 2 (s, z), which reads b 2 (s, z) = 1 (s<65) 1 (z>1/4) , and we may therefore set η = 65.The aggregate Markov models fitted in [2] have the reset property and, consequently, the associated macrostate processes are (timeinhomogeneous) semi-Markov and with transition rates according to (4.2).In Figure 2, we provide for illustrative purposes a selection of the estimated transition rates that are particularly relevant for the coverage of interest, but only for initial duration u = 0, 1 years.We refer to the numerical example in [2] for further details.
We emphasize that the particular simple type of duration dependence of the payments allows for simplifications in the computation schemes similar to those from the durationindependent case of Subsection 5.3.Indeed, the vector of state-wise expected cash flows now reads a a a u (40, s) = 1 (u+s−40>1/4) p(40, u, s, where the elements of p(40, s, u, 1/4) can be calculated using Corollary 5.2, confer also with Algorithm 1.Since we only need the transition probabilities at a single (and rather small) end duration z = 1/4, the computational complexity is comparable with that of the duration-independent case, where only z = 0 is needed.
The corresponding vector of state-wise prospective reserves is obtained by discounting and accumulating the vector of state-wise expected cash flows: For the interest rate, we use the forward rate curve published on the 3rd of November, 2022, by the Danish Financial Supervisory Authority.We are implicitly assuming that the financial market is independent of the state process of the insured.We calculate the vectors of state-wise expected cash flows and corresponding prospective reserves across initial states and durations and for the various fits.We also calculate these quantities for the underlying true semi-Markov model, where we need to use Kolmogorov's forward integro-differential equation of [10,Section 3], since this model is not an aggregate Markov model.Figure 3 shows the resulting expected cash flows in the disability state for initial duration u = 0, 1, while Figure 4 shows the prospective disability reserves as functions of the duration since the onset of disability.Since it is unable to capture the duration effects that are present, the Markov chain corresponding to d 2 = 1 performs, as anticipated, very badly.Maybe more surprisingly, the addition of just a single additional disability microstate corresponding to d 2 = 2 leads to significant improvements; this model may already be competitive, depending also on the trade-off between accuracy and computational load.Furthermore, and consistent with our expectations, the accuracy appears to further improve as the number of disability microstates, d 2 , increases.The latter point might be less clear from a purely visual standpoint, especially for d 2 = 7, but it is confirmed by the final log-likelihood values from the statistical analyses in [2], confer also with Table 1.Another interesting observation is the differences in accuracy across regions of initial duration, u.For instance, while the model corresponding to d 2 = 3 outperforms the model corresponding to d 2 = 5 for u ∈ [1, 2], it is far from competitive for smaller u.As discussed in [2], this is mainly due to the irregular (non-smooth) nature of the true transition rates in that region, see also Figure 2.

Concluding remarks
In this paper, we introduce a new class of models for multi-state life insurance, known as aggregate Markov models, with a specific focus on the valuation of liabilities.The use of matrix analytic methods lead to a unifying and transparent treatment of both the probabilistic as well as the computational facets of the models.In contrast, the companion paper [2] is devoted to the statistical aspects.
In summary, aggregate Markov modelling offers a versatile approach for modelling multistate life insurance contracts, retaining the advantageous analytical features of Markov chain modelling while providing added flexibility, including the incorporation of duration effects.Our theoretical results suggest that aggregate Markov modelling may have a competitive advantage over semi-Markov modelling when the contractual payments exhibit a sufficient level of simplicity in their duration dependence.This assertion is further supported by a numerical example.
Future research endeavors could address several open questions.For instance, the resolution of Conjecture 4.7, which proposes that any (smooth) semi-Markovian model can be approximated arbitrarily well by an aggregate Markov model with the reset property, remains unresolved.An affirmative solution, achieved through an efficient construction, would perhaps encourage the wider adoption of aggregate Markov models.
Throughout this study, our focus has been solely on smooth models, meaning models with absolutely continuous sojourn time distributions.However, to unify modelling in Collecting results confirms the identity for n = 0. Now suppose the identity holds for n ∈ N 0 .We want to establish the identity also for n + 1.By assumption, P t n+1 < T n+1 ≤ t n+1 + h, X X X(T n+1 ) = (y n+1 , yn+1 ) S n = s n = α α α(s n , t n+1 , y n+1 )e e e yn+1 α α α(s n )1 1 1 dy n h + o( h).
In particular,

Remark 4 . 2 .
The first statement of Proposition 4.1 corresponds to and ∆(b b b) is a diagonal matrix with the vector b b b as diagonal.This is similar to equations (3.8)-(3.11) in [6].
40 r(v) dv a a a u (40, s) ds.

Figure 3 :
Figure 3: Expected cash flow s → a 2,u (40, s) in the disability state with initial durations u = 0 (left) and u = 1 (right) for different numbers of disability microstates, d 2 , along with the true expected accumulated cash flows.The case d 2 = 1 corresponds to a Markov chain.

Figure 4 :
Figure 4: Prospective disability reserve as a function of duration since onset of disability, u → V 2,u (40), for different numbers of disability microstates, d 2 , along with true prospective reserve.The case d 2 = 1 corresponds to a Markov chain.