Calibration of Transition Intensities for a Multistate Model: Application to Long-Term Care

: We consider a non-homogeneous continuous time Markov chain model for Long-Term Care with ﬁve states: the autonomous state, three dependent states of light , moderate and severe dependence levels and the death state. For a general approach, we allow for non null intensities for all the returns from higher dependence levels to all lesser dependencies in the multi-state model. Using data from the 2015 Portuguese National Network of Continuous Care database, as the main research contribution of this paper, we propose a method to calibrate transition intensities with the one step transition probabilities estimated from data. This allows us to use non-homogeneous continuous time Markov chains for modeling Long-Term Care. We solve numerically the Kolmogorov forward differential equations in order to obtain continuous time transition probabilities. We assess the quality of the calibration using the Portuguese life expectancies. Based on reasonable monthly costs for each dependence state we compute, by Monte Carlo simulation, trajectories of the Markov chain process and derive relevant information for model validation and premium calculation.


Introduction
Population ageing is a phenomenon which, undoubtedly, is inevitable in the future, in all regions of the world, see (OECD 2019). Forecasts point to a severe ageing of the world population and Portugal is not an exception as pointed out, e.g., in (OECD 2013). This social/economic problem of increasing proportions carries many difficulties to solve, namely the dependence of elderly and convenient provision of Long-Term Care (LTC)that is, the health and well-being support needed in the later stages of life-so it becomes imperative to give special attention to this worldwide problem. This is an issue that actually is on the agenda of many countries, especially the developed ones, where the ageing phenomenon, and consequently the elderly dependence, is more evident. Many countries have already implemented various measures of social protection to elderly dependence, including the USA, Germany, Spain, among others. Some studies have already been published, using real data, such as (Guillen 2008;Fuino and Wagner 2018). In some of these countries the insurance market is already providing LTC products. A substantial work is yet to be done in Portugal. With this paper, besides contributing with a general calibration methodology for transition rates, we aim to provide information on Portuguese data that allows insurers and national care providers to address this problem.
One of the classic approaches to model Long-Term Care Insurance takes advantage of multi-state models as, for instance, the works of (Waters 1984;Haberman and Pitacco 1998;Cordeiro 2002b;Christiansen 2012;Dickson et al. 2013;Fong et al. 2015). Several authors have contributed with important research on this subject. For instance, in (Cordeiro 2002b), approximations to the transition intensities defined for a multi-state model for Permanent Health Insurance are obtained by estimating a generalized linear model for the number of claim inceptions. In (Cordeiro 2002a), using the approximations to the transition intensities and numerical algorithms, which allows for an efficient evaluation of basic probabilities, the average duration of a claim and the claim inception rates by cause of disability are calculated. In (Fleischmann 2015), a method to extract state-change intensities from published empirical observations of prevalence rates is presented. The authors in (Xie et al. 2005) focus on a continuous time Markov model for the length of stay of elderly people moving within and between residential home care and nursing home care. A procedure to determine the structure of the model and to estimate parameters by maximum likelihood is also presented. A Markovian multi-state model in order to calculate premiums for a given LTC Insurance plan is the research object of (Helms et al. 2005), where the approach is a direct estimation of the transition probabilities. In (Kessler 2008), the author presents an analysis of the LTC line of business for insurance companies, an important overview for commercial purposes. An analysis of LTC needs and consequent costs, together with some proposals for the insurer's intervention, is presented in (de Castries 2009). Although dated from a decade ago, in (Mitchell et al. 2008) there is a fairly complete report on the development of LTC in Japan. The subject has also been recently object of significant work, as in (Fong et al. 2015), where generalized linear models are used. In (Rolski et al. 1999, p. 349) we find a construction of the non-homogeneous continuous time Markov chain from an initial distribution and an intensity matrix, allowing for simulation. Further on, in (Rolski et al. 1999, p. 354), there is a description of the pension insurance model with emphasis placed on the obtention of the net prospective premium reserve. In (Haberman and Pitacco 1998, p. 47-61) there is a very important presentation of actuarial values of benefits which can be used to determine the premiums to be paid. A most valuable source on financing issues of Long Term Care Insurance is presented in (Costa and Courbage 2012).
The Portuguese National Network for Continuous Care (RNCCI, in Portuguese), was created in 2006 from a partnership between the Ministry of Health and the Ministry of Labour and Social Solidarity. In this public system, all inhabitants are eligible for LTC care and it provides different forms of care for individuals who need some type of LTC support. For more details on the Portuguese National LTC system, see (Lopes et al. 2018a) and (Lopes et al. 2018b).
At this moment in Portugal there is no private insurance market providing LTC contracts for the general population, mostly due to high risk and lack of data. The data and modelling framework presented in this paper could be used as a starting point for the development of LTC policies in Portugal, since the theoretical framework on these types of insurance is widely studied in the literature.
The primary goal of this paper, therefore, is to develop a multiple state model using official Portuguese data. In (Oliveira et al. 2017), a discrete time transition matrix was estimated from this dataset. In this paper, we aim to relax the assumption of discrete model, allowing to develop a continuous-time Markov model framework, in order to further develop well known pricing and reserving techniques for LTC insurance products, as well as to provide information to the Portuguese National Network for Long-Term Care that allows for cost estimation with this public service.
In this paper, we consider a five state non-homogeneous continuous time Markov model with one autonomous state, three dependent states and one exit state, namely states 1≡autonomous, 2≡light dependence, 3≡moderate dependence, 4≡severe dependence and 5≡death. The choice of the five state model is loosely justified by the widespread use of a reduced Barthel index, see (Mahoney and Barthel 1965), allowing for a general classification of elders in roughly three states of dependence, according to the performance achieved on their daily tasks. In order to obtain a general model, we allow for the transitions between all dependence states, as illustrated in Figure 1. Naturally, our methodology may be applied to other classification schemes and/or a multi-state model with a different set of transitions.
As for the mathematical approach, we follow the Transition Intensities Approach (TIA) which is vehemently advocated in (Waters 1984) and followed by several of the above cited authors. The transition intensities are key parameters in multi-states models as they reflect the dynamics of the underlying processes that shape an individual trajectory through dependence states during elderly ages. All the other corresponding characteristics of the Markov process and the evolution through the LTC model can, generally, be determined if the transition intensities are known, see for instance, (Haberman and Pitacco 1998) and (Waters 1984).
The first goal of this paper is to propose a method of obtaining transition intensities between dependence states by calibration with discrete time transition probabilities estimated from available data. This first goal intends to provide a method for obtaining transition intensities when only discrete time transition probabilities are available, as is the case in this paper. We illustrate the proposed methodology, obtaining transition intensities for the Portuguese population, using data from 2015, of the Portuguese National Network of Continuous Care Database. We believe this will bring valuable information for insurers, practitioners and institutions. For that, we adopt intensities of the form of the hazard function for the Gompertz-Makeham law, as in (Haberman and Pitacco 1998), and we propose a loss function that allows us to calibrate transition intensities from transition probabilities previously estimated from data.
Following TIA, the transition probabilities of the multi-state model are determined from the intensities by numerical integration of the Kolmogorov forward differential equations. The numerical integration is mandatory as, for the model presented in Figure 1, there is no evident closed form solution to the Kolmogorov equations (see Remark 2 ahead). We note that the numerical integration is not a severe limitation to our purposes and methods.
As a second goal we intend to simulate the cost of LTC in Portugal and for that, we use the continuous time transition probabilities derived from the Kolmogorov equations and we use Monte Carlo simulation in order to obtain, for a representative sample of trajectories of the Markov chain, the sojourn times in each state and the associated costs of LTC. These costs, in turn, will allow for the determination of insurance premiums and sojourn times in the system. These sojourn times may give us a a posteriori validation of the methodology by comparing them with the life expectancy of the official mortality table of the Portuguese Institute of Statistics, see (INE 2015), and the healthy life expectancy in Portugal, see (Pordata 2016b). Other authors, see (van den Hout et al. 2014) for instance, have already used multi-state models to predict healthy life expectancy.
The general methodology proposed in this paper will allow for future calibration of intensities if the population's characteristics change or to obtain transition intensities for smaller populations, such as LTC insurance portfolios or even for different schemes of LTC products.

Context and Notations
Let {S(t) , t ≥ 0} denote the continuous time Markov process describing, in each realization, the path of an individual through the different states i of his health condition, with i ∈ S = {1, 2, . . . , r}. For each t ≥ 0, the random variable S(t) takes one of the values in S and the event {S(t) = i} indicates that the individual is in state i at time t.
For states i and j of the continuous time Markov chain and for x, t ≥ 0, we have: 1.
The transition probabilities , for an individual aged x, defined by: i.e., t p ij x is the probability that an individual aged x in state i will be in state j at age that is, the matrix t p ij x i,j=1,...,r is stochastic.

2.
The sojourn probabilities , for an individual aged x, given by i.e., t p ii x is the probability that an individual aged x in state i stays in state i throughout the whole period from age x to age x + t. As a consequence of the Markov property, the transition probabilities satisfy the Chapman-Kolmogorov equations and, also, the sojourn probabilities satisfy In the so called Transition Intensity Approach (TIA) we want to determine the transition probabilities from the transition intensities µ ij x , for i = j, with Remark 1. Let us stress that if, for i = j, we have P ij (x, x + h) ≡ 0 then obviously µ ij x ≡ 0; but, more important for our purposes ahead, we may have µ ij x ≡ 0 and P ij (x, x + h) = 0 for some h > 0; in this case, we only have that the slope at zero of P ij (x, x + h), as a function of h, is null.
We may define the (total) intensity of decrement from state i , i ∈ S, as reflecting the conditional probability of leaving state i over an infinitesimal interval [x, x + dx] given that the process is in state i at time (age) x, by: The Kolmogorov forward ordinary differential equations (with corresponding boundary conditions), for all states i, j ∈ S and 0 ≤ x ≤ t, given by are a consequence of the Chapman-Kolmogorov equations (1). For the sojourn probabilities there is also an equation that can be easily derived by from where we can obtain: In matrix notation, see (Ross 1996), the system of equations (3) may be written as: by considering the intensity matrix For that, just observe that the term with indexes ij , i, j = 1, . . . , r , i = j, corresponding to the matrix product P(x, t)Q(t) is given by: that is, exactly, the Kolmogorov forward Equations (3).

Remark 2.
In the case of constant intensities, the solution of equations (4) are obtained by a direct integration amounting to the computation of a matrix exponential. In the general case, for every x, the Kolmogorov equations (4) are a system of linear equations with an r × r matrix, Q(t), with regular coefficients. For intensities in which we do not have, for every x, y ∈ I, Q(x) · Q(y) = Q(y) · Q(x), we cannot guarantee the existence of a solution obtained by direct integration (see Ver Eecke 1985, p. 146), and that is the case of the model treated in this work. Nevertheless, general results ensure the existence of a unique solution with good properties, see (Teschl 2011) or (Ver Eecke 1985, and by numerical integration, as performed in this work, these solutions can be determined in a form suitable for many other computational purposes.

Modeling Transition Rates
In the absence of a developed market of LTC insurance products, as is the case in Portugal, it may be hard to find quality data allowing the statistical estimation of the parameters of the model, namely in the TIA approach, the intensities parameters. Nevertheless, in the case of having some information on the state transitions allowing for the estimation of a discrete time homogeneous Markov chain model, we will now show that it is possible to calibrate the non-homogeneous continuous-time model to the discrete time one.
Our methodology of modeling transition rates relies on the following idea: (i) Using real data and a multi-state model, estimate discrete time probabilities, which will be used to calibrate continuous-time transition rates; (ii) Adopting a hazard type transition intensities and choosing a given starting age x, derive a numerical solution for the Kolmogorov differential equations, in order to obtain the transition probabilities given by the continuous time multi-state model; (iii) Using an optimization criteria, calibrate the transition intensities parameters in order to meet the n-step discrete time transition probabilities estimated from data; (iv) Finally, using sojourn times of the multi-state model, assess the quality of the calibration with the obtained life expectancy for age x, by comparison with the official life expectancy and by the sojourn time of the autonomous state by comparison with official healthy life expectancy at age x. In our case, these indicators are both available for Portuguese population, see (INE 2015) and (Pordata 2016a).
In our formulation, following (Haberman and Pitacco 1998, p. 25), we consider variable intensity rates of Gompertz-Makeham type, of the form: for some sets Γ, A and B. For each pair of states (i, j) , i = j, we define the set of parameters to calibrate by Although specific data from this regard is lacking, such hypothesis seems reasonable. We note that, in the proposed methodology, other types of transition intensities may be efficiently adopted, such as in (Olivieri and Pitacco 2000), for instance.
We note that assumption (5) is equivalent to consider variable intensity matrices Q(t, θ) in Equation (4), with (non diagonal) entries of the form for some sets Γ, A and B. For Γ, A and B closed bounded intervals we have that the parameter set defined by: being a product of compact real intervals is compact. We will identify the set of all sets of admissible intensities of interest, say I, with Θ. Therefore, to any such element θ of such a set of admissible intensities Θ, there corresponds univocally an intensity matrix Q = Q(t, θ).
Since we only have access to estimated transition probabilities, we consider (5) a reasonable and general assumption. We note, however, that the proposed methodology applies to other types of intensity rates.

Model Calibration via an Optimization Problem
In this section, we will show the existence and uniqueness result of the calibration of transition intensities with data. Theorem 1. Let, for n ≥ 1 and a fixed initial age x n , P x n = p (n) i,j i,j=1,...,r be a sequence of transition matrices for homogeneous Markov LTC models and for any set of intensities µ(t, θ), the correspondent matrix of transition probabilities P(x, t, θ) for the same reference age x. Consider the loss function Then, for the optimization problem inf θ∈Θ O(θ) there exists θ 0 ∈ Θ such that, the unique minimum being attained at possibly several points θ 0 ∈ Θ.
Proof. As we suppose that Θ is compact, we just have to show that, for every fixed t, P(x, t, θ) is a continuous function of θ ∈ Θ. Recall that considering the dependence of the matrix Q on its entries and, of those on the parameters, Equation (3) may be written as an equation that can be read in integral form as: We will now follow the general idea of the proof of the Picard-Lindelöf theorem for proving existence and unicity of solutions of ordinary differential equations as exposed for the backward Kolmogorov equations in (Rolski et al. 1999, pp. 348-49). By replacing P(x, s, θ) in the righthand member of Equation (8) by this righthand member we obtain, P(x, t 2 , θ)Q(t 1 , θ)Q(t 2 , θ)dt 2 dt 1 and, by mathematical induction, we have P(x, t k , θ)Q(t 1 , θ)Q(t 2 , θ) · · · Q(t k , θ)dt k · · · dt 1 . Now, let M(t) = sup θ ij ∈Θ |µ(t, θ ij )| = max θ ij ∈Θ |µ(t, θ ij )| < +∞. As M(t) being continuous is integrable over any compact set, by Lemma 8.4.1 in (Rolski et al. 1999, p. 348), we have that we have that the series which sum represents P(x, t, θ), that is Q(t 1 , θ)Q(t 2 , θ) · · · Q(t n , θ)dt n · · · dt 1 is a series of continuous functions of the parameter θ ∈ Θ converging uniformly and so the sum is a continuous function of the parameter θ.
As expected, having Θ as a compact set, the previous theorem ensures the existence of an optimal calibration under the objective function in Formula (6).
We note that the proposed methodology corresponds to a minimum square error criteria, a common approach in this kind of problems. We will choose the set of intensities that minimizes the loss function (6), thus ensuring a calibration of the intensities to real data.
Remark 3. Despite the fact that the optimization problem defined by Formula (6) has a unique solution, the effective computation of this solution most probably is a hard problem for two different reasons. The first reason is that, in our model, there are three parameters for each intensity and a total of 16 intensities, which amounts to an optimization problem in a space of 48 variables. Such high dimensionality has an implicit a priori challenge (see, for instance (Grosan and Abraham 2009)). The second reason, coupled with the first reason just mentioned, is that the function to be optimized depends on the parameters by the computation of the transition probability functions, using numerical integration. Being so, it seems justified to employ a simple calibration approach to the determination of the parameters, such as the one presented below.
Setting an initial age x n , a time horizon N and considering transition intensities ruled by Equation (3), and considering the model framework represented in Figure 1, we calibrate the transition intensities for each pair of states (i, j) , i = j , i = 1, . . . , 4 , , j = 1, . . . , 5, in order to fit the one step transition probabilities estimated from data.
In a first step, consider baseline parameters for intensities µ 12 x+t and µ 15 x+t , representing the transition intensities from autonomous to light dependence and from autonomous to death, respectively. We propose, as a starting point, the ones assumed in (Haberman and Pitacco 1998, p. 100), namely: In a second step, considering p ij , i = 1, . . . , 4, j = 1, . . . , 5, the one step probability of an individual aged x moving from state i to state j (estimated from data) and γ, α, β the fine tuning parameters chosen to reduce the loss function in Formula (6), we perturb the parameters illustrated in (9) in order to reflect the structure of the one step transition estimated matrix, by the transformations presented in the following:
We note that, from a practical point of view, and for the model framework presented in Figure 1, it means that we consider N × r × (r − 1) = 40 × 5 × 4 = 800 transition probabilities in order to obtain a small value for the loss function O(θ), defined in (6).
The goal at this stage is to have the intensities parameters to reflect the structure of the data inferred matrices. By reducing the complexity of the optimization problem, it becomes more robust and easier to manage.

Outcomes from Portuguese Data
Using data from the Portuguese National Network of Continuous Care, we illustrate the obtained results for the methodology presented in this paper. The database for 2015 contains 23,894 users and a total of 93,134 records of observations. For some more details on the Portuguese LTC database, see again (Lopes et al. 2018a) and (Lopes et al. 2018b).

The Transition Probability Matrix
In this section, we present the discrete time one step probability matrix, obtained from the RNCCI database, by a clustering analysis, as described in (Oliveira et al. 2017), for a set of states and transitions defined in Figure 1. This one step transition matrix will be used in the calibration process for obtaining intensity rates.
Observing the estimated homogeneous matrix (13), for ages over 60 years old, we highlight that: (i) for each dependence state, the higher probabilities correspond to maintaining the same dependence level in one step (one year); (ii) the rates of mortality increase with the dependence level; (iii) we observe non-null transition probabilities between all dependence levels.
For illustration purposes, in Table 1, we present the transition matrices obtained for different sets of ages, estimated from the same database in (Oliveira et al. 2017), and matrix P, already presented in (13), obtained by a convex combination of the other ones, according to the proportion of individuals in RNCCI for each age group.
The above comments over the structure of the one step transition matrix for the whole database remain valid when we focus only on a shorter set of ages. Moreover, analysing the set of transition matrices of Table 1, we can conclude that non-homogeneous transition intensities are more adequate to model the LTC phenomena, since there are significant differences between the magnitude of transition probabilities over different levels of dependence, as age increases. Figure 2 illustrates the mortality rates for each age x , x ≥ 60 for the Portuguese population in 2015, see INE (2016), in comparison with the mortality rates for the population under study on this section. These mortality rates were estimated from matrices on Table 1 and from distribution of RNCCI patients through dependence levels, presented in Table 2.  With P = 0.1969 P 1 + 0.19558 P 2 + 0.17481 P 3 + 0.22991 P 4 + 0.2028 P 5 . Table 2 illustrates the initial distribution of individuals, for the whole population and at age 65, through the dependence states, in the RNCCI database of 2015.

Calibrated Intensity Rates
Setting an initial age x = 65, a time horizon N = 40 and considering transition intensities ruled by equation (3), we calibrated the transition intensities for each pair of states (i, j) , i = j , i = 1, . . . , 4 , , j = 1, . . . , 5, in order to fit the Portuguese transition probabilities illustrated in transition matrix (13).
In Table 3 we present the values of the mean error for each probability corresponding to some choices of the fine tuning parameters (γ, α, β). Taking the set (γ, α, β) = (0.00001, −0.002, −0.000001), as our best choice of parameters, we present, in Table 4, the calibrated parameters for the intensities used in the simulation, which reflects a mean error of 9.68% in the approximation of the transition probability functions by the estimated discrete transition probabilities.  Figure 3 illustrates the transition probability functions for a 65 year old individual, obtained by numerical integration of the Kolmogorov forward differential equations.

Remark 4.
A natural question related to the proposed calibration method is to determine if the relation between the transition matrix used for calibration and the continuous time probability transition matrix obtained by integrating the Kolmogorov equations with the calibrated intensities corresponds-in some sense-to the solution of the imbedding problem for finite Markov chains (see, for instance, (Kingman 1962), (Johansen 1973) and (Johansen 1974)). In the imbedding problem we seek to determine the existence of a continuous time Markov chain transition probability matrix corresponding to a discrete time probability transition matrix that represents the observation of the continuous time Markov chain at unit intervals of time. We observe that in the case where P(x, t, θ) is a solution to the imbedding problem of P = p i,j i,j=1,...,r a one step Markov chain probability transition matrix, then with P (n) = p (n) i,j i,j=1,...,r the n-step probability transition matrix we should have in Formula (6), O(θ) = 0.

Simulation Procedure and Results
After calibrating the transition intensities and integrating the Kolmogorov forward differential equations, we simulate 10,000 trajectories of the continuous time Markov chain in order to obtain significant information over the LTC population in Portugal.
The simulation algorithm used, which is tantamount to a simulation of a continuous time Markov chain, is presented in Algorithm 1. With the simulation process, we obtained information on the sojourn times distribution, illustrated in the left side of Figure 4; the life expectancy at age x = 65, see the right side of Figure 4 and obtained statistical measures for the whole life cost of Long Term Care, by considering a set of reasonable monthly costs for each dependence level, see Figure 5.

Sojourn Times and Life Expectancy
The simulation procedure requires the probability distribution functions of the sojourn in each transient state i , i = 1, . . . , 4, which are obtained by the calibrated intensity rates, by numerical integration of the Kolmogorov forward differential equations. The cumulative distribution functions are displayed in the left side of Figure 4. Table 5 illustrates the average sojourn times for each dependency level, for a 65 year old individual. The "life expectancy" in Table 5 corresponds to the average sum of all sojourn times, ∑ 5 i=1 T ii x , before the individual deceases. We note that the life expectancy given by the official 2015 mortality tables of the National Institute of Statistics in Portugal, see (INE 2015), for a 65 year old PortugPortugueseuese, is 19.19, which reflects that our estimate has a relative difference of −20.74% compared to the official data. According to the recent results on (Bravo et al. 2020), the 95% confidence interval for the life expectancy of Portuguese individuals that, in 2015, were 65 year old is [19.49; 22.66]. In our opinion, the magnitude of this relative difference is explained by the the observation of Figure 2. We note that the population in RNCCI system is mostly a population with comorbidities, which induces higher mortality and, as a consequence, reduces the life expectancy. We remark that we are considering a sub-population of Portuguese individuals that, during 2015, used the National Network of Continuing Care. In fact, as presented in Table 2, 55.22% of the individuals aged 65 were registered as light dependent level and 22.64% as severe dependent. We believe that in the general Portuguese population these percentages are not so high, and this justifies that we observe a lower life expectancy. As a second benchmark, we note that the official Healthy Life years at age 65, in Portugal, in 2015, was 7 years, see (Pordata 2016a). In that sense, we consider our results a reasonable approximation and it gives us a validation for the calibration procedure and results. We also remark that during the simulation procedure, different scenarios have been tested (different ages structure in the population, different initial distributions for the 65 year old individuals) and all of them producing results consistent with the defined scenarios. Our option, in the paper, relied on presenting only the results obtained from the RNCCI database.
We present, in the right side of Figure 4, the empirical distribution for the life expectancy for a 65 year old Portuguese, obtained in the simulation process. We remark that it is possible to obtain similar results for other initial ages, rather than x = 65, given that transition rates µ ij x+t , i = 1, . . . , 4, j = 1, . . . , 5 are calibrated for other ages x as well as for other age structure populations, obtaining a new one step transition matrix from an appropriate convex combination of the matrices presented in Table 1.
We note that, for the population under study, if we had chosen to calibrate the intensities and present the results for the mean age of RNCCI database (80 years old), the differences would be smaller. However, for instance, the Healthy Life years at age 80 is not available for comparison with our results. In that sense, our choice relied on calibrating the intensities for age 65.

Simulating the Costs of LTC
The protection of individuals in later stages of life can be financed by social protection, by insurance policies or a combination of both.
At this point, we intend to obtain, by simulation procedures, the distribution of costs of LTC for a population with the characteristics of RNCCI users, as well as some basic illustrative statistics.
The presented results were achieved by admitting that the monthly cost for users in each of the dependent states is, respectively, 500, 1 500 and 3 000 euros. These values are of the same magnitude of the expenditures per level of care reported in (Zuchandke et al. 2012, p. 221), for LTC in Germany in 2010. Figure 5 illustrates the distribution of LTC costs, under the monthly costs assumed, for a 65 year old individual. The simulated average whole life cost of Long Term Care may be seen as the Risk Premium of an LTC insurance, if the same rates are used for increasing costs and discounting, or as an average amount per patient in social assistance institutions. Regarding the skewness of the distribution, we note that, for insurance purposes, the risk measurement and the premium calculation should reflect this fact. Table 6 highlights the Life Expectancy, the Risk Premium and the correspondent monthly equivalent investment, for a period of 30 years-i.e. starting to invest at the age of 35 years old, with a technical interest rate of 3%-necessary to face the global cost amount. We remark that, naturally, the methodology proposed in this paper allows for the simulation for other populations with different characteristics, namely for the general Portuguese population or for some insurance portfolios. Just for illustration purposes, we show, in Table 7, the advantages of an early subscription over some kind of investment allowing for covering the amount of the Risk Premium of an expensive LTC contract. We illustrate, in line, different subscription ages x and, in column, different initial ages for entering in a LTC system. For each of these pairs (x, x + n) we calculated the monthly investments needed to face the global whole life LTC cost, according to the obtained simulations. The last line illustrates the "life expectancy" after age x + n and the previous one shows the global average cost of LTC after that age. We remark that, as previously referred, the values presented in Table 7 are not perfectly accurate for ages different than 65, since we used, for all ages x + n, the calibrated transition intensities for a generic 65 year old individual. Despite less accurate, the results are still interesting and highlight the need for early investment on LTC coverages.
To conclude this section, we discuss the effect of the number of repetitions in the simulation process. We can appreciate in both Tables 8 and 9 that there is no appreciable difference of the results both on the Risk Premiums computations and on the sojourn times among samples of 10,000 and 10,000 repetitions. As so, our choice of a sample of 10,000 trajectories for the discussion of the simulation results seems justified.

Conclusions
The calibration-simulation method introduced, allows for a study of insurance products and risk management of Long Term Care via a continuous time Markov chain model. The expected life span after age 65 (total life span and healthy life expectancy)-which is an output of the continuous time Markov chain simulation-is a good a posteriori validation tool for the proposed methods when comparing it to data from the official Portuguese mortality tables. We observe that, in our work, we obtained results of the same order of magnitude as those observed in the German LTC system, for similar input values of monthly costs. The calibration procedure proposed is efficient although it does not necessarily ensure a global optimum of the natural loss function. We point out some limitations of our study that opens opportunities for future work. Namely, taking advantage of having data matrices for five age classes would allow for using, in the simulation process, different probability distributions according to the age at a certain point in time of the trajectory being simulated. A second determinant line of work would be to achieve a performant algorithm to the calibration global optimization problem referred in Remark 3.
Some new developments may also arise by considering semi-Markov process modeling, provided the available data to perform such modelling.
Author Contributions: These authors contributed equally to this work. All authors have read and agreed to the published version of the manuscript.
Funding: This work was partially supported through the project of the Centro de Matemática e Aplicações, UID/MAT/00297/2020 financed by the Fundação para a Ciência e a Tecnologia (Portuguese Foundation for Science and Technology). The APC was funded by Future HealthCare.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.