Spatiotemporal Econometrics Models for Old Age Mortality in Europe

In the past decade, panel data models using time-series observations of several geographical units have become popular due to the availability of software able to implement them. The aim of this study is an updated comparison of estimation techniques between the implementations of spatiotemporal panel data models across MATLAB and R softwares in order to fit real mortality data. The case study used concerns the male and female mortality of the aged population of European countries. Mortality is quantified with the Comparative Mortality Figure, which is the most suitable statistic for comparing mortality by sex over space when detailed specific mortality is available for each studied population. The spatial dependence between the 26 European countries and their neighbors during 1995–2012 was confirmed through the Global Moran Index and the spatiotemporal panel data models. For this reason, it can be said that mortality in European population aging not only depends on differences in the health systems, which are subject to national discretion but also on supra-national developments. Finally, we conclude that although both programs seem similar, there are some differences in the estimation of parameters and goodness of fit measures being more reliable MATLAB. These differences have been justified by detailing the advantages and disadvantages of


Introduction
At present, spatial econometric methods are evolving very rapidly at both a theoretical and practical level. Spatial econometrics studies the spatial autocorrelation and spatial structure in regression models for cross-sectional and panel data [1]. These models have been applied in various fields such as sociology, epidemiology, geology, biology, economics, and criminology [2] but less in the field of actuarial science. Nonetheless, its importance in assessing the variation of risk in non-life insurance as a function of a geographical area has been recognized [3]. Further extensions of the spatial modeling related to frequency and severity of claims include models such as [4][5][6].
In the context of life insurance, the disaggregation of overall mortality at some levels, for example, by spatial or socio-economic factors, of mortality measures should be considered in modeling and forecasting according to [7]. In this sense, joint mortality models have been used, such as [8] and references therein.
Even though Debón et al. [9] used some spatial techniques on residuals to forecast dynamic life tables, the importance of the spatial effect has not been well recognized in modeling and forecasting dynamics of mortality. In recent years, there has been research on multi-population mortality models, for example, using a group of countries with similar social-economic status, or males and females in the same population [10]. However, these models do not measure spatial dependence and temporal and spatial effects obtained with data panel models. These methods consider panel regression models with spatially autocorrelated dependent variables and errors with fixed effects (FE) and random effects (RE) specifications.
On the other hand, Debón et al. [11] concluded that there are notable differences in mortality between the Eastern and Western of Europe, a significant disadvantage for Eastern Europe, and especially for males in Baltic countries. Therefore, it is interesting to see if we can explain these differences in mortality through a model including covariables.
For that task, there is extensive literature supporting software that estimates crosssectional spatial models in R and MATLAB; however, few methods estimate panel data spatiotemporal models. Nevertheless, the increased availability of data, and in particular panel data, has stimulated tremendous methodological progress in studies on the specification and estimation of spatiotemporal panel data models. This interest has gone hand-in-hand with developing new packages to model spatial dependence over time in panel data such as [12,13]. At present, the splm R-package [12] is available in R and the [14] code in MATLAB. They both offer access to math functions, language, statistics, and a community of users. Therefore, a quick look at MATLAB [15] and R [16] might suggest that they are relatively similar. However, a closer look at each one's technical capabilities and assessing other vital factors, such as documentation and quality could lead to a different conclusion.
For this reason, this study aims to provide an updated comparison of estimation techniques between the implementations of spatiotemporal panel data models across MATLAB and R software in order to fit real mortality data. The case study used concerns the male and female mortality of the aged population of European countries. Mortality is quantified with the Comparative Mortality Figure, which is the most suitable statistic for comparing mortality by sex over space when detailed specific mortality is available for each studied population [17]. In this context, the use of Generalized Linear Models (GLM) can be traced back to the early 1980s thanks to the work of [18] being suitable for the distribution of the number of deaths.
The article is structured as follows. In Section 2, we describe the data, Comparative Mortality Figure, Global Moran's Index, and panel data models. Next, in Section 3, we present the main results, and models implemented in R and MATLAB. Section 4 provides a discussion about the most important results. Finally, in Section 5 the main conclusions of the study are drawn.

Materials and Methods
This section first describes the structure of the data used and how they were obtained. Secondly, it reviews the methodology employed, mainly relying on data panel models.

Data
The database used comprises deaths and populations of 26 European countries from the Human Mortality Database (HMD) [19] [20] quantified the dynamics of mortality in Europe and detected two significant clusters for ages older than 64. In this paper, to explain the behavior of mortality as a function of socioeconomic variables, information about four variables for these 26 countries and 18 years was obtained from The World Bank Database [21]. These variables are the Gross Domestic Product (GDP), public health expenditure, CO 2 emissions, and education expenditure based on a literature review of this area [22][23][24] and data availability for European countries from [21]. In addition, there are recent studies examining the impact of CO 2 emissions, health expenditure, and education on life expectancy [24]. Furthermore, it is essential to point out that the new impacts of environmental degradation on human health affect today's society both in terms of loss of quality of life and spending on health care. The definitions of the variables can be found below: • Gross Domestic Product per capita (current U.S. dollars) (GDP) is an economic variable used as a measure of income. It is defined as the sum of gross value added by all resident producers in the economy plus any product taxes and minus any subsidies not included in the products' value. GDP has traditionally been used to show the economic and social development of countries [25]. In recent years, there has been significant interest in the relationship between health, proxied by life expectancy, and income, as explained by [26], preceded by [27]. It has generally been well accepted that populations in countries with higher GDP levels have better health and longer life expectancy, as higher living standards lead to enhanced prevention and treatment of disease [27]. It should be noted that all variables are expressed in per capita values. • Public health expenditure per capita (% of total health expenditure) (PHE) is a social variable that consists of recurrent and capital spending from government (central and local) budgets, external borrowings and grants (including donations from international agencies and nongovernmental organizations), and social (or compulsory) health insurance funds. In countries with high income per capita, the contributions to social security are essential and sustain, to a great extent, the financing of the health system. Consequently, the lower the mortality in a country, the healthier its population [28].
Several studies, such as [29,30], among others, have shown that health expenditure has a significant negative impact on mortality rate and a positive impact on life expectancy. It should be noted that the health expenditure variable is reported until 2012; for this reason, the study is limited to that year. • CO 2 emissions per capita (metric tons) is an environmental variable that is used to indicate the effect of air pollution on mortality. Carbon dioxide emissions are those stemming from the burning of fossil fuels and the manufacture of cement. They include carbon dioxide produced during the consumption of solid, liquid, and gas fuels and gas flaring. Countries with higher carbon dioxide emissions levels are at higher risk of their citizens having health problems [24,31]. • Education Expenditure per capita (% Gross National Income) (EE) refers to the current operating expenditure on education, including wages and salaries and excluding capital investments on buildings and equipment. This variable is an essential factor that determines health as a measure of educational level. People with higher educational levels have better jobs, higher incomes, and lower-risk behavior [32].
Gavurova [33] is the closest study to our analysis, but we should highlight three distinctive features in our study. First, we use a suitable statistic to compare mortality for each country by sex, the Comparative Mortality Figure. Second, our sample is contextualized in older European people. Hence, the reliability of our results is enhanced for retired people in Europe. In addition to studying the relationship between mortality and the covariates GDP and health expenditure, we consider the impact of CO 2 emissions and education expenditure per capita over mortality. Third, we consider the neighborhood relationships between countries over time. This point is vital since studies such as Carracedo et al. [20] show that European mortality has a spatial dependence.

Comparative Mortality Figure
The Comparative Mortality Figure (CMF) was used in this study to compare the mortality experience over time by sex and country. There are two principal standardization methods-direct and indirect-which produce Comparative Mortality Figure (CMF) and Standardized Mortality Ratio (SMR), respectively. To compare the mortality experience over time by sex and country, the CMF was used in this study for two reasons. First, the same denominator applies in the calculations for each country, which permits comparison of the mortality experience by sex for different places, and second, the age-specific mortality rates are available for each country [17] required in the expression of CMF but not in the expression of SMR.
The CMF allows for the direct standardization method of mortality rates across age groups to permit comparison of the free effects of differences in the size of subgroups of the populations [17]. In this study, we used males (m) in the set of European countries in the year 2012 as the standard population. Therefore, the CMF is defined as follows, where E i,t,s represents the expected deaths in country i, in year t, and sex s and O 2012,m is the number of observed deaths in the set of European countries in year 2012 for males. As we have d x,i,2012,m the number of observed deaths for the age x, Then E i,t,s are obtained by using the age-specific death rate of each country i for the age x, year t, and sex s, m x,i,t,s , where P x,2012,m is the population in the set of European countries at age x in year 2012 for male sex, which can be obtained by adding all p x,i,t,s the population of each country i with the following expression The death rate of the study population (m x,i,t,s ) are obtained as where d x,i,t,s and p x,i,t,s are the number of deaths and the study population, respectively. If the CMF is greater than 1, there is a higher expected number of deaths than observed; in this case, there are "excess deaths" in the studied countries. On the contrary, a CMF value below 1 indicates a lower expected number of deaths than observed relative to the standard population. According to [17,34], the numbers of observed deaths in the countries is a random variable that follows a Poisson distribution; in the case of CMF that distribution can be used for modelization.

Spatial Dependence of CMF
Global Moran's Index (GM) is widely used to test for the presence of spatial dependence between adjacent locations [35,36]. In this study, it is a summary measure that shows the intensity of spatial dependence of all countries' CMF considered.
where CMF t,s is the mean of the CMF of all countries for each year t and sex s and W N is the neighborhood or spatial weights matrix. The definition of proximity is given by W N , which is an N × N (26 × 26) positive matrix that provides a weight ω i,j to each pair of spatial units (i, j) [37]. In this study, the weights matrix is used in row standardized form.
In addition, two countries are considered neighbors when they have a common border (first order of neighborhood [38]), and therefore, if a country is not a neighbor of itself, it takes the following values, where n i is the number of neighbors of country i and V(i) the set of neighbors of the country i. The interpretation for GM t,s is the following: • GM t,s > 0-Positive spatial autocorrelation between countries. The CMF of countries and their neighbors goes in the same direction. • GM t,s > 0-Negative spatial autocorrelation between countries. The CMF of countries and their neighbors varies in a different direction. • GM t,s = 0-Absence of spatial autocorrelation between the 26 European countries, meaning a random spatial pattern.
Moran's test for spatial autocorrelation was calculated to test the significance of the GM index. The null hypothesis establishes that the CMF i,t,s is randomly distributed among the spatial units of the study area (H 0 : GM t,s = 0) [39].

Spatiotemporal Panel Data Models
The analysis of data from the spatiotemporal panel is currently a field of econometrics undergoing significant methodological advances [12]. Panel data consist of a cross-section of observations (individuals, countries, regions) followed through time. Specifically, the data from this study are panel data that combine a spatial dimension N (26 countries) and temporal dimension T (18 years). In a spatiotemporal panel, there may be dependency or correlation between the close observations (spatial units) over time (temporal units).
Econometrics panel data models offer advantages over cross-section regression or time series as they control for unobserved heterogeneity produced by both spatial and temporal units. As a result, these models reduce the issues related to multicollinearity problems between the variables by building more efficient estimates in the parameters of the panel data models [40]. Panel data usually contain more degrees of freedom and more sample variability, combining both cross-sectional and time-series data, hence, improving the efficiency of econometric estimates [41]. The spatiotemporal panel data models are used as the regression models that employ the panel structure's temporal and spatial heterogeneity to estimate parameters of interest [42].
Due to the nature of our variables, we propose a log-log spatiotemporal panel data model that assumes the log transformation of the CMF and the explanatory variables to achieve approximate normality and symmetry about the distribution of CMF and provide straightforward interpretability of the results [43]. In this model, the explanatory variables' coefficients represent the elasticity of CMF for the explanatory variables [44].
The spatial lag term considers that the value of the CMF in a country depends on the value of the CMF in its vicinity. This fact will be confirmed in Section 3.2. In addition, the fixed effects model is generally more appropriate than the random effects model when the sample used is fixed, i.e., the countries have not been drawn randomly from a very large population [37]. In our case, we would like to model the space-time data of adjacent spatial units (countries) or "located in an unbroken study area" [45]. For these reasons, we fitted spatial lag models with fixed effects. This model fits the behavior of CMF by country and time as a function of explanatory variables assuming that the differences between spatial units, time units, or both are constant [46]. For this reason, spatial and temporal dummy variables were created to account for the unobserved characteristics of cross-sectional units (not changing with time but affecting the dependent variable) and the unobserved characteristics of temporal units (not changing with countries but affecting the dependent variable). These models use the notation that has become usual in this context and adapted to explain the log transformation of CMF, namely, • Spatial Lag Model (SLM): • Spatial Lag Model with time fixed effects (SLMTFE); • Spatial Lag Model with spatial and time fixed effects (SLMSTFE); where log(CMF it ) represents an ordered vector of dimension (26)(18) × 1 corresponding to observations of the dependent variable for each country i and year t, α is the intercept, λ is the spatial parameter associated with the dependent variable, log(X it ) represents a matrix of dimension (26)(18) × 4, the log transformation of the explanatory variables ordered by spatial units first and then by time period, which are related to the parameters β of dimension 4 × 1. The fixed effects considered in the model are as follows: • µ i is the spatial fixed effect (not spatially autocorrelated), which captures the unobservable characteristics that change across countries but remain constant over time. • ν t is the temporal fixed effect (not temporally autocorrelated), which captures the unobservable characteristics that change over time but remain constant across countries.
λ denotes the spatial autoregressive parameter on the spatially lagged dependent variable to follow the econometrics literature [47,48]. SLMSTFE has the conditions that the sum of the spatial and temporal effects are zero [49]. These conditions are achieved using 26 spatial dummies and 18 time dummies because 26 and 18 are the total number of countries and years considered in this study. Thus, the spatial effect represents the deviation of the spatial unit i from the mean α, and the time effect represents the deviation of the time unit t from the mean α.
On the other hand, GLMs are an extension of linear models for response variables with non-normal distributions and nonlinear transformations [18]. GLM provides a method for estimating a function of the dependent variable's mean, also called link function, as a linear combination of a set of explanatory variables. Poisson regression is a GLM model where the dependent variable is a count that follows a Poisson distribution. The canonical link results in a log-linear relationship between mean and linear predictor. The Poisson model variance is identical to the mean; thus, in the case the variance is larger than the mean, the data are over-dispersed.
A way of modeling over-dispersed count data is to use the quasi-Poisson model, which produces the same coefficient estimates as the standard Poisson model, but the inference is adjusted over-dispersion. Consequently, the quasi-Poisson family does not correspond to a model with fully specified likelihood, and therefore, statistical tests and goodness-of-fit measures such as AIC, BIC, likelihood-ratio are unavailable in the output of model [16]. In R, the quasi-Poisson model with estimated dispersion parameter can also be fitted with the glm function and setting family = quasipoisson.
We suggest adapting log-log SLMSTFE to GLM, the model proposed for log(CMF) is, In the next section, these models are estimated by maximum likelihood with splm and glm functions in R software and the sar_panel_FE function in MATLAB. In this way, the glm function takes into account the value of the quasiPoisson-likelihood and not the log-likelihood value.

Exploratory Analysis of the CMF
A preliminary study of the temporal evolution of the CMF of European countries for males and females is discussed in this section. Firstly, the CMF of men and women in Europe from 1995 and 2012 was represented on the maps. These maps are all comparable to each other as the color scale is the same for all these figures using the corresponding common quartiles of all years and both sexes. Figure 1 shows the CMF in Europe for the years 1995 and 2012 for males. It should be noted that the CMF of men for the year 1995 in Europe has a value higher than 1.12, which is exceptionally high for the eastern countries, followed by Ireland, Denmark and, finally, the rest of the western ones. Its evolution for 2012 shows lower values than in 1995, with a minimum of 0.462 reached in the western countries, followed by some eastern countries and incredibly high for Belarus.   The same analysis for females is shown in Figure 2. The value of CMF for 1995 in Europe ranges between 0.462 and 1.45, lower than in males. In 2012, the female CMF decreased so that it reached at most 1.12 for the second quintile of the scale in eastern countries.   Secondly, CMF variability was explored for both sexes versus the countries and the years. Figure 3 shows the comparison of the CMF range for males and females in all countries. The countries were ordered alphabetically, and a higher average value and variability were confirmed for men, especially in the eastern countries. Figure 4 shows the comparison of the CMF range for males and females in all years. Again, higher variability in men than in women was detected, while a decreasing trend was observed in both sexes.   After the exploratory analysis of the CMF, the next step was to test the significance of mortality dependence on space where the null hypothesis is no spatial correlation.

Spatial Dependence of CMF
The Moran tests for each year are shown in Table 1 for males and females. The p-values obtained for all years are significant (p-values < 0.05), indicating a spatial dependence between the mortality of the 26 countries in the period 1995-2012 in both sexes.

Spatiotemporal Panel Data Models
After testing the significance of mortality dependence, the next step was fitting spatiotemporal panel data models.

Fitting Spatiotemporal Panel Data Models
Parameter estimation for the spatiotemporal panel data models described in Section 2.4 were performed using the splm and glm functions of splm and stats R-packages, respectively. It is important to note that the version of the splm R-package used was 1.5-2, which produces the same parameter estimates as version 1.3-7. At present, versions other than the two mentioned are not correct because model fitting does not work well. Regarding MATLAB software, the sar_panel_FE function of demopanelscompare file was used.
It is essential to mention that the dataset's structure for fitting the spatiotemporal panel data models is different depending on the software used. In R, both the dependent and the independent variables are ordered by spatial units first and then by time. On the other hand, in MATLAB, all variables are ordered by time first, then by spatial units.
A summary of the results of the spatiotemporal panel data models together with the corresponding goodness-of-fit measures by sex are shown in Table 2 for splm function in R, Table 3 for sar_panel_FE function in MATLAB and Table 4 using the glm function in R. Parameter estimations with p-values less than 5% were considered significant.   The SLMSTFE is the one with the highest R 2 and likelihood values and lowest value of σ 2 in R and MATLAB. However, the spatial effect in the model has more weight than the temporal, and the inclusion of the time effect increases both goodness-of-fit measures, but not significantly. For this reason, and following the principle of parsimony, SLMSFE is the best model for both sexes.
The SLMSTFE has the lowest residual deviance value, but, as before, the SLMSFE is the best model for both sexes. In this case, only the education expenditure variable is significant for SLMSFE in males with the glm function, which simplifies the model.
Next, the results of the estimation of the spatial fixed effects of SLMSFE with MATLAB are shown in Table 5. This summary includes the model's fixed spatial parameters, µ i , the t-value of these parameters, and the corresponding p-values. All spatial effects µ i are not significant, but the global spatial effect is significant according with Table 3 where the inclusion of µ i in SLM model increases the explanatory power of the model.

Validation Spatial Lag Model with Spatial Fixed Effects
The SLMSFE model was assessed using the R 2 , the value of the log-likelihood, and σ 2 in MATLAB. The R 2 for males and females are 98.60% and 96.40%, respectively. The values of the likelihood for males and females are 904.11 and 735.63, respectively. Considering σ 2 , its value for males is 0.0008 and females 0.0023. Since most of the variability in these models is explained, we conclude that these two models fit well. It is important to note that the corresponding R 2 for the ordinary least squares models by sex is 84.02% for males and 79.46% for females. These results indicate the importance of modeling the spatial correlation since the inclusion of the model's spatial effect improves its fit by approximately 15% and 17% in males and females.
Next, residuals were analyzed to check if the SLMSFE model explains the whole spatial dependence of the 26 European countries detected by the Global Moran Index (see Section 2.3). Table 6 shows the result of applying the Global Moran test to the SLMSFE model residuals in 1995-2012 by sex. The p-values obtained for all years and both sexes are not significant (p-values > 0.05), indicating that non-spatial dependence remains in the residuals. Then, in order to check that the residuals were homoscedastic in the SLMSFE model, the Breush Pagan test [50] was applied using the function bptest of lmtest package [51]. Table 7 shows the output of Breusch-Pagan statistics and the corresponding p-values by sex. It can be observed that the p-values in the residuals of the SLMSFE model are greater than the 5% significance level for both sexes. Therefore, there was no evidence to conclude that the SLMSFE model residuals are not homoscedastic. On the other hand, the SLMSFE model residuals have an average around zero; specifically, for males and females they are 1.232 × 10 −17 and −4.285 × 10 −17 , respectively.
Next, the normality of the SLMSFE model residuals by sex was analyzed with the function qqPlot of car R-package. Figure 5 represents the quantile-quantile (Q-Q) plot of the SLMSFE model residuals, which suggests that the error terms are normally distributed for both sexes. It should be noted that, for females, two high leverage observations were detected: observation 417, which corresponds to Slovakia in the year 1997, and observation 454, which corresponds to the UK in the year 1998.  Tables 2 and 3 show very similar parameter estimations: only the second or third decimal is different between coefficients. Nevertheless, different goodness-of-fit measures were obtained in the output of the models. While the splm function comprises only the value of log-likelihood (logLik) setting quiet = FALSE, the sar_panel_FE function contains three measures by default: the value of the coefficient of determination (R 2 ), the value of the log-likelihood (log-likelihood) and the residual variance (σ 2 ). The log-likelihood obtained within MATLAB versus R is more reliable as some values are negative for the latter one. It is noteworthy that the SLM model is not obtained in Table 2 because the non-fixed effect option is not possible in the argument effect of the splm function. Table 4 shows that the values of the parameter estimation are different from those in Tables 2 and 3. The differences in sign are marked in bold, which correspond to the variables log(CO 2 ), log(EE), and log(PHE), which are almost always not significant. The signs of the spatial parameter and log(GDP) variable remain constant using the three functions. This is due to the fact that glm assumes Poisson distribution, while the splm and sar_panel_FE assumes Normal distribution in the maximum likelihood estimation. It should be noted that there is only one measure of goodness-of-fit in the output glm, which is the residual deviation.
To incorporate the fixed effects in the function glm a reference level is set; by default, it is the first value in alphabetical order. Therefore, α represents the fixed effect of that country for the first year. Moreover, this reference can be changed thanks to the function contr.treatment. On the contrary, the functions splm and sar_panel_FE consider α as the mean because they have the conditions that the sum of the spatial and the sum of the temporal effects are zero [49]. Thus, the spatial effect represents the deviation of the spatial unit i from the mean α, and the temporal effect represents the deviation of the temporal unit t from the mean α.
Despite there being differences in parameter estimation and goodness-of-fit measurements using functions splm, sar_panel_FE, and glm, SLMSFE is the best for all of them. The results suggest that with sufficient social and economic variables for each country and considering the neighborhood structure, we can adequately explain any geographical structures in European mortality of old age.

Interpretation Spatial Lag Model with Spatial Fixed Effects
In this subsection, the SLMSFE is interpreted as being the best model using both pieces of software and the three functions. The Elhorst program code [14] in MATLAB was selected as it is considered a pioneer in including spatial dependence in panel data, and millions of engineers and scientists have tested its algorithms [52].
In the SLMSFE, all the independent variables are significant for females, with p-values less than 0.05 indicating that they all contribute significant information in explaining log(CMF); in the case of males, log(PHE) is not significant. Moreover, coefficients represent the percentage change in Y when the percent change in X is 1% as Equation (4) is a log-log model.
The estimate for α represents the average value of log(CMF) when the fixed effects and explanatory variables take values equal to 0. This value is positive for males (1.479) and females (2.901). The spatial lag estimations (λ) are highly significant in both models, supporting the conjecture that the value of log(CMF) in one country is influenced by the values of the log(CMF) in the neighboring countries. For males, the positive values of λ indicate that the CMF in one country would be 0.806% higher if its neighbors had a CMF average increase of 1%. In females, the λ positive value is lower than in males; the CMF in one country would be increased by 0.524% if its neighbors had a CMF average increase of 1%. These results indicate a greater spatial dependence on male than female mortality for the aging population. Therefore, mortality in European population aging depends not only on differences on the health systems but also on supra-national developments [20].
The log(GDP) coefficient indicates that with a 1% increase in GDP, the elasticity of CMF is −0.040% for males and −0.103% for females, on average. This relationship is confirmed in several studies such as the recent study by Gavurova et al. [33], among others. These authors studied the influence of socioeconomic determinants on mortality by sex in the European Union (EU) using panel data models. They confirm that GDP per capita, representing countries' economic growth, decreases mortality in EU countries for both sexes.
With respect to log(CO 2 ) emissions, the coefficient indicates that, with a 1% increase in CO 2 , the elasticity of CMF is 0.057% for males and 0.123% females, on average. This result is consistent with Balan [24] among others, who studies whether there is a causal relationship between environment and health, measured by life expectancy, for EU countries. The author considers CO 2 emissions as the logical consequence of polluting industrial activities, which implies strong economic growth. The authors of [24] confirms that CO 2 emissions decrease life expectancy, thereby increasing mortality.
Finally, concerning the coefficients of two social variables log(EE) and log(PHE) following the literature such as [24,33], a negative relationship would be expected between overall mortality and both variables. In our study, these coefficients are positive, and in the case of log(PHE) for males are not significant. There is controversy about the positive sign of both variables, which would require further investigation. However, we think that this occurs as we are studying old-age mortality.
With respect to spatial fixed effects in Table 5, it is worth mentioning that although the countries do not differ from the average α at least, the extremes differ between them. For both sexes, estimates with a negative sign indicate countries with lower log(CMF) than the mean α; on the contrary, the estimations of spatial effects with a positive sign belong to countries with higher log(CMF). This sign means that these countries' unobserved characteristics decrease or increase mortality compared to the all countries' average.

Limitations
In this subsection, it is worth mentioning two critical limitations of this study. The first is the database used, which could be expanded with other socio-economics, demographic, environmental, and health variables from The World Bank Database [21]. These variables should be related to mortality and are perfectly justified with previous studies. This justification is fundamental to avoid irrelevant variables and simplify interpretability in the independent variables. On the other hand, the choice of the World Bank Database is crucial in this spatial study because it allows us to have more information on neighboring spatial units for a longer time horizon than other databases such as Eurostat.
The second limitation is that our objective is specific to mortality data. We fit spatiotemporal panel data models to real mortality data. It is not general, fitting to any data. Studying mortality has always been of interest in the actuarial and epidemiological context. However, the coronavirus disease 2019 (COVID-19) has produced unknown effects on the overall population mortality [53]. Therefore, now more than ever, mortality is attracting considerable research attention, even more so in advanced ages.

Conclusions
This study reproduces the sar_panel_FE function written by Elhorst [14] in MATLAB and in free software R using the splm and glm functions. R software is free and accessible to everyone, while MATLAB is commercial software with a very high price. Nevertheless, millions of engineers and scientists have tested the spatial-temporal algorithms of MATLAB; on the contrary, some functions such as splm in the splm R-package are more recent.
Therefore, this study compares MATLAB and R implementations of spatiotemporal econometrics panel data models: 1.
The main difference lies in the second or third decimal of estimated parameters.

2.
Signs and values estimated parameters differ when using the glm function in R. This is because the glm function assumes Poisson distribution, while the splm and sar_panel_FE take on Normal distribution in the maximum likelihood estimation.

3.
Concerning the splm function, at the present time, we have found that versions other than 1.3-7 and 1.5-2 do not correctly estimate the model parameters.

4.
The goodness-of-fit measures in the output are different depending on the function used. The splm function only gives the value of the log-likelihood, the sar_panel_FE function offers the values of coefficient of determination, log-likelihood, and residual variance, and the glm function the value of residual deviance.

5.
The log-likelihood values obtained in MATLAB are more reliable than in R because negative values appear in the latter. The tree function showed that the SLMSFE is the model that best fits the European old-age mortality data where the spatial effect is essential and the temporal one does not appear. 6.
An important advantage of the glm function compared to the rest is that the reference level for fixed effects can be changed. On the contrary, the splm and sar_panel_FE functions consider the reference level for fixed effects as the mean of fixed effects.
Considering other authors' works, we should highlight the methodology's distinctive features presented here. As far as we know, the spatiotemporal panel data models implemented by the splm R-package [12] and MATLAB [14] have not been used to study the CMF of European countries. The spatiotemporal panel data models used by [54] model mortality quantified by the Standardized Mortality Ratio (SMR) in the United States. We must point out that the CMF reflects better than the SMR the changing situation of mortality in the countries as its expression requires the age-specific mortality rates for each country. Our work could be complementary to [55], since we show the advantages and disadvantages of each function used in R and MATLAB software. Bivand and Piras [55] constitutes an exciting comparison of a generalized method of moments and maximum likelihood implementations for spatial econometrics models using MATLAB, Stata, Python, and R.
This article shows that mortality modeling should take into account both the socioeconomic characteristics underlying the modeling process and the neighborhood relationships of geographic locations. Therefore, having good, quality, socioeconomic data and correctly setting the neighborhood criterion data panel models can quantify the vast majority of the variation in European old-age mortality by country. In the changing environment in which we live, it is necessary for the authorities responsible for population planning to quantify those changes that transcend the borders of any one spatial unit. Therefore, it is of vital importance to study the mortality in Europe, especially in advanced ages. Such understanding aims to ensure the correct formulation and implementation of policies for sustainable development in Europe. In this sense, we would like to point out that although this study only applies to mortality for European countries, this methodology can be extended to the comparison of mortality of subpopulations in any geographic area when age-specific mortality rates are available for each spatial unit and a panel of covariates with the same spatial and temporal dimension as mortality.

Funding:
The authors are grateful for financial support through a grant from Mapfre (ayudas a la investigación Ignacio H. de Larramendi 2017 Seguro y Previsión Social).

Acknowledgments:
We would like to express our sincere thanks to the anonymous reviewers for their careful review of the manuscript and valuable remarks.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: