This paper analyzes yearly death counts after World War II of 8 countries in three regions, North America, Europe, and Asia-Pacific region to identify and to compare stochastic structures of death counts. The 8 countries are the United States and Canada in North America, United Kingdom, France, Italy, and Spain in Europe, and Taiwan and Australia in Asia-Pacific region. Death counts are from the year 1946 or 1970 (depending on availability) and are separated by gender and total counts to see whether gender influences the stochastic structures. The structural time series models (Harvey, 1981, 1989) assume that a time series can be formulated directly with the unobserved components such as trend, slope, seasonal, cycle, and daily effect. The random effect of each unobserved component is characterized by its stochastic structure and a distribution of its irregular component. Structural time series models that this paper entertained are a local level with a random walk model, a fixed local linear trend model, and a local linear trend model. These models are sensible choices based on the preliminary examination of death counts data. Structural time series models use the Kalman filter (Kalman, 1960) to estimate unknown parameters of the entertained model, forecast future values of time series, and estimate unobserved components in the stochastic model by filtering and smoothing. To apply the Kalman filter, the structural time series models need to be converted to a state space form (Durbin and Koopman, 2012) that is the standardized form for the Kalman filter. To check the validity, two diagnostic procedures are used for a fitted model: One is for checking the normality of residuals, for which the Shapiro-Wilk test is used and the normal QQ plot and density plot are used to confirm the results of the Shapiro-Wilk test. Second is for checking the independence of residuals, for which the Run test is used. To find the best-fitted model among valid models, Akaike information criterion (AIC), Bayesian information citerion (BIC) and sum of square of one-step-ahead prediction errors (SSPE) are used. The best-fitted valid models of death counts for each country by female, male and total are fully examined to see any differences or similarities among countries and regions. The organization of this paper is as follows. In Section 2, the structural time series model is presented. It also shows how to set up the state space form for each structural time series model entertained. In Section 3, the Kalman filter is introduced. This section also shows which R packages to use for the Kalman filter. In Section 4, results of analyzing death counts of 8 countries by female, male and total are presented. Finally, Section 5 concludes the paper.

The structural time series models assume that a time series can be formulated directly with the unobserved components that are characterized by its own stochastic structure and an irregular term. By varying its own stochastic structure and the distribution of an irregular term of an unobserved component, structural times series models can fit a variety of time series in many fields. Outcomes of a structural time series model are estimates of unknown parameters, forecasts of future values of time series, and estimates of unobserved components in the model. Estimates of unobserved components of a model give an insight to fully understand a stochastic structure of a time series of interest. There are many papers regarding the structural time series model. For example, Harvey and Todd (1983) compared the structural time series model with Box and Jenkins’ ARIMA model. Harvey and Peters (1990) showed the number of methods to compute the maximum likelihood estimator of unknown parameters of the structural time series model.

For death counts data, three structural time series models are entertained: local level with random walk model, fixed local linear trend model, and local linear trend model. Choices of these models are based on the examinations of plots of death counts of 8 countries by female, male and total. Plots of death counts of all 8 countries by female, male and total are presented in Section 4.

2.1. Local level with a random salk model

First structural time series model entertained in this paper is the local level with a random walk (LLRW) model:

where y_t is an observed time series data, μ_t is an unobserved trend component that represents the long-term movement in a time series, η_t is an irregular component that shows the stochastic behavior of the trend of time series, and ε_t is an irregular component that shows the stochastic behavior other than the trend defined in the second equation in (2.1).

2.2. Local linear trend model

Second model entertained is the local linear trend (LT) model:

where y_t is an observed time series data, μ_t is an unobserved trend component, β_t is an unobserved slope component, ζ_t is an irregular component that shows the stochastic behavior of the slope of the trend of time series. η_t and ε_t are defined as (2.1). This model shows that not only a trend of series but also a slope of the trend is stochastic with its own irregular component.

2.3. Fixed local linear trend model

Third model is the fixed local linear trend (FT) model that is a variation of LT model:

where y_t is an observed time series data, β is a deterministic slope of the trend, μ_t, η_t, and ε_t are defined as (2.1) and (2.2). This model assumes that a series has a deterministic slope, rather than a stochastic.

2.4. State space form

Structural time series model uses the Kalman filter to estimate unknown parameters, to forecast future values and to estimate unobserved components. To use the Kalman filter, a state space form (Durbin and Koopman, 2012) is required. The state space form is a standardized form of a stochastic model as an input for the Kalman filter. For the LLRW model (2.1), it is the state space form itself. Thus, the state space form (2.1) has 1 × 1 state vector, μ_t and two irregular components, η_t and ε_t. The state space form of the LT model (2.2) is:

The state space form (2.4) has 2 × 1 state vector, (μtβt) and three irregular components, η_t, ζ_t and ε_t. The state space form of the FT model (2.3) is:

The state space form (2.5) has 2 × 1 state vector, (μtβt) and two irregular components, η_t and ε_t.

Kalman filter has been applied to many fields since Kalman’s first paper (Kalman, 1960) was published. Harrison and Stevens (1971) was the first paper to apply the Kalman filter to a time series analysis. Since then, the Kalman filter has been used to analyze several time series models such as ARIMA models (Box and Jenkins, 1976), structural time series models (Harvey, 1989), and ARMAX models (Hannan and Deistler, 1988). Some areas where it has been applied include disease control (Gove and Houston, 1996), actuary claim reserves forecasting (Chukhrova and Johannssen, 2017), rain fall forecasting (Asemota et al., 2016; Zulfi et al., 2018), and machine learning (Nobrega and Oliveira, 2019).

Kalman filter can be applied to either univariate and multivariate time series, and to either time variant structure or time invariant structure of time series. Time series data, Y_t, for t = 1, 2, …, T denotes the observed values of a time series of interest. Y_t could be a univariate or multivariate. α_t denotes the unobserved component vector, called the state vector. For the Kalman filter, the observed data Y_t and the unobserved state vector α_t have the linear relationship as

where Y_t is N × 1 (N = 1 for the univariate and N > 1 for the multivariate), Z_t is N × m matrix, α_t is m × 1 vector, and ε_t is N × 1 vector of irregular components. For equation (3.1) , Z_t is assumed to be a known quantity that shows the relationship between Y_t and α_t, and ε_t is assumed to have a multivariate normal distribution with mean of N × 1 vector of zero and covariance of N × N matrix, h_t. Equation (3.1) is called the observation equation. For a univariate Y_t, the observation equation (3.1) can be written as

where Y_t is a scalar, zt′ is 1 × m vector, α_t is m × 1 vector, ε_t is a scalar irregular component whose distribution is a normal with mean of zero and a scalar variance, h_t. The other equation required in the Kalman filter is called the transition equation, which shows the stochastic behavior of the state vector, α_t. The transition equation is also linear as

where T_t is m × m matrix, R_t is m × g matrix, γ_t is g × 1 vector of irregular terms whose distribution is assumed to be a multivariate normal with mean of g × 1 vector of zero and covariance of g × g matrix, Q_t. The matrices Z_t, T_t and R_t, and covariance matrices, h_t and Q_t may or may not change over time. If the matrices do not change over time, the model is called the time-invariant Kalman filter. Otherwise, it is called the time-variant Kalman filter. Equation (3.1) and (3.3) together is called the multivariate state space form of the Kalman filter, and (3.2) and (3.3) together is called the univariate state space form of the Kalman filter. By comparing (3.2) and (3.3), the state space form of the LLRW model (2.1) has the following identities:

The state space form of the LT model (2.4) has the following identities:

The state space form of the FT model (2.5) has the following identities:

The Kalman filter has three assumptions: A1) the initial state vector, α₀ has a mean of a₀ and Covariance matrix of P₀. A2) the disturbances η_t and γ_t are independent of each other. This assumption could be relaxed. A3) the disturbances η_t and γ_t are independent with the initial state, α₀.

Given the information of initial state variable, α₀, the Kalman filter starts off the recursive algorithm. The recursive algorithm provides two estimates of the state variable α_t. First is the filtered estimate. It is the estimate of α_t given Y_t = (y_t, …, y₁). Second is the forecast estimate. It is the estimate of α_t for t = T + 1, T + 2, … given Y_T = (y_T, …, y₁). The filtered estimates of the state variable α_t provide estimates of unobserved components in the structural time series models. These estimates are minimum mean square estimators (MMSE).

The information required to start off the Kalman filter is mean and covariance matrix of the initial state vector, α₀, which are a₀ and P₀, respectively. If the state vector α_t is a nonstationary, then a distribution of α₀ is given as a diffuse prior whose mean is given as the first data in the series and covariance matrix is given as kI where k is a very large scalar and I is an m × m identity matrix where m is the dimension of the state vector, α_t. For the LLRW model, α_t = μ_t and it is a nonstationary and m = 1. Therefore, a₀ is set as y₁ and P₀ is set as 1.2 e + 10. For the LT model, αt=(μtβt) is also a nonstationary and m = 2. Therefore, a₀ is set as (y10) and P₀ is set as (k00k) where k = 1.2e + 10. For the FT model, αt=(μtβt) is also a nonstationary and m = 2 and slope β_t is a deterministic. Therefore, a₀ is set as (y1β) and P₀ is set as (k000), where k = 1.2e + 10 and the value of β is the estimate of slope obtained by doing a simple linear regression of yearly data on time.

In the state space form in (3.1) and (3.3) for a multivariate time series or in (3.2) and (3.3) for a univariate time series, there are some unknown parameters in matrices, Z_t, T_t, h_t and Q_t. Before running the recursive algorithm, these unknown parameters should be estimated. From (3.4), (3.5) and (3.6), unknown parameters of LLRW, LT and FT models are ψ_LLRW = (VarN, VarE), ψ_LT = (VarN, VarK, VarE) and ψ_FT = (VarN, VarE), respectively. If the disturbances of ε_t and η_t are normally distributed, the likelihood function of the observations could be obtained from the Kalman filter via the prediction error decomposition (Harvey and Peters, 1990). Unknown parameters are estimated by maximizing the likelihood function with respect to the unknown parameters.

To analyze our data by the Kalman filter, a R function, StructTS in the R package, stats is used. The function StructTS, first, estimates unknown parameters in matrices, Z_t, T_t, h_t and Q_t, and then given the estimates of unknown parameters, provides predictions and filtered estimations for univariate time series based on the state space form of (3.2) and (3.3). To estimate unknown parameters in matrices, Z_t, T_t, h_t and Q_t, StructTS calls a R function, optim in the R package, stats. For the method for optimization in optim, L-BFGS-B is used. L-BFGS-B is described by Byrd et al. (1995) which allows box constraints, which means that each variable can be given a lower and/or upper bound. For our models, since unknown parameters are variances of irregular terms, the lower bound is set as 0 and the upper bound is set as infinity. The initial values to start the function optim are any nonnegative numbers. This method uses a limited-memory modification of the BFGS, which is a quasi-Newton method. Outputs of StructTS are estimates of unknown parameters, loglikelihood, standardized residuals, and filtered estimates of state vector α_t. The standardized residuals are used for the diagnostics of fitted models. From the filtered estimates of α_t, estimates of unobserved components of LLRW, LT and FT are obtained.

4.1. Death counts data

Yearly death counts data of 8 countries analyzed in this paper are extracted from the Human Mortality Database (HMD). HMD provides several data such as death counts, census counts, birth counts, and population estimates for calculations of death rates and life tables. The main goal of the HMD is to document the longevity revolution of the modern era and to facilitate research in its causes and consequences. HMD includes relatively wealthy and highly industrialized countries since it is based on design to populations where death registration and census data are virtually complete. In this paper, death counts data of 8 countries are analyzed: Two countries in North America (U.S. and Canada), four countries in Europe (U.K., France, Italy, and Spain), two countries in Asia-Pacific (Taiwan and Australia). From death counts from HMD, three death counts data are generated for each country: Female death counts, Male death counts, and Total death counts. The annual periods of data for countries are for U.S. (1946–2017), Canada (1946–2016), U.K. (1946–2016), France (1946–2017), Spain (1946–2018), Italy (1946–2017), Australia (1946–2018), and Taiwan (1970–2014). HMD keeps death counts data for Taiwan only from the year 1970. This paper uses data after World War II (1939–1945) since there are serious outliers during World War II, especially for data of countries in Europe.

4.2. Diagnostics and Best fitted model selection

All three models, LLRW, FT, and LT are fitted to three types of death counts (female, male and total) of all 8 countries. For each model fitted, two diagnostic procedures based on residuals are applied to see the validity of the fitted model. The residuals in the Kalman filter are obtained from one-step-ahead-prediction errors (also called innovations). One-step-ahead-prediction errors are obtained by the Kalman filter with the estimated values of unknown parameters of the fitted model. Standardized residuals are obtained from one-step-ahead-prediction errors divided by the standard deviation of one-step-ahead-prediction errors. If a fitted model is valid, then the standardized residuals are independent and identically distributed by a normal distribution (Harvey, 1989). Thus, the first diagnostic procedure is to check the normality of the standardized residuals and the second diagnostic procedure is to check the independence of the standardized residuals. Shapiro-Wilk test is used to check the normality, and normal QQ plots and density plots are used to confirm the conclusions of Shapiro-Wilk test. Run test is used to check the independence of the standardized residuals. Models that pass two diagnostics procedures are treated as the valid ones. Once the valid models are identified for each type of data after the diagnostic procedures, a best-fitted valid model is selected based on AIC (Akaike Information Criterion), BIC (Bayesian Information Criterion) and SSPE (Sum of Square of one-step-ahead-prediction errors). If these three criteria do not recommend a model unanimously, then the most recommended one is selected as the best-fitted valid model.

4.3. North American countries, U.S. and Canada

Figure 1 shows plots of death counts of female, male, and total for U.S. and Canada. Both U.S. and Canadian data show a similar linear trend for death counts of female, male and total. Also, both countries’ data show that death counts for males are larger than those of females up to the year 2000. Since then, death counts of females and males are close to each other. Since both countries’ data show a clear linear trend, it is sensible to fit the LT and FT model. The LT model assumes that the slope of linear trend is stochastic, and the FT model assumes that it is deterministic. Outcomes of the LLRW model are also provided to compare with them. Table 1 shows p-values of Shapiro-Wilk tests for testing normality of the standardized residuals of three models fitted to three types of death counts in U.S. and Canada. The p-values of two countries show that none of the three models failed for normality on the level of 0.01. Figure 2 and 3 show normal QQ plots for U.S. and Canada, respectively. These plots confirm the results of Shapiro-Wilk tests. Table 2 shows the p-values of Run tests for independence of standardized residuals. The p-values of two countries show that none of the three models failed for independence on the level of 0.01. Thus, for two countries, none of the three models failed to be valid. Table 3 presents the best model for each type of data based on AIC, BIC and SSPE for the valid models of U.S. and Canada. If a model is not selected unanimously by all three criteria, then the most commonly recommended one is selected as the best one. For example, for total deaths of U.S., the FT is the selected as the best one. Thus, for the U.S., the FT models are the best for both female and total deaths and the LT model is the best for male deaths. For Canada, the LT models are the best for both female and total deaths and the FT model is the best for male deaths. Table 4 presents the estimated values of unknown parameters of the best model for each type of deaths in U.S. and Canada. It is noted from Section 3 that unknown parameters of LLRW, LT and FT model are, respectively, ψ_LLRW= (VarN, VarE), ψ_LT = (VarN, VarK, VarE) and ψ_FT = (VarN, VarE). Thus, for FT and LLRW models, the parameter estimate for VarK is NA.

4.4. European countries, U.K. and France

Figure 4 and 7 show plots of death counts of female, male and total for U.K. and France, and Italy and Spain, respectively. Deaths of European countries do not show a clear linear trend that deaths of North American countries show. Table 5 shows p-values of the Shapiro-Wilk tests for testing normality of the standardized residuals of models fitted to three types of death counts in U.K. and and France. The p-values of two countries show that none of the three models failed for normality on the level of 0.01. Figure 5 and 6 show normal QQ plots of U.K. and France, respectively. These plots confirm the results of Shapiro-Wilk tests. Table 6 shows the the p-values of Run tests for independence of standardized residuals. The p-values of two countries show that none of the three models failed for independence on the level of 0.01. Thus, for two countries, none of the three models failed to be valid. Table 7 presents the best model for each type of data based on AIC, BIC and SSPE for the valid models for each type for the U.K. and France. As North American countries, if a model is not selected by unanimously, then the most commonly recommended one is selected. For the U.K., the LT models are the best for all three types and for France, the LLRW models are the best for all three types. Table 8 shows estimates of unknown parameters of the best model for each type of deaths in U.K. and France. As Table 4, for the LLRW models, parameter estimates for VarK are NAs.

4.5. European countries, Italy and Spain

Table 9 shows p-values of Shapiro-Wilk tests for testing normality of the standardized residuals of models fitted to three types of death counts in Italy and Spain. The p-values show that for Italy, none of models failed for normality on the level of 0.01, and for Spain, however, the LLRW models for both male and total deaths failed for normality on the level of 0.01. Figure 8 and 9 show the normal QQ plots for Italy and Spain, respectively. For Spain, normal QQ plots of both male and total deaths of the LLRW model show clearly the non-normality (left skewed shape) of data which confirms the results of Shapiro-Wilk tests. Figure 10 shows the normal QQ plots of both male LLRW and total LLRW model of Spain again and corresponding density plots for standardized residuals. residuals. Density plots show clearly the left skewed shape of densities. Table 10 shows the p-values of Run tests for Italy and Spain. The p-values of two countries show that none of the three models failed for independence on the level of 0.01. Thus, for countries, all models except male LLRW of Spain and total LLRW of Spain, pass both diagnostic procedures and thus are valid ones. Table 11 presents the best model for each type of data based on AIC, BIC and SSPE for the valid models for Italy and Spain. For Italy, the LT model is the best for female and the LLRW models are the best for both male and total deaths, and for Spain, the LT models are the best for all three types. Table 12 shows estimates of unknown parameters of the best-fitted model for each type of deaths in Italy and Spain.

4.6. Asia-Pacific countries, Taiwan and Australia

Figure 11 shows plots of death counts of female, male and total deaths for Taiwan and Australia. Deaths of these two countries show a clear linear trend similar to deaths of the two North American countries mentioned above. Table 13 shows p-values of Shapiro tests for testing normality of the standardized residuals of models fitted to three types of death counts in Taiwan and Australia. The p-values show that none of models for all three types of deaths in Taiwan and Australia failed for normality on the level of 0.01. Figure 12 and 13 show the normal QQ plots for Taiwan and Australia, respectively. These plots confirm the results of Shapiro-Wilk tests. Table 14 shows the p-values of Run tests for Taiwan and Australia. The p-values of two countries show that none of the three models failed for independence on the level of 0.01. Thus, for two countries, none of the three models failed to be valid. Table 15 presents the best model for each type of data based on AIC, BIC and SSPE for the valid models for Taiwan and Australia. For Taiwan, the LT models are the best for both female and total deaths and the FT model is the best for male deaths, and for Australia, the FT models are the best for all three types of deaths. Table 16 shows estimates of unknown parameters of the best-fitted model for each type of deaths in Taiwan and Australia.

4.7. Signal extraction for unobserved components

One advantage to fit a structural time series model using the Kalman filter is that unobserved components of the time series can be estimated. Both the LLRW and the FT model have an unobserved trend component, μ_t. The LT model has two unobserved components, trend and slope, μ_t and β_t, respectively. The Kalman filter provides filtered estimates for unobserved components. For U.S. deaths, it is noted from Table 3 that the best-fitted model for female, male and total deaths is the FT, the LT and the FT model, respectively. Figure 14 shows filtered estimates of μ_t of U.S. female deaths, filtered estimates of μ_t and β_t of U.S. male deaths, and filtered estimates of μ_t of U.S. total deaths.

This paper analyzes female, male and total death counts of 8 countries in several regions in the world using three structural time series models with the Kalman filter. Three structural models are a local level with a random walk (LLRW) model, a fixed local linear trend (FT) model and a local linear trend (LT) model. LLRW model implies that the level of data moves stochastically based on a random walk model. Thus, this model is good to fit data without a clear linear pattern. FT model implies that the level of data moves stochastically with a deterministic slope. Thus, this model is good to fit data with a clear linear pattern. LT model implies that the level of data moves stochastically with a stochastic slope. That is, both level and slope move stochastically and thus this model is the most flexible. Thus, this model is good to fit data both with and without a clear linear pattern. Death counts of all three types of deaths in two North American countries and two Asia-Pacific countries show similar linear trends. Best fitted stochastic models for both North American countries and Asia-Pacific countries are either FT or LT model. Table 17 shows the signal to noise ratios for trend (VarN/VarE) and those of slope (VarK/VarE) of the best-fitted stochastic models for two North American countries. For the U.S., female deaths show higher signal to noise ratio for trend than that of male deaths. For Canada, however, female deaths show lower signal to noise ratio for trend than that of male deaths. Signal to noise ratios of slope are small for both countries. Note that NA in the table is for the FT model where VarK is not presented. Table 18 shows the signal to noise ratios for two Asia-Pacific countries. It shows that two countries do not have high signal to noise ratios for trend. Male deaths of both countries have higher signal to noise ratios than those of female death. Taiwan has low signal to noise ratios of slope as well. European countries do not show any clear linear trend. Best fitted stochastic models for European countries are either LLRW or LT. Table 19 and 20 show that European countries have low signal to noise ratios for both trend and slope.

Table 4, 8, 12 and 16 show the estimates of values of variances of irregular components, VarN, VarK and VarE. These variances show how much each component in the structural model moves stochastically up and down. For example, value of VarN shows how much the level of the trend, μ_t moves stochastically up and down over time, value of VarK shows how much the slope, β_t moves stochastically up and down over time and value of VarE shows how much the observation, y_t moves stochastically up and down around the unobserved trend, μ_t over time. Thus, the large variances do not imply that data has the large magnitudes of the trend and the slope, but imply that data has the large stochastic movements of the trend and the slope. Thus, large values of VarN and VarK of our data imply that our data has the trend and the slope which move a lot stochastically. To see the point, it is noted that for Australia, the best models for both female deaths and total deaths are the FT and from Table 18, ratios of VarN/VarE are similar each other, but from Table 16, VarN of total deaths is a lot larger than that of Female. Thus it implies that the level of the trend for total deaths moves stochastically more than that of female deaths. Figure 15 shows the estimates of the trend for both female deaths and total deaths from the Kalman filter. The plots of two estimated trends show that the trend of total deaths shows more volatility than that of female. Standard deviation of the trend of total deaths is 23064.75 and that of female deaths is 12681.7.

In Table 18 through 20, value of signal to noise ratio for trend (VarN/VarE) of LT model is zero for female model of Taiwan, male and total model of U.K., female and total model of Italy, and all models of Spain. For these cases, LT model (2.2) is adjusted as

This adjusted model implies that the level of the trend at time t moves stochastically based on only the level of trend and slope at time t − 1 without an irregular component of the trend. In this case, the stochastic movement of the level of the trend is purely based on the stochastic movement of slope, VarK.