Seasonal mortality forecasting based on influenca epidemics in the Netherlands

(1)

on Influenza Epidemics in The

Netherlands

Emanuela Festa

Master’s Thesis to obtain the degree in Actuarial Science and Mathematical Finance University of Amsterdam

Faculty of Economics and Business Amsterdam School of Economics Author: Emanuela Festa Student nr: 10746994

Email: emanuela festa@hotmail.com Date: October 29, 2016

Supervisor: Andrei Lalu

(2)

Statement of Originality

This document is written by Student Emanuela Festa who declares to take full responsibility for the contents of this document.

I declare that the text and the work presented in this document is original and that no sources other than those mentioned in the text and its references have been used in creating it.

The Faculty of Economics and Business is responsible solely for the supervision of completion of the work, not for the contents.

(3)

Abstract

It is important for pension funds and life insurance companies to be able to determine longevity and mortality risk with great accuracy, especially due to Solvency II regula-tions. Current mortality forecasting methods provide estimates of future mortality rates using observed annual data. Mortality, however, does not stay constant throughout the year. Studies have shown that economically developed countries such as The Nether-lands tend to have a much higher mortality rate in winter than in summer. Possible causes include influenza epidemics and cold winter temperatures. A higher degree of accuracy in mortality forecasting could perhaps be achieved if future mortality rates could be predicted using higher frequency recordings, for example, weekly observations instead of annual ones. To this end, this thesis focuses on weekly mortality rates in The Netherlands in the time span of 1995 to 2015. In order to explain their seasonal behaviour, correlations between mortality and temperature, and between mortality and influenza epidemics are studied. The standard Lee-Carter model is then used as the basis for designing seasonal mortality forecasting models. Different types of seasonal ARIMA models to estimate Lee-Carter’s time-varying parameter, the mortality index, are considered. One of the model specifications incorporates a flu indicator to estimate the excess winter mortality more accurately. The other models are seasonal ARIMA structures and ARIMA with seasonal dummy variables. The models are then evaluated according to their performance in out-of-sample forecasts and according to measures of residual variance. The results show that a flu indicator increases the out-of-sample forecast performance of seasonal ARIMA models significantly and residual variance is reduced the most when using an ARIMA with seasonal dummy variables. Finally, the thesis argues that a flu indicator model could be a useful tool for testing different future flu epidemics scenarios.

Keywords Lee-Carter, Seasonal ARIMA, Seasonal mortality, Excess winter mortality, Fore-casting mortality, Dutch mortality data, Influenza epidemics

(4)

Acknowledgements

First of all I would like to thank my supervisor Andrei Lalu for his availability and the discussions on seasonal mortality forecasting. My gratitude goes also to my husband Pedro Borges for his loving support during my studies and for contributing with the reviewing of my thesis. Finally, a big thank you to my family for their unconditional support, which made writing this thesis possible.

(7)

Introduction

Calculating longevity and mortality risk is important for pension funds and insur-ance companies providing pension plans and life insurinsur-ance. In Europe, in order to be compliant with Solvency II regulations, these institutions need to have enough capital reserved to be able to finance the liabilities they will incur over the next year with a probability of at least 99.5%. It is therefore very important to be able to estimate and forecast mortality as accurately as possible.

Existing mortality forecasting methods use observed annual mortality rates as input and hence yield annual forecasts. However, in certain countries, mortality does not stay constant throughout the year. In economically developed countries, mortality rates tend to increase significantly in the winter season, and in some years more than in others. A study byReichert et al.(2004) suggests that influenza epidemics play a very important role in causing this phenomenon. Another factor that could be responsible for higher mortality rates in winter are cold spells (Huynen et al.,2001).

In order to increase accuracy of mortality forecasts and therefore improve the de-termination of mortality risk, one could study the death rate observed on a monthly or weekly basis. Moreover, if a strong correlation between winter excess mortality and either average temperature or influenza epidemics can be found, one of these variables could be added as an indicator in mortality forecasting models based on higher frequency mortality rates. Creating a model that can forecast the mortality on a weekly basis in-stead of yearly might allow pension funds and insurance companies to manage reserves more effectively and estimate their risk for the coming years with more accuracy.

The forecasting method chosen for this research is the standard model by Lee and Carter(1992) due to its simplicity and empirical accuracy. The Lee-Carter model is one of the most frequently used models for mortality forecasting and various extensions, improvements and applications have been found since its introduction in 1992 (see for exampleLee,2000).

The aim of this thesis is to find appropriate models for forecasting the seasonal behaviour of the time-varying parameter of the Lee-Carter model, and hence of the central mortality rate. To this end, various specifications of seasonal autoregressive integrated moving average (ARIMA) models are considered. Furthermore, an improved ARIMA specification that uses information from a flu indicator to model excess winter mortality is proposed.

As input for the models weekly mortality rates in The Netherlands spanning 20 years (from 1995 to 2015) are used. The performance of these models is evaluated using sample forecasts and by comparing measures of residual variance. Finally, the out-of-sample forecasts of the seasonal models are compared to the out-of-out-of-sample forecast of a

(8)

2 Emanuela Festa — Seasonal Mortality Forecasting

yearly model by aggregating the weekly forecasts of the mortality rate into yearly ones. This is an interesting comparison to make because it could provide some evidence in support of the hypothesis that by aggregating higher frequency forecasts a higher degree of accuracy of yearly mortality forecasts can be obtained. If this is the case, modelling higher frequency mortality rates and then aggregating the forecasts could become a more accurate method for determining future yearly mortality, and hence mortality and longevity risk.

Here is a brief outline of how the thesis is structured. Chapter2provides a description of the types of data that were studied and their sources. Also, a discussion on the choice of variables to use for modelling excess winter mortality is given. Chapter 3 focuses on the theory of the models considered for forecasting yearly and seasonal mortality. The first section describes the standard Lee-Carter model and the estimation of its parameters, whereas the second section presents the non-age specific version of the model. The third section describes a method for forecasting the time-trend parameter of the Lee-Carter model, also called the mortality index. Sections 4 and 5 describe more complex methods that can be applied in order to model and forecast the seasonal behaviour of the weekly mortality index, and the last section describes how the research on flu epidemics can be applied in order to model excess winter mortality. Chapter 4 presents the results of the forecasts obtained with the models described in chapter 3 and evaluates their performance. Finally, in chapter5, the conclusions are presented.

(9)

Data description

Three types of data are considered for the research. The first type gives an indication of how many people are affected by influenza in The Netherlands each week and the second is a measure of the average temperature in The Netherlands per week. The third type concerns demographics: more specifically number of deaths and population size in The Netherlands per week. These statistics have been used in order to calculate the central mortality rate per week in The Netherlands over the last 20 years. The next three sections explain in detail what kind of data was obtained and from which sources. In the fourth section examples are shown of how the central mortality rate correlates with data on flu epidemics and with average temperature during winter seasons.

2.1 Influenza epidemics

The World Health Organization (WHO) talks in terms of “influenza seasons”. An influenza season is a time period that begins each year at week 40 and ends at week 20 in the following year. The advantage of talking in terms of influenza seasons instead of years is that normally an influenza season includes at most one flu epidemic whereas a year could possibly include two. This system is more useful because it makes it clear which flu epidemic one is referring to.

The Dutch government established that a flu epidemic is present in The Netherlands when, for two consecutive weeks, at least 51 out of 100.000 people per week have reported an influenza-like illness (ILI) to their general practicioner. One of the organisations that gathers this type of statistics in The Netherlands is the Netherlands Institute for Re-search of Healthcare (NIVEL). The NIVEL brings out weekly bulletins called “NIVEL Zorgregistraties eerste lijn” (Hooiveld et al.,2016) during the influenza seasons, report-ing, among other things, the number of ILI patients that have been registered by a certain number of general practitioners, referred to as the sentinel doctors, who are working in sentinel station practices. The sentinel stations are spread all over the coun-try but the total number and their distribution across The Netherlands may vary per year. This kind of registration is called Continuous Morbidity Registration (CMR). The NIVEL also publishes extensive yearly reports called “Continue Morbiditeits Registratie Peilstations Nederland”, which gather and analyse all the information collected in the previous year and compare it to information collected in earlier years. See for example

Donker(2012) andBartelds(2003). The reports are not limited to ILI but cover a wide range of pathologies.

The data from the weekly NIVEL bulletins and from the yearly NIVEL reports published in the last 20 years has been used in this research to derive an indicator for the intensity and length of flu epidemics from 1995 to 2015.

(10)

2.2 Average temperature

Daily and monthly average temperatures are available from the European Climate Assessment & Dataset (ECA&D) (Klein Tank,2002). In order to obtain weekly average temperatures the daily average temperatures have been averaged to form weekly ones. The location in The Netherlands where the temperatures were measured is De Bilt. This location is chosen because it is in the centre of The Netherlands and therefore gives more or less an average of the temperatures in the whole country. In the North and more inland towards the South and the East, temperatures tend to be slightly lower during winter than in the centre of The Netherlands, whereas in the West they tend to be slightly higher. Overall, however, average temperatures in different regions of The Netherlands are fairly similar.

2.3 Central mortality rate

The age specific central mortality rate mx,t is defined as the average death rate

experienced by a population group of age x at time t. It is calculated as follows:

mx,t =

dx,t

ex,t

, (2.1)

where dx,t represents the number of people aged x that have died at time t and ex,t,

also known as the exposure to mortality risk, represents the average number of people aged x living at time t.

In general time t is expressed in terms of years and it is unusual for actuaries to work with non-annual demographic data. Death rates of a higher frequency (for example monthly or weekly) are not as common and therefore harder to obtain for a specific time-span. The time-span that can be studied is therefore determined by the availability of weekly data from Statistics Netherlands (CBS). This is a time-span ranging from 1995 up to and including 2014. For this time-span the CBS has provided the number of deaths in The Netherlands per week and the number of people living in The Netherlands per month (measured at the beginning of the month). Weekly population figures were estimated from the monthly figures using a method described in Appendix A. The data on number of deaths and population is also available for female and male population subgroups and (for shorter time-spans) for different age groups.

Due to data availability constraints age differentiation will be excluded in the model. Furthermore, instead of the yearly behaviour of the mortality time trend the seasonal behaviour will be studied. To this end the non-age specific central mortality rate, de-noted by mtwill be used. Here t represents the week number belonging to a particular

year. For simplification it is assumed that a year consists of exactly 52 weeks. The non-age specific central mortality rate is defined as follows:

mt=

dt

et

, (2.2)

where dtrepresents the number of deaths in week t and et, represents the average number

(11)

2.4 Correlation of m

t

with flu data and average

tempera-tures in winter

The first part of this section evaluates the correlation between average temperature and the central mortality rate mt, whereas the second part addresses the correlation

between flu epidemics and mt. Based on these evaluations, the most suitable variable

for fitting seasonal mortality is chosen. Only data from the winter seasons (starting on week 40 in year y and ending on week 20 in year y + 1) is shown because measurements regarding influenza are only published for this period of time by organisations such as the NIVEL.

2.4.1 Correlation between mt and average winter temperatures

In order to study the correlation between mtand weekly average temperature, the R2

values of these variables were calculated for all winter seasons occurring between 1995 and 2015. The plots in figure 2.1 illustrate four examples (taken in five year intervals) of the correlations found. The regression lines and R2 values belonging to the plots are given in table 2.1.

Figure 2.1: Correlation between mt and average temperature for four specific winter

seasons.

Winter Season Linear Regression Equation R2 1997/1998 y = 9.2 × 104x − 8 6.2% 2002/2003 y = −1.8 × 105x + 37 17.9% 2007/2008 y = −2.0 × 105x + 41 19.5% 2012/2013 y = −2.2 × 105_{x + 43} _41.1%

Table 2.1: Linear regression equations and R2values from the scatter plots in figure2.1. Three of the plots in figure2.1suggest a negative correlation between average temper-ature and mt, whereas in one of the plots (winter 1997/1998) it appears to be positive.

Out of the 20 plots made, however, only one finds a positive correlation. Therefore, the correlation between average temperature and mt is predominantly negative, as

(12)

The histogram in figure2.2illustrates the distribution of the 20 R2 values calculated. They range from 6% to 61% and the median was found to be 31%.

R2 Frequency 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 0 1 2 3 4 5 6

Figure 2.2: Histogram showing the R2values between average temperature and mortality for all the winter seasons between 1995 and 2015.

2.4.2 Correlation between mt and number of recorded ILI patients

The correlation between mt and the number of ILI patients recorded by sentinel

doctors was calculated for all winter seasons occurring between 1995 and 2015. In figure 2.3four examples of correlation plots between these two variables are shown for the same winter seasons considered in figure2.1. The regression lines and R2 _{values belonging to}

the plots are given in table2.2.

Figure 2.3: Correlation between mtand number of recorded ILI patients for four specific

winter seasons.

The correlations in the plots of figure 2.3 were found to be positive. Positive cor-relations were also found for the remaining 16 years in the time range from 1995 to 2015.

(13)

Winter Season Linear Regression Equation R2 1997/1998 y = 3.0 × 106_{x − 473} _59.0%

2002/2003 y = 9.2 × 105x − 134 29.4% 2007/2008 y = 5.6 × 105x − 53 11.3% 2012/2013 y = 3.0 × 106x − 383 75.3%

Table 2.2: Linear regression equations and R2values from the scatter plots in figure2.3. The histogram in figure2.4illustrates the distribution of the 20 R2 _{values calculated.}

The values range from 7% to 87% and the median was found to be 57%.

R2 Frequency 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 0 1 2 3 4 5 6

Figure 2.4: Histogram showing the R2 values between recorded ILI patients and mor-tality for all the winter seasons between 1995 and 2015.

The histograms show that the correlation strength between flu epidemics and mt

varies considerably from year to year but is overall stronger than the correlation found between average temperature and mortality. In fact, the correlation between flu epi-demics and mt was found to be stronger for 16 out of the 20 years considered, that is

to say in 80% of the cases. Furthermore, the median R2 _{value was found to be 57% for}

flu epidemics versus mt, and only 31% for average temperature versus mt. Following

these observations, the decision was made to use the recorded ILI patients as exogenous variable for modelling seasonal mortality in The Netherlands.

(14)

Chapter 3

Modelling of weekly Dutch

population mortality

The modelling of weekly Dutch population mortality from 1995 to 2015 in this re-search is based on the standard Lee-Carter model (Lee and Carter,1992). As mentioned in the introduction, this model is chosen due to its simplicity and empirical accuracy.

Section 3.1 briefly explains the theory of the model presented in Lee and Carter’s paper from 1992, Modeling and Forecasting US Mortality. Section 3.2 describes how the non-age specific version of the Lee-Carter model is applied in order to conduct research on seasonal mortality in The Netherlands. The remaining sections describe various ARIMA models that were used to fit the κtparameter of the Lee-Carter model,

also known as the mortality index.

In section3.3the base specification for the ARIMA model used to fit κtis described.

It is a simple ARIMA(0,1,0) with constant, which can also be referred to as a random walk with drift (RWD). However, the mortality index determined per week over a 20 year time span has a strong seasonal behaviour due to excess mortality occurring in the winter season. This seasonal behaviour cannot be modelled using a simple RWD model and therefore more sophisticated models need to be found.

In order to capture seasonality, other types of ARIMA models are explored. Section 3.4describes seasonal ARIMA (SARIMA) structures where both the seasonal and non-seasonal behaviours of κtcan be modelled. Section3.5describes an ARIMA(0,1,0) model

with 51 seasonal dummy variables and section3.6introduces a SARIMA model in which an influenza indicator based on the number of reported flu cases in the winter season is used as an exogenous regressor to fit seasonal behaviour of the mortality index.

3.1 The Lee-Carter Model

Lee and Carter published in 1992 a very easy to use model for fitting and forecast-ing the central mortality rate. Accordforecast-ing to this model, the age specific log-mortality ln(mx,t) is given by

ln(mx,t) = αx+ βxκt+ σxt. (3.1)

The αx and βx parameters are age specific, whereas the κt parameter, a stochastic

process also known as the mortality index, represents the change in mortality rates over time. The αx parameter can be interpreted as the average mortality belonging to a

certain age group and the βx parameter as the age sensitivity of the mortality index.

(15)

The error term t is assumed to be independent and identically distributed with mean

0 and variance 1.

Without restrictions the model is undetermined. One can easily verify that for any constant c, the following set of parameters will lead to the exact same model: {αx, βx, κt},

{αx, βx/c, c · κt} or {αx − c · βx, βx, κt+ c}. To make the model identifiable, Lee and

Carter imposed the following constraints on βx and κt:

X x βx= 1; (3.2) X t κt= 0. (3.3)

From these restrictions it follows that the αx parameter is equal to the age specific

log-mortality averaged over time:

αx= 1 T T X t=1 ln(mx,t). (3.4)

3.1.1 Estimation of the age specific and time-varying parameters Equation 3.4 can be used to estimate the αx parameter:

ˆ αx= 1 T T X t=1 ln(mx,t). (3.5)

The estimation of βx and κtcan be obtained by finding the least squares solution to

X

x,t

(Zx,t− βxκt)2, (3.6)

where Zx,t = ln(mx,t) − ˆαx (Koissi and Shapiro, 2008). This problem was solved by

Lee and Carter through Singular Value Decomposition (SVD). According to SVD, ma-trix Zx,t can be decomposed into a series of singular column-vectors Ux,i multiplied by

singular row-vectors VT

t,i, multiplied by factors

√

ηi, where i = 1, ..., r and r = rank[Zx,t]:

Zx,t = r X i=1 √ ηiUx,iVt,iT. (3.7)

Zx,t can then be approximated using just the first component, which corresponds to

the highest eigenvalue:

Zx,t≈

√

(16)

The estimation of βxand κthence follows from the minimisation of equation3.6and

from equation 3.8: ˆ βx = Ux,1 P iUx,i ; (3.9) ˆ κt= √ η1Vt,1T X i Ux,i. (3.10)

Note that in expressing ˆβx the first constraint imposed byLee and Carter, given by

equation 3.2, has been taken into account, and hence it is normalised.

The mortality time trend in this model is the same for all ages. This implies that the forecasting of future central mortality rates for all age groups simplifies to forecasting future values of κt. The modelling and forecasting of the mortality index κtis discussed

from section3.3 onwards.

3.2 Application of the Lee-Carter Model to weekly

non-age specific mortality

The non-age specific version of the Lee-Carter model can be seen as a special case of the model described in the previous section. In this case, the log-mortality is given by

ln(mt) = α + κt+ σt, (3.11)

where the α parameter is equal to the average log-mortality over time and the β param-eter is equal to 1 because this model does not include the age sensitivity of the mortality time trend κt. Both α and β can be expressed as one single term instead of as vectors.

The error term t is again assumed to be independent and identically distributed with

mean 0 and variance 1. This error term represents the variability that is not explained by the model, and which this thesis aims to decrease by looking into seasonal patterns. Lee and Carter used yearly measurements of the central mortality rate in order to de-termine κt but in order to model seasonality weekly measurements are used. Therefore,

t denotes the week number belonging to a particular year instead of a year.

Note that the restrictions given in equations 3.2 and 3.3 still apply to this model. These constraints leads to

α = 1 T T X t=1 ln(mt). (3.12)

3.2.1 Estimation of the parameters

For the estimation of α the result obtained in equation 3.12can be used. It is equal to the average log-mortality over time:

(17)

ˆ α = 1 T T X t=1 ln(mt). (3.13)

The β parameter is equal to 1 in the one dimensional case (no age subgroups) and therefore no estimation is needed.

κt can simply be estimated using equation 3.11:

ˆ

κt= ln(mt) − ˆα. (3.14)

Note that the one dimensional case has been formulated using the same notation as in section 3.1.1in order to be consistent with the general specification.

3.3 Non-seasonal ARIMA models

When κtis calculated using the central mortality rate in The Netherlands per week

a strong seasonal variation is found. A time series that is measured at this frequency can be decomposed in a trend component Tt, a seasonal component St and an error

component Rt. An illustration of the seasonal behaviour of κt is shown in figure 3.1.

The figure shows values for all weeks of the year from 1995 up to and including 2014. The trend component, which was obtained using a moving average (method: symmetric window of length 52 with equal weights) is shown in red.

Year κt 1995 2000 2005 2010 2015 −0.1 0.0 0.1 0.2 0.3 κt Tt

Figure 3.1: Plot of κtper week obtained from Dutch mortality rates from 1995 to 2015.

The red line denotes the trend Tt, obtained using a moving average.

In order to model the stochastic process κt various types of ARIMA(p,d,q) models

can be used. The non-seasonal type is described in this section and more complex models that include seasonality are described from section 3.4onwards.

(18)

3.3.1 Model specification

ARIMA models are named this way because they may include autoregression (AR), differencing and moving averages (MA). In an ARIMA(p,d,q) time series, p represents the order of the autoregressive part, d is the order of integration and q represents the order of the moving average part (Heij et al., 2004). The general formulation of this type of time series is

φ(L)(1 − L)dyt= φ0+ θ(L)t, (3.15)

where

• L is the lag operator which shifts the data back by one or more periods. For example: Lyt= yt−1 and Ldyt= Ld−1yt−1= yt−d;

• (1 − L)d _{is an expression for the differencing process, which takes place in order}

to remove trends or seasonality and make the time series stationary. A differenced time series is said to have integration order d;

• φ(L), which can also be written as AR(p), is the autoregressive equation with parameters φi where i = 1, ..., p:

φ(L) = 1 − φ1L − φ2L2− ... − φp−1Lp−1− φpLp; (3.16)

• θ(L), which can also be written as MA(q), is the moving average equation with parameters θj where j = 1, ..., q:

θ(L) = 1 + θ1L + θ2L2+ ... + θq−1Lq−1+ θqLq; (3.17)

• φ0 is a constant that is equal to the drift of the time series when d > 0;

• t is a normally independent, identical distributed error term with mean 0 and

variance σ2, also known as a white noise process.

3.3.2 Random Walk with Drift model

Lee and Carter (1992) used a random walk with drift model (RWD), which is an

ARIMA(0,1,0) with constant, in order to fit and then forecast the mortality index κt.

The process can be written as

κt= φ0+ κt−1+ t, (3.18)

where φ0 is the drift of the time trend and t the error term with mean 0 and variance

σ2 (Vellekoop,2016). Note that the above equation can also be written as

κt= κ0+ tφ0+ t

X

s=1

s. (3.19)

(19)

ˆ φ0=

κt− κ0

t , for t > 0. (3.20)

Here the error term vanishes since it has zero mean. Furthermore, the residual variance can be estimated with the following equation:

ˆ σ2= 1 t − 1 t X j=1 (κj − κj−1− ˆφ0)2, for t > 1. (3.21)

This means that the drift and residual variance parameters can be directly obtained from the data.

For the yearly data used byLee and Carter(1992), this model provided a fairly good fit and gave reasonable forecasts. However, when using weekly data the drift should not be determined with equation 3.20 because the first and last values of κt can vary

significantly depending on which year (and on which week of the year) the time series ends and begins. In this case, these values represent the mortality index in week 1 (of year 1995) and week 52 (of year 2014) but the mortality index in any specific week of the year can vary considerably from year to year. The difference between the first and last values of the time series can be equal to any in a wide range of values, both positive and negative, and the resulting drift would therefore, very probably, not represent the actual, observed drift of the time series. Another method for determining the drift of the weekly κt time series should therefore be used.

The drift can alternatively be determined by fitting a regression line to the time trend Ttof κt, which is shown in red in figure3.1. The regression line of Ttis illustrated

in figure3.2, and through this method the drift is found to be ˆφ0 = −1.00 × 10−4. This

is the value that will be used for the drift in all models discussed in this thesis for fitting the mortality index.

Figure 3.2: Trend Tt of the mortality index per week (red line) and linear regression

(black line) for determining the drift of the time series.

3.3.3 Mortality forecasting with RWD

In order to forecast mortality using the Lee-Carter model it suffices to forecast future values of κt. When κt follows a simple RWD model and all the values of κt are known

(20)

up to and including time T , the forecast κT +f can be written as

κT +f = κT −1+f + φ0+ T +f. (3.22) This is equivalent to κT +f = κT + f φ0+ f X h=1 T +h. (3.23)

Since the error terms have zero means, the expected value of the forecast for κT +f will

be simply equal to κT+f φ0. Furthermore, κT +f is normally distributed and has variance

σ2

f . Hence the 95% confidence interval will be equal to ±1.96 standard deviations, or

±1.96σ

√ f .

This result can then be used to find the expected value of the forecasted central mortality rate:

E[(mT +f)] = exp{ ˆα + ˆκT + f ˆφ0}. (3.24)

3.4 Seasonal ARIMA structures

One way to take seasonality into account in ARIMA models is to apply seasonal differencing. A seasonal ARIMA model can be written as SARIMA(p,d,q)(P ,D,Q)m,

where m = number of periods in a year, and can include both a non-seasonal part, denoted by (p,d,q), and a seasonal part denoted by (P ,D,Q)m. The general equation of

a SARIMA with constant is of the following form:

φ(L)Φ(L)(1 − L)d(1 − Lm)Dyt= φ0+ θ(L)Θ(L)t, (3.25)

where, in addition to the non-seasonal terms described in section3.3.1

• (1 − Lm₎D _{is an expression for the seasonal differencing process of order D;}

• Φ(L) is the seasonal autoregressive equation AR(P ) with parameters Φi, where

i = 1, ..., P :

Φ(L) = 1 − Φ1Lm− Φ2L2m− ... − ΦP −1L(P −1)m− ΦPLP m; (3.26)

• Θ(L) is the seasonal moving average equation MA(Q) with parameters Θ_j, where j = 1, ..., Q:

Θ(L) = 1 + Θ1Lm+ Θ2L2m+ ... + ΘQ−1L(Q−1)m+ ΘQLQm. (3.27)

Note that in the seasonal terms, all lags are elevated to the power m. For a more in depth description of this type of seasonal ARIMA structures the reader may refer to

(21)

3.4.1 SARIMA(0,1,0)(0,1,0)52 with drift model

As equation 3.25 shows, SARIMA specifications can become complex. To illustrate how κt can be fitted using a SARIMA model, a relatively simple model will be taken

as an example: the SARIMA(0,1,0)(0,1,0)52 with drift. This model takes into account

a general yearly trend in the non-seasonal part and applies differencing in the seasonal part (with m = 52 for weekly data) in order to model the seasonal trend as well.

Equation 3.25 reduces to:

(1 − L)(1 − L52)κt= φ0+ t, (3.28)

which is equivalent to

κt= κt−1+ κt−52− κt−53+ φ0+ t. (3.29)

The forecasted values of κt in this model are obtained with:

κT +f = κT −1+f + κT −52+f− κT −53+f+ φ0+ T +f. (3.30)

Another way of expressing this is:

κT +f − κT −52+f = κT − κT −52+ f φ0+ f

X

h=1

T +h. (3.31)

Therefore, the expected values of the forecasts can be found through:

E[(κT +f)] = κT −52+f+ κT − κT −52+ f φ0. (3.32)

Note that when forecasting κT +f values where f > 52, the previous forecasts (where

f < 53) need to be computed first.

3.5 ARIMA model with weekly dummy variables

Another way to model the seasonal behaviour shown in figure 3.1, is to assume that every single week in a year follows its own specific time trend. In this case an ARIMA time series can be expressed in the most general form as

φ(L)(1 − L)dyt= βt+ θ(L)t, (3.33)

where the drift βttakes on 52 different values, depending on which week of the year the

(22)

This model can also be described using so-called dummy variables. For a detailed explanation on the use of dummy variables for modelling parameter variation the reader may refer toHeij et al.(2004, p. 303-319). Dummy variables Dh (where h = 1, 2, ..., 52)

are defined as follows: Dht is equal to 1 when the tth observation falls in week h and is

equal to zero otherwise. Equation 3.33can subsequently be re-written as

φ(L)(1 − L)dyt= β1D1t+ β2D2t+ ... + β52D52t+ θ(L)t. (3.34)

However, since it is preferable to have a constant term in the model (denoting the drift), the first dummy variable D1 is eliminated and the first term on the right hand

side of equation3.34 is replaced by φ0:

φ(L)(1 − L)dyt= φ0+ 52

X

i=2

βiDit+ θ(L)t, (3.35)

where φ0 describes the drift that the time trend would have if only the first week of

every year was taken into account. Then β2 is equal to the difference one would need to

add or subtract to φ0 if only week 2 was taken into account every year instead of week

1. Similarly for β3 up until β52.

To fit the mortality index κt a simple ARIMA(0,1,0) is used in combination with

dummy variables and equation 3.35reduces to

κt= κt−1+ φ0+ 52

X

i=2

βiDit+ t. (3.36)

This is an interesting model to investigate because it models each week of the year separately. One negative aspect of this model is that there are a lot of parameters to estimate. The second is that the mortality rate can vary significantly in the same week from one year to the next and therefore the estimated parameters can be inaccurate.

The forecasts of κt are given by:

κT +f = κT +f −1+ φ0+ 52 X i=2 βiDi(T +f )+ T +f, (3.37) which simplifies to κT +f = κT + f φ0+ f 52 X i=2 βiDi(T +f )+ f X h=1 T +h, (3.38)

so that the expected values of the forecasts can be found through

E[(κT +f)] = κT + f φ0+ f 52

X

i=2

(23)

3.6 SARIMA with flu indicator as exogenous regressor

This section describes how to use an indicator representing the intensity of an in-fluenza epidemic to act as an exogenous regressor for the mortality index κt. The purpose

of the indicator is to take into account the excess mortality in The Netherlands during the winter season and hence improve the quality of the fit and of the mortality forecast. The main advantage of using a predictor in the model such as a flu indicator is that, when information on imminent influenza epidemics is available, the mortality forecasts for the immediate future can be improved.

In general the model can be expressed as

φ(L)Φ(L)(1 − L)d(1 − Lm)Dyt= cIt+ φ0+ θ(L)Θ(L)t, (3.40)

where Itis the influenza indicator and c is a parameter to be estimated.

For a SARIMA(0,1,0)(0,1,0)52with drift, with flu indicator, for example, the model

is given by

κt= κt−1+ κt−52− κt−53+ cIt+ φ0+ t. (3.41)

The forecasts for this model are

κT +f = κT −1+f + κT −52+f− κT −53+f + cIT +f + φ0+ T +f. (3.42) This simplifies to κT +f− κT −52+f = κT − κT −52+ c f X h=1 IT +h+ f φ0+ f X h=1 T +h, (3.43)

and the expected value of the forecasts are given by:

E[(κT +f)] = κT −52+f + κT − κT −52+ c f

X

h=1

IT +h+ f φ0. (3.44)

3.6.1 Constructing a flu indicator

One way to construct the flu indicator Itin the model is to simply setting it equal to

the number of flu cases registered by the NIVEL. In figures3.3and3.4the flu indicator was plotted together with a multiple of κt for comparison. Figure 3.3 compares these

two variables from 1995 to 2005 and figure 3.4 shows how they compare from 2005 to 2015. The figures show that during the winter season the two are very highly correlated. The peaks of the mortality index κtcoincide quite well with the peaks of the number

of ILI cases registered by the NIVEL, with the exception of the 2009/2010 winter season, when a different kind of influenza epidemic (the A(H1N1), also known as “swine flu”)

(24)

broke out. According to the NIVEL, the A(H1N1) influenza affects mostly the younger part of the population, causing less deaths than the A(H2N2) and A(H3N2) viruses, which affect mostly the elderly. In order to improve the quality of the indicator one could choose to exclude registered influenza cases affected by the A(H1N1) virus. See also the paper by Reichert et al. (2004) whose study finds that excess mortality in winter had a far higher correlation with the occurrence of A(H2N2) and A(H3N2) virus epidemics than with A(H1N1) and B virus epidemics.

Occasionally, the mortality index peaks slightly precede the peaks of the number of ILI cases (see, for example, the winters of 1997, 1999 and 2009). This seems illogical but could perhaps be caused by delays in registration of the ILI cases.

Figure 3.3: Number of flu cases per week registered by NIVEL and κt× 1000 from 1995

to 2005.

Figure 3.4: Number of flu cases per week registered by NIVEL and κt× 1000 from 2005

to 2015.

The graphs also show that the summer seasons during time period 2005 to 2015 are not well described by the indicator. The troughs of κt, unlike the peaks, differ

considerably from the “troughs” of the flu indicator, which are simply equal to zero. Sharp peaks occurring in summer can also be observed (albeit smaller than the winter

(25)

peaks). These could possibly be caused by heat waves but a correlation study between temperature and mortality in the summer season would have to be conducted in order to answer this question with more certainty. The different summer behaviour will decrease the quality of the fit and one could consider modifying the indicator such that negative values are inserted to indicate the summer seasons, or even include an extra indicator for the summer. However, modelling summer troughs goes beyond the scope of the thesis and it will therefore be limited to forecasting winter peaks.

(26)

Chapter 4

Results

This chapter presents the results obtained from forecasting mortality using the four ARIMA models described in chapter 3. Section 4.1 shows three year forecasts of the mortality index κtper week and compares the relative performance of the models.

Sec-tion4.2compares how well the models perform in out-of-sample forecasts, by measuring the difference between forecasted and observed values of κt and mt. The section

con-cludes with a comparison between the forecasting accuracies of the seasonal models and that of a yearly RWD model. And finally, section 4.3 shows the results obtained for long term forecasts (twenty years projection). The parameters of all the ARIMA mod-els mentioned in this chapter can be found in Appendix B and the R code used for the calculations can be found in Appendix C.

4.1 Full sample forecasts

Three year forecasts of κthave been made with the following models: Random Walk

with Drift, SARIMA(p,d,q)(P,D,Q)52 with drift, ARIMA(0,1,0) with weekly dummy

variables and SARIMA(p,d,q)(P,D,Q)52 with drift, with flu indicator as exogenous

re-gressor. The models were estimated with observations from 1995 to 2015. At the end of the section a table summarises the main properties of the forecasts.

In the forecast graphs presented in this chapter, the black solid line represents ob-served values of κt whereas the solid blue line represents forecasted values. Two

confi-dence intervals are shown in the graphs of the first section: the smaller one with darker shade shows the 80% confidence interval and the larger one, shaded light grey, shows the 95% confidence interval.

4.1.1 Random Walk with Drift

When the model used for fitting κtis an ARIMA(0,1,0) with constant, only the trend

can be forecasted since the seasonal component is not modelled. Figure 4.1 shows the three year forecast of the trend of κt. The residual variance ˆσ2 is 1.97 × 10−3 and the

confidence intervals are quite large, indicating the imprecision of this model.

The drift ˆφ0 was calculated in chapter3 as −1.00 × 10−4. The calculation was done

using the trend line in figure3.2rather than using the Arima function from the forecast package (Hyndman,2015) in R, in order to avoid obtaining an incorrect drift value. The Arima function calculates the drift using equation 3.20, thereby only considering the first and last values of the time series (κ1 and κ1040). As explained in chapter 3, the

difference between the first and last values of the time series can be equal to any in a 20

(27)

Forecasts from ARIMA(0,1,0) with drift Year κt 1995 2000 2005 2010 2015 −1.0 −0.5 0.0 0.5 1.0 κt forecast

Figure 4.1: Three year forecast of the mortality index κtper week using the RWD model.

wide range of values, both positive and negative, and the resulting drift would therefore be quite different from the general drift of the time series.

In order to obtain the drift estimated in chapter 3 a correction was applied to the first and final values of the observed time series. κ1 was assigned a value equal to the

first point of the regression line and κ1040 to the last point (times a factor to account for

the different number of data points). Following this correction, when the Arima function was used to fit the time series a drift ˆφ0 of −1.00 × 10−4 was obtained, with a standard

error (s.e.) of 1.40 × 10−3, implying the drift estimation is not significant.

Since the forecast provides no information whatsoever on the future seasonal be-haviour of κt, the model is clearly more suitable for fitting yearly rather than weekly

time series. Therefore, the rest of this chapter will focus on the other three types of ARIMA models.

4.1.2 SARIMA(p,d,q)(P,D,Q)52 with drift

If seasonal effects are taken into account by fitting κtusing a SARIMA(p,d,q)(P,D,Q)52

with drift structure, as described in section3.4, a seasonal trend appears in the forecast. Different SARIMA(p,d,q)(P,D,Q)52 models were tested using the standard

Box-Jenkins methodology for ARIMA model selection (Box et al., 2015). Amongst these the SARIMA(2,0,1)(0,1,1)52 with drift was found to be the best fit for κt. Figure 4.2

shows the three year forecast of κtusing this model. The forecast has a residual variance

ˆ

σ2 of 1.68 × 10−3. This is a large improvement with respect to the RWD model. The confidence intervals are also much narrower than the ones shown in figure 4.1.

When the model is chosen using the Box-Jenkins methodology, a correction for the drift, like the one that was applied in the RWD model, is no longer necessary since higher orders of p, q and Q render the model more accurate. The drift estimated with the Arima function in R is equal to the drift that was estimated in chapter3using figure 3.2(φ0= −1.00 × 10−4). See Appendix B for the estimation of the remaining parameter

(28)

Forecasts from SARIMA(2,0,1)(0,1,1)[52] with drift

Year κt 1995 2000 2005 2010 2015 −0.2 −0.1 0.0 0.1 0.2 0.3 κt forecast

Figure 4.2: Three year forecast of the mortality index κt per week using the

SARIMA(2,0,1)(0,1,1)52 with drift model.

4.1.3 RWD with weekly dummy variables

Forecasts from RWD with Dummy Variables

Year κt 1995 2000 2005 2010 2015 −1.0 −0.5 0.0 0.5 1.0 κt forecast

Figure 4.3: Three year forecast of the mortality index κt per week using a RWD with

seasonal dummy variables.

Figure 4.3 shows the three year forecast obtained by using a RWD with dummy variables for each week of the year (see section 3.5). In Appendix B the values of the estimated parameters (including 51 dummy variables) can be found.

Similarly to the forecast obtained using the SARIMA(2,0,1)(0,1,1)52with drift model,

each period (of length 1 year) looks exactly the same. The shape of a single period in the forecast is the result of fitting each of the 52 weeks in a year with separate parameters.

(29)

The residual variance ˆσ2 is lower than the one in the RWD model but higher than in the SARIMA model: 1.75 × 10−3.

In a RWD model the drift is determined by the first and last values of the time series so also in this case a correction needs to be made in order to ensure that the forecast uses the correct drift ( ˆφ0 = −1.00 × 10−4). The correction was made by first setting κ1

and κ1040 (the first and last values) equal to the average of all the first and last values

of the years in the sample. Then it was calculated that in order to obtain the correct drift the difference between the first and last value needs to be -0.104. Half of this value was therefore subtracted from κ1 and added half of it to κ1040.

Similarly to the SARIMA models, the fit and drift with this type of model are also improved when higher orders of p and q are considered. Appendix B also includes the parameters estimated for an ARIMA(5,1,1) with dummy variables, which is the best model with dummy variables according to Box-Jenkins methodology (Box et al.,2015). The coefficients of the estimated ˆβ parameters (the so-called dummy variables) are ap-proximately equal to the coefficients of the ˆβ parameters belonging to the ARIMA(0,1,0) with dummy variables model. The standard errors, however, are reduced considerably (up to a factor of 2). Also ˆσ2 is considerably reduced.

What all the models so far fail to predict are exceptional spikes in the winter season, indicating higher than usual excess mortality. The peaks (and troughs) in the forecast are all very regular while the time series clearly shows that there is great variation in the height of the peaks. This aspect of the time series is considered in the model described in section 3.6, which takes into account the intensity of influenza epidemics during the winter season. The results of the forecasts obtained with the flu indicator model are shown in the following subsection.

4.1.4 SARIMA(p,d,q)(P,D,Q)52 with drift, with flu indicator

Forecasts from SARIMA(2,0,1)(0,1,1)[52] with drift, with Flu Indicator

Year κt 1995 2000 2005 2010 2015 −0.2 −0.1 0.0 0.1 0.2 0.3 κt forecast

Figure 4.4: Three year forecast of the mortality index κt per week using a

SARIMA(2,0,1)(0,1,1)52 with a flu indicator as exogenous regressor.

If the intensity of influenza epidemics in the coming winter season is known approx-imately, the prediction of excess mortality in winter could be made more accurately

(30)

using the model described in section3.6. This method is therefore mainly applicable for short-term mortality forecasting.

Figure 4.4shows a hypothetical three year forecast of the mortality index in winter. The flu indicator predicting influenza intensity for years 2015 to 2018 was chosen to be equal to the one used for years 1995 to 1998. From the graph it can be seen that the three predicted peaks follow the same pattern of relative heights as the peaks in years 1995 to 1998, however they don’t reach the same level of mortality because of the declining mortality rate observed in the second decade relative to the first one.

Similarly to the SARIMA(2,0,1)(0,1,1)52 model without flu indicator, there was no

need to correct the drift since the model, selected through Box-Jenkins methodology, manages to capture the correct drift by including higher orders of p, q and Q. Fur-thermore, the residual variance ˆσ2 is a lot lower than in the other models presented: 1.54 × 10−3.

The ˆc parameter in equation3.41was estimated as 6 × 10−4 (with a s.e. of 1 × 10−4) and the coefficients of the remaining parameters can be found in Appendix B. It’s interesting to note that all the parameter coefficients, apart from the one belonging to Θ1 (which is already highly significant), improve their significance considerably with

respect to the parameter coefficients of the SARIMA(2,0,1)(0,1,1)52 model without flu

indicator.

4.1.5 Model evaluation

This subsection shows how the models discussed in this section perform in terms of the “Akaike Information Criterion” (AIC), the “Bayesian Information Criterion” (BIC) and the residual variance ˆσ2. These criteria are useful in selecting the best ARIMA models because they all represent a measure for the variance (or information loss). Therefore, the lower their values, the better the model. However, model selection should also take into account the relative complexity of the models. When considering models that yield a similar goodness of fit, less complex models should be preferred to more complex ones. To this end, the AIC and BIC also include a penalty term which is proportional to the number of parameters fitted.

The AIC and BIC can be calculated as follows:

AIC = −2 ln(L) + 2k, (4.1)

BIC = −2 ln(L) + k ln(n), (4.2)

where

• L is the maximum value of the likelihood function for the model;

• k denotes the total number of parameters used (including the variance parameter ˆ

σ2);

• n denotes the number of observations used for estimation. Table4.1shows that the residual variance ˆσ2

is smallest when using an ARIMA(5,1,1)

model with weekly dummy variables. Despite the large number of parameters involved, this model gives a better performance than the ARIMA(0,1,0) model with weekly dummy variables also in terms of the AIC and BIC.

(31)

Model AIC BIC σˆ2

RWD -3516 -3507 0.00197

SARIMA(2,0,1)(0,1,1)52 with drift -3267 -3238 0.00168

RWD with weekly dummy variables -3536 -3274 0.00175 ARIMA(5,1,1) with weekly dummy variables, with drift -3746 -3454 0.00141 SARIMA(2,0,1)(0,1,1)52 with drift, with flu indicator -3358 -3324 0.00154

Table 4.1: An evaluation of the models used for forecasting the mortality index.

The RWD model is the model that performs least well in terms of residual variance. As noted earlier, the model can be used for forecasting the trend (although a correction for determining the right drift is needed) of κt but gives no information on its seasonal

behaviour.

The SARIMA(2,0,1)(0,1,1)52with flu indicator is the second best model when

consid-ering the residual variance. It has the advantage that it has less parameters to estimate than the ARIMA(5,1,1) model with weekly dummy variables, and produces a more accurate forecast of the winter mortality when information on flu epidemics is available.

4.2 Out-of-sample forecasts of κ

t

and m

t

In this section the performances of the SARIMA(2,0,1)(0,1,1)52 with drift, RWD

with weekly dummy variables and ARIMA(2,0,1)(0,1,1)52 with drift, with flu indicator

models are measured through computation of out-of-sample forecasts. The out-of-sample forecasts are based on the following data: the mortality index per week from the first week of 1995 to the last week of 2011. In other words, the full sample except for the last three years is used.

The graphs below show how well the forecasts from 2012 to 2015 (represented by the blue lines) match the observed values of κt (represented by the red lines) in the same

three year period according to the various models that were tested. Confidence intervals are not included with the forecasts in order to be able to zoom in better on the graphs. In table 4.2 the out-of-sample forecast performance of the models is compared by calculating the Mean Absolute Error (MAE) and the Root Mean Square Error (RMSE). Further in this section, out-of-sample forecasts of mtare shown. The performance of

the models is measured again using the MAE and RMSE. These values are shown in table 4.3.

At the end of the section the accuracy of the out-of-sample forecasts of the seasonal models is compared to the accuracy of out-of-sample forecasts of a yearly RWD model by aggregating the weekly forecasts into yearly ones and then calculating the MAE and RMSE (see table 4.4).

4.2.1 SARIMA(2,0,1)(0,1,1)52 with drift

Figure 4.5 shows the out-of-sample forecast of the mortality index per week using the SARIMA(2,0,1)(0,1,1)52with drift model. The forecast is a good approximation for

the actual values of κtin the summer seasons but excess mortality in winter is less well

forecasted. The height of the peaks only matches in the 2013/2014 winter season. For the remaining winter seasons, the forecast of the peaks is too low.

(32)

Figure 4.5: Three year forecast of κtusing the SARIMA(2,0,1)(0,1,1)52with drift model

(blue line) and actual values of κt from 1995 to 2012 (black line) and from 2012 to 2015

(red line).

4.2.2 RWD with weekly dummy variables

Figure 4.6: Three year forecast of κt using an ARIMA(0,1,0) with seasonal dummy

variables (blue line) and actual values of κt from 1995 to 2012 (black line) and from

2012 to 2015 (red line).

Figure 4.6 shows the out-of-sample forecast of the mortality index per week using the ARIMA(0,1,0) with weekly dummy variables model. Years 2012 to 2015 appear to be fairly well forecasted. Similarly to the SARIMA model, the summer seasons match reasonably well. What the model fails to capture are the high peaks at the end of year

(33)

2012 and 2014 as this model will always yield equal peak heights. Peaks appear to be slightly higher than in the SARIMA model.

4.2.3 SARIMA(2,0,1)(0,1,1)52 with drift, with flu indicator

Figure 4.7: Three year forecast of κt using a SARIMA(2,0,1)(0,1,1)52 with drift, with

a flu indicator as exogenous regressor (blue line) and actual values of κt from 1995 to

2012 (black line) and from 2012 to 2015 (red line).

The out-of-sample forecast of the mortality index per week using a flu indicator as exogenous regressor with a SARIMA(2,0,1)(0,1,1)52with drift model is shown in figure

4.7. The forecast yields similar results to the SARIMA(2,0,1)(0,1,1)52 with drift model

without flu indicator but the forecast of the winter peaks is improved. Especially the highest peak occurring in the 2012/2013 winter season is improved.

4.2.4 Out-of-sample forecasts of the central mortality rate

It is also interesting to look at how the out-of-sample results for κttranslate to the

forecasted central mortality rate. Using equation 3.11 and the obtained forecasts of κt

the forecasted central mortality rate can be estimated as follows:

ˆ

mT +f = exp{ ˆα + ˆκT +f}. (4.3)

The following three graphs (figures 4.8,4.9 and 4.10) show how ˆmT +f compares to

(34)

Forecasts from SARIMA(2,0,1)(0,1,1)[52] with drift

Year centr al mor tality r ate m t 2012.0 2012.5 2013.0 2013.5 2014.0 2014.5 2015.0 0.00014 0.00015 0.00016 0.00017 0.00018 0.00019 mt 2012 − 2015 forecast (out−of−sample)

Figure 4.8: Three year forecast of mtusing the SARIMA(2,0,1)(0,1,1)52with drift model

(blue line) and actual values of mt from 2012 to 2015 (black line).

Year centr al mor tality r ate m t 2012.0 2012.5 2013.0 2013.5 2014.0 2014.5 2015.0 0.00014 0.00015 0.00016 0.00017 0.00018 0.00019 mt 2012 − 2015 forecast (out−of−sample)

Figure 4.9: Three year forecast of mt using an RWD with seasonal dummy variables

(blue line) and actual values of mt from 2012 to 2015 (black line).

4.2.5 Model evaluation

By defining Yt as the observation of the mortality index at week t and Ft as the

forecast of the mortality index at week t it is possible to compare the models based on the forecast error et= Yt− Ft. According toHyndman and Koehler(2006), the MAE =

mean(|et|) and the RMSE =

p

mean(e2

t) are suitable measures for comparing different

methods applied to the same data set. Measures involving percentage errors (such as the Mean Absolute Percentage Error), on the other hand, should be avoided when there are values of Yt close to 0 because these will have a highly skewed distribution.

(35)

Forecasts from SARIMA(2,0,1)(0,1,1)[52] with drift, with Flu Indicator Year centr al mor tality r ate m t 2012.0 2012.5 2013.0 2013.5 2014.0 2014.5 2015.0 0.00014 0.00015 0.00016 0.00017 0.00018 0.00019 mt 2012 − 2015 forecast (out−of−sample)

Figure 4.10: Three year forecast of mt using a SARIMA(2,0,1)(0,1,1)52 with drift, with

a flu indicator as exogenous regressor (blue line) and actual values of mt from 2012 to

2015 (black line).

Mortality index

In table 4.2 the MAE and RMSE for the three ARIMA models used in this section are shown. The calculations include all of the out-of-sample observations and forecasts. .

Model MAE RMSE

SARIMA(2,0,1)(0,1,1)52 with drift 0.0513 0.0644

RWD with weekly dummy variables 0.0400 0.0503 SARIMA(2,0,1)(0,1,1)52 with drift, with flu indicator 0.0364 0.0463

Table 4.2: A comparison of the out-of-sample performance of the models used for fore-casting the mortality index.

From the results shown in the table it can be concluded that the ARIMA(2,0,1)(0,1,1)52

with flu indicator model yields the best results as it produces the most accurate forecast of κt for the years 2012 to 2015. Compared to the SARIMA(2,0,1)(0,1,1)52 with drift

it improves the MAE by 29% and the RMSE by 28%, which is significant. The model that performs worse is the SARIMA(2,0,1)(0,1,1)52 with drift.

Central mortality rate

Table 4.3 shows the the MAE and RMSE obtained when the models are used to forecast the central mortality rate mt:

Similarly to the forecasts of the mortality index κt, the model that performs best in

out-of-sample forecasting is the SARIMA(2,0,1)(0,1,1)52 with drift, with flu indicator

(36)

30 Emanuela Festa — Seasonal Mortality Forecasting .

Model MAE RMSE

SARIMA(2,0,1)(0,1,1)52 with drift 7.9 × 10−6 1.0 × 10−5

RWD with weekly dummy variables 6.4 × 10−6 8.2 × 10−6 SARIMA(2,0,1)(0,1,1)52 with drift, with flu indicator 5.7 × 10−6 7.4 × 10−6

Table 4.3: A comparison of the out-of-sample performance of the models used for fore-casting the central mortality rate.

A comparison between seasonal and yearly models

Finally, a comparison of the forecasting performance in terms of MAE and RMSE was made between the seasonal models with weekly data and a RWD model with annual data. In order to compare weekly and annual models, each 52 weekly out-of-sample forecasts of mt from the seasonal models were aggregated to form yearly forecasts of

mt. The results are shown in table4.4.

Model MAE RMSE

RWD 3.0 × 10−4 3.0 × 10−4

(input: annual data; forecast: annual mt)

SARIMA(2,0,1)(0,1,1)52 with drift 3.9 × 10−4 3.9 × 10−4

(input: weekly data; forecast: aggregated weekly mt)

RWD with weekly dummy variables 1.8 × 10−4 1.8 × 10−4 (input: weekly data; forecast: aggregated weekly mt)

SARIMA(2,0,1)(0,1,1)52 with drift, with flu indicator 2.0 × 10−4 2.0 × 10−4

(input: weekly data; forecast: aggregated weekly mt)

Table 4.4: A comparison of the out-of-sample forecasting performance of a yearly and three seasonal models used for forecasting the central mortality rate.

The results show that for two out of three models, the aggregated seasonal forecasts of mt are more accurate than forecasts from a standard RWD model that uses annual

data. The only model that performs worse is the SARIMA with drift model (without flu indicator). Although the sample used for these calculations is quite small, the results suggest that some seasonal ARIMA models are more accurate than the yearly RWD model.

4.3 Long-term forecasts

In this section the forecasts of κt in the long term are shown. The graphs were

obtained by extending the forecasts in the first section of this chapter from three to twenty years. Once again confidence intervals are left out in order to zoom in better on the graphs.

The SARIMA(2,0,1)(0,1,1)52 with drift and RWD with seasonal dummy variables

models in figures4.11and4.12show very similar forecasts. Both models yield very reg-ular forecasts since in each one all periods repeat exactly the same pattern. A difference between the two models is that the RWD with seasonal dummy variables model predicts slightly higher peaks. Both models project very regular forecasts, which do not reflect the great variation in peak height shown by the observed κt values from 1995 to 2015.

The flu indicator model is suitable for testing various future influenza epidemics scenarios. As an example, figure 4.13 shows a twenty year forecast made with the flu

(37)

Forecasts from SARIMA(2,0,1)(0,1,1)[52] with drift Year κt 2000 2010 2020 2030 −0.2 −0.1 0.0 0.1 0.2 0.3 _κ t forecast

Figure 4.11: Twenty year forecast of the mortality index per week using the SARIMA(2,0,1)(0,1,1)52 with drift model.

Year κt 2000 2010 2020 2030 −0.2 −0.1 0.0 0.1 0.2 0.3 _κ t forecast

Figure 4.12: Twenty year forecast of the mortality index per week using a RWD model with seasonal dummy variables.

indicator model. The flu indicator used for forecasting years 2015 to 2035 was the same as the one used for fitting the model from year 1995 to 2015. As the example of figure 4.13shows, the winter peaks follow roughly the pattern of the peaks from 1995 to 2015.

(38)

Forecasts from SARIMA(2,0,1)(0,1,1)[52] with drift, with Flu Indicator

Year κt 2000 2010 2020 2030 −0.2 −0.1 0.0 0.1 0.2 0.3 κt forecast

Figure 4.13: Twenty year forecast of the mortality index per week using a SARIMA(2,0,1)(0,1,1)52 with drift, with a flu indicator as exogenous regressor. The

flu indicator used for the forecast uses the same flu values observed from 1995 to 2015. However, it is possible to see that the accuracy of the forecast decreases with time as the winter peaks tend to become more and more regular instead of following the peak pattern of the observed time series.

(39)

Conclusions

The first conclusion that can be made from the results obtained in this research is that in order to model the seasonal component of the time-varying parameter κt, which

was determined per week, a more complex type of ARIMA structure is needed than the RWD used by Lee-Carter in 1992 to model yearly mortality. Therefore, to find a suitable model for forecasting seasonal mortality three types of seasonal ARIMA models were investigated: a SARIMA(p,d,q)(P,D,Q)m, an ARIMA with m − 1 dummy variables and

a SARIMA(p,d,q)(P,D,Q)m combined with an indicator based on influenza intensity.

After a careful model selection procedure following Box-Jenkins methodology, it was found that a SARIMA(2,0,1)(0,1,1)52and an ARIMA(5,1,1) with dummy variables can

indeed significantly approximate the seasonal behaviour belonging to the weekly mortal-ity index (and hence of the weekly mortalmortal-ity rate). When comparing the two models, it was found that the ARIMA with dummy variables can produce acceptable results even when a simple RWD form is considered (although model selection through Box-Jenkins methodology will of course improve the goodness of fit). The SARIMA(p,d,q)(P,D,Q)m,

however, needs to be more complex than a RWD SARIMA (a SARIMA(0,1,0)(0,1,0)52

with constant) in order to achieve similar results.

The out-of-sample forecasts presented in the “Results” chapter show that the models perform well in the short term since the forecasts, in general, do not differ considerably from the actual observations (except for at the “winter peaks”). The models, however, produce a very regular type of forecast that assumes that the mortality index will behave exactly the same way each year, apart from a light overall drift that tends to decrease the mortality over the years.

In order to model the great variation in excess mortality during the winter season shown by κt, a flu indicator was added to the SARIMA(2,0,1)(0,1,1)52 model. The

out-of-sample forecasts results show that the addition of this exogenous regressor in the model does improve the goodness of fit significantly (29% improvement of the MAE and 28% improvement of the RMSE with respect to the SARIMA without indicator). What the flu indicator does not model, however, is the mortality variation in the summer season. Also, the long term forecasts show that, given a specific flu indicator scenario, the ability to forecast excess mortality using the flu indicator decreases with time. So the further in the future the forecast is, the lower the influence of the flu indicator.

When out-of-sample forecasts of the central mortality rate obtained with a standard annual RWD model were compared with aggregated seasonal out-of-sample forecasts, two of the three seasonal models provided more accurate forecasts. These were the SARIMA with flu indicator and the RWD with dummy variables. This suggests that aggregating higher frequency forecasts could improve the accuracy of annual forecasts. However, more tests should be made (and with larger samples) in order to further

(40)

investigate this claim.

Influenza was chosen as a predictor for excess mortality in winter because flu epi-demics occur almost every year and studies suggest that they are very likely the singular cause of winter excess mortality in economically developed countries (seeReichert et al.

(2004)). Moreover, the histograms shown in chapter 2show that it is a better variable to use than average temperature due to its stronger correlation with the weekly mortal-ity rate (in the time period of 1995 to 2015). The results show that, if information on influenza epidemics is available and it is incorporated in the SARIMA model described in this thesis, seasonal mortality forecasts can become considerably more accurate than when influenza is not taken into account.

Although predicting the intensity of future flu epidemics is not yet possible, research on this topic is taking place in the medical field and may become feasible in the future. See for example Viboud et al. (2003), which describes a method that proved suitable for retrospectively predicting temporal and geographical diffusion of influenza up to 10 weeks in advance, and Shaman and Karspeck(2012), for more recent studies using web-based estimates of local influenza infection rates.

Another interesting development is the increasing use of flu vaccines in The Nether-lands by elderly people, which could be one of the causes of the decreasing flu epidemic intensity observed over the last decades (Dijkstra et al., 2009). If a clear correlation between flu vaccinations and excess winter mortality exists, the flu indicator model de-scribed in this thesis could be used to study the impact of different vaccination policies on seasonal mortality.

Improvements of the forecasts could be achieved once a larger sample is available. Since the current flu indicator only models winter excess mortality, it would also be a good idea to try and investigate the correlation between heat waves and summer mortality (see for exampleHuynen et al.,2001) so that an indicator for excess summer mortality can be constructed. Moreover, to improve forecasts of excess winter mortality another indicator based on temperature could be added to the model. Furthermore, the estimation of the drift ˆφ0 parameter could, perhaps, be improved by using yearly data

instead of applying the estimation procedure described in chapter3.

In order to continue the research on seasonal mortality forecasting, one could dif-ferentiate the data (if available) according to gender and age. Finally, it would also be interesting to model seasonal mortality using another frequency, for example, using monthly data instead of weekly.

Seasonal mortality forecasting based on influenca epidemics in the Netherlands

on Influenza Epidemics in The

Netherlands

Emanuela Festa

Contents

Acknowledgements

Introduction

Data description

2.1

Influenza epidemics

2.2

Average temperature

2.3

Central mortality rate

2.4

Correlation of m

with flu data and average

tempera-tures in winter

Chapter 3

Modelling of weekly Dutch

population mortality

3.1

The Lee-Carter Model

3.2

Application of the Lee-Carter Model to weekly

non-age specific mortality

3.3

Non-seasonal ARIMA models

3.4

Seasonal ARIMA structures

3.5

ARIMA model with weekly dummy variables

3.6

SARIMA with flu indicator as exogenous regressor

Chapter 4

Results

4.1

Full sample forecasts

4.2

Out-of-sample forecasts of κ

and m

4.3

Long-term forecasts

Conclusions