Parsimonious Seasonal Model

As noted in Sect. 3.5.5, series that are observed at many times within the main period, such as weeks in a year, require a very large number of start-ing values. Further, only a fairly short series may be available for estimation. However, many such series have a simple state space structure, such as a nor-mal level of sales outside certain special periods. The details for the additive error versions of such models are given in Sect. 3.5.5. The multiplicative error versions follow, using the same modification as elsewhere in this chapter.

66 4 Nonlinear and Heteroscedastic Innovations State Space Models

4.4.5 Other Heteroscedastic Models

The variance structure may be modified in several ways. Perhaps the sim-plest way is to incorporate an additional exponent in a manner reminiscent of the method employed by Box and Cox (1964) but without a transformation. We use the local trend model, ETS(M,A,N) in (4.4) by way of illustration. We separate out the error terms and modify them by using a power of the trend term, 0 ≤ θ ≤ 1:

y_t = (__{t− 1}+ b_{t− 1}) + (__{t− 1}+ b_{t− 1}) ^θ ε _t, __t = (__{t− 1} + b_{t− 1}) + α (__{t− 1} + b_{t− 1}) ^θ ε _t,

b_t = b_{t− 1}+ β (__{t− 1}+ b_{t− 1}) ^θ ε _t.

For example, θ = 1/3 would produce a variance proportional to the 2/3rds power of the mean, much as the cube-root transformation does. The present formulation enables us to retain the linear structure for the expectation, which in many ways is more plausible than the transformation.

4.5 Exercises

Exercise 4.1. Verify the variance expressions for ETS(M,N,N) given inSect. 4.2.1.

Exercise 4.2. Use the approach of Sect. 4.3.1 to derive the recursive relationsfor the state vector of the ETS(A,M,M) model, given in Table 2.2. Extend the argument to include the ETS(A,M_d,M) model.

Exercise 4.3. Use the approach of Sect. 4.3.1 to derive the recursive relationsfor the state vector of the ETS(M,M,M) model, given in Table 2.2. Extend the argument to include the ETS(M,M_d,M) model.

Exercise 4.4. Evaluate the h -step-ahead forecast mean squared error for thelocal trend model and for the local level model with drift, given in Sect. 4.4.1. Compare the two for various combinations of h, α and β.

Exercise 4.5. Evaluate the h -step-ahead forecast mean squared error for thedamped trend model ETS(M,A_d,N) and compare it with that for the local trend model for various combinations of h, φ, α and β.

Exercise 4.6. Show that the stability conditions for the heteroscedastic modelare exactly those for the ETS(A,A,N) model.

Exercise 4.7. The data setdjiclosecontains the closing prices for the DowJones Index for the first day of each month from October 1928 to December 2007, along with the monthly returns for that series. Fit a heteroscedastic ETS(M,A,N) model to these data for a selected part of the series. Compare your results with the random walk with drift model for the returns series.

Estimation of Innovations State

Space Models

For any innovations state space model, the initial (seed) states and the parameters are usually unknown, and therefore must be estimated. This can be done using maximum likelihood estimation, based on the innovations representation of the probability density function.

In Chap. 3 we outlined transformations (referred to as “general expo-nential smoothing”) that convert a linear time series of mutually dependent random variables into an innovations series of independent and identically distributed random variables. In the heteroscedastic and nonlinear cases, such a representation remains a viable approximation in most circumstances, an issue to which we return in Chap. 15. These innovations can be used to compute the likelihood, which is then optimized with respect to the seed states and the parameters. We introduce the basic methodology in Sect. 5.1. The estimation procedures discussed in this chapter assume a finite start-up; consideration of the infinite start-up case is deferred until Chap. 12.

Any numerical optimization procedure used for this task typically requires starting values for the quantities that are to be estimated. An appro-priate choice of starting values is important. The likelihood function may not be unimodal, so a poor choice of starting values can result in sub-optimal estimates. Good starting values (i.e., values that are as close as possible to the optimal estimates) not only increase the chances of finding the true optimum, but typically reduce the computational loads required during the search for the optimum solution. In Sect. 5.2 we will discuss plausible heuristics for determining the starting values.

5.1 Maximum Likelihood Estimation

The unknown model parameters and states must be estimated. Maximum likelihood (ML) estimators are sought because they are consistent and asymptotically efficient under reasonable conditions; for a general discus-sion see Gallant (1987, pp. 357–391). Hamilton (1994, pp. 133–149) derives

68 5 Estimation of Innovations State Space Models

the convergence properties of various numerical algorithms for computing ML estimates.

The likelihood function is based on the density of the series vector y. It is a function of a p -vector θ of parameters such as the smoothing parameters and damping factors. The likelihood also depends on the innovations variance σ ², but for reasons that will become clearer soon, it is convenient to separateit from the main parameter vector θ. Finally, the likelihood depends on the k -vector x ₀of seed states.

Under the more traditional assumptions employed in time series anal-ysis, the generating process is presumed to have operated for an extensive period of time prior to the period of the first observation, in which case the seed state vector must be random. A likelihood function must only be based on observable random quantities; unobserved random variables must be aver-aged away. We sidestep the need for averaging, and hence simplify the task of forming the likelihood function, by assuming that the process has had no life prior to period 1, in which case the seed state vector x ₀ is fixed and may be treated as a vector of parameters. The case of random seed states will be considered in Chap. 12.

It was shown in Chap. 3 that any time series {y_t } governed by a linear state space model with Gaussian innovations has a multivariate Gaussian distribution (3.2). In Sect. 4.1, the same basic result was derived as an approx-imation for the nonlinear version of the model, a remarkable conclusion that depends critically on the assumption of a fixed seed state. In essence, the joint density of the series was shown to be the weighted product of the densities of the individual innovations:

n

p (y | θ, x 0, σ ^2) =∏ p (ε t)/ |r (xt− 1) |.

t =1

So the Gaussian likelihood can be written as

, σ ^2 y) = (2π σ ^2) ^−^n ^/2

n

_

− 1

n

ε ^2

/ σ ^2

(θ, x

∏

r (x

t

)

exp

∑

,(5.1)

L

|

−

−

t

^_ t =1

t =1

_

_

_

_

and the log likelihood is

_

_

n

n

log L = − ^n 2 log(2π σ ^2) − ∑ log |r (xt− 1) | − ^12 ∑ ε ^2 t / σ ^2.

(5.2)

t =1

t =1

Then taking the partial derivative with respect to σ ² and setting it to zero gives the maximum likelihood estimate of the innovations variance σ ² as

n

σ ˆ²= n^{− 1}∑ ε ² _t. t =1

This formula can be used to eliminate σ ² from the likelihood (5.1) to give the concentrated likelihood

5.1

Maximum Likelihood Estimation

n

− 1

L

(θ, x 0

|

y) = (2π e σ ˆ^2) ^−n ^/2

∏ r (xt

−

1).

^_ t =1

_

_

_

_

_

_

_

Thus, twice the negative log-likelihood is given by

n

− 2 log L (θ, x 0 | y) = n log(2π e σ ˆ^2) +2∑log |r (xt− 1) |

t =1

n

n

= cn + n log

∑

ε ^2 t

+ 2

∑ log

|

r (xt

−

1)

,

t =1

t =1

|

where c_n is a constant depending on n but not on θ or x ₀. Hence, maximum likelihood estimates of the parameters can be obtained by minimizing

		n		n
L	^∗ (θ, x 0) = n log	∑	ε ^2 t	+ 2 ∑ log	r (xt	−	1)	.	(5.3)
	t =1		t =1	|		|

Equivalently, they can be obtained by squared errors criterion:

S (θ, x ₀) =^_exp(L^∗ (θ, x ₀))^{_1/ n}

minimizing the augmented

	n			_	2/ n	n
=	∏ r (xt	−	1)		∑ ε ^2 t.
^_ t =1			t =1
_				_
_				_
_				_

sum of

(5.4)

In homoscedastic cases, r (x_{t− 1}) = 1 and (5.4) reduces to the traditional sum of squared errors.

Use of (5.4) criterion in place of the likelihood function means that the optimizer does not directly select the best value of σ ². The number of vari-ables being directly optimized is reduced by one, with consequent savings in computational loads. More importantly, however, it avoids a problem that sometimes arises with the likelihood function, when the optimizer chooses a trial value of the variance that is quite out-of-kilter with the errors, and a consequent numerical stability issue occurs.

For particular values of the parameters and seed states, the value of the innovation is found with ε _t = [ y_t − w (x_{t− 1})]/ r (x_{t− 1}). The state is revised with the transition

x_t = f (x_{t− 1}) + g (x_{t− 1}) ε _t.

In the case of homoscedastic errors and linear functional relationships, as assumed in Chap. 3, this simplifies to

ε t = yt − w xt− 1,	(5.5)
xt = F xt− 1+ gε t,	(5.6)

which is the general linear form of exponential smoothing (Box et al. 1994). Although expression (5.4) was derived from the likelihood (5.1), we could

start the whole process by directly specifying that the objective is to mini-mize S (θ, x ₀). This approach, known as the Augmented Least Squares (ALS)

70 5 Estimation of Innovations State Space Models

per capita	6500 7000
Quarterly GDP	5500 6000
				Year

Fig. 5.1. Plot of Australian quarterly gross domestic product per capita from theSeptember quarter of 1971 to the March quarter of 1998.

method, does not require us to make any assumptions about the distri-butional form of the errors. More generally, when the ML estimates are computed from (5.4) without any assumption of Gaussianity, we refer to the results as quasi-maximum likelihood estimators (Hamilton 1994, p. 126). Such estimators are often consistent, but the expressions for the standard errors of the estimators may be biased, even asymptotically.

5.1.1 Application: Australian GDP

To illustrate the method, we use the Australian quarterly real gross domestic product per capita¹ from the September quarter of 1971 to the March quar-ter of 1998. The deseasonalized series, which consists of 107 observations, is depicted in Fig. 5.1. We fitted a local linear trend model (3.13):

yt = _t− 1+ bt− 1	+ ε t	(5.7a)
_t = _t− 1+ bt− 1	+ αε t	(5.7b)
bt = bt− 1+ βε t		(5.7c)

by minimizing (5.4).

The results obtained depend on the constraints imposed on the parameter values during estimation. As seen in Sect. 3.4.2, the constraints 0 < α < 1 and 0 < β < α are typically imposed. They ensure that the states can be

¹ Source: Australian Bureau of Statistics.

5.2 A Heuristic Approach to Estimation

Table 5.1. Maximum likelihood results: the local trend model applied to theAustralian quarterly gross domestic product.

	Constraints
	Conventional	Stable
α	1.00	0.61
β	1.00	2.55
_ 0	4571.3	4568.7
b 0	36.5	35.1
MSE
MAPE	0.36%	0.24%

interpreted as averages. However, another set of constraints arises from the stability criterion (Sect. 3.3.1, p. 36). They ensure that the observations have a diminishing effect as they get older. This second set of constraints (derived in Chap. 10) is α ≥ 0, β ≥ 0 and 2 α + β ≤ 4. It contains the first constraint set and is much larger.

The results are summarized in Table 5.1. The parameter estimates with the conventional constraints lie on the boundary of the parameter space. In the second case, the estimates lie in the interior of the parameter space.

The MSE is more than halved by relaxing the conventional constraints to the stability conditions. The MAPEs indicate that both approaches provide local linear trends with a remarkably good fit. The lower MAPE of 0.24% for the stability approach is consistent with the MSE results.

The optimal value of 2.55 for β in the second estimation may seem quite high. It ensures, however, that the growth rate is very responsive to unantic-ipated changes in the series. A plot of the estimated growth rates is shown in Fig. 5.2. The effect is to ensure that the local trend adapts quickly to changes in the direction of the series values.

5.2 A Heuristic Approach to Estimation

The method of estimation described in the previous section seeks values of the seed vector x ₀ and the parameters θ that jointly minimize the aug-mented sum of squared errors. The inclusion of the seed state vector can be a source of relatively high computational loads. For example, if it is applied to a weekly time series, a linear trend with seasonal effects has 54 seed states, but only three smoothing parameters. In such cases, it is common practice to approximate the seed values using a heuristic method, and simply minimize with respect to the parameters θ alone.

Heuristic methods are typically devised on a case by case basis. For exam-ple, the seed level in a local level model is often approximated by the first value of a series, or sometimes by a simple average of the first few values of a series (see Makridakis et al. 1998).

72 5 Estimation of Innovations State Space Models

Growth rate
	−50
	−100
			Year

Fig. 5.2. Estimation of the local trend model for the Australian quarterly gross domes-tic product per capita from the September quarter of 1971 to the March quarter of 1998: growth rate estimates.

In Sect. 2.6.1 we described a heuristic method that works well for almost all series. For non-seasonal data, we fit a straight line a + bt to the first ten observations and set _ ₀ = a. We use b ₀ = b when the model assumes an additive trend, and for a model with a multiplicative trend we set b ₀ = 1 + b / a. For seasonal data, the seed seasonal components are obtained using aclassical decomposition method (Makridakis et al. 1998) applied to the first few years of data. Then, as for non-seasonal data, we fit a straight line to the first ten deseasonalized observations to obtain _ ₀ and b ₀.

These are simply starting values for the full optimization. We do not gen-erally use them in forecasting, unless the seasonal period m is very large (as with weekly data).

Such heuristics can be very successful, but they are not failsafe. The expo-nential decay of the effect of the seed vector in (3.6) is supposed to ensure that the effect on the sum of squared errors of any additional error introduced by the use of a heuristic is negligible. This is only true, however, when the smoothing parameters are relatively large, and so the states change rapidly over time. Then the structure that prevailed at the start of the series has little impact on the series close to the end. However, when the smoothing param-eters are small, the end states are unlikely to be very different from those at the start of the observational period. In this case, early seed values are not discounted heavily, so the effects of any extra error arising from a heuristic are unlikely to disappear. Heuristics may potentially lead to poor results in some circumstances. Hyndman et al. (2002) report that optimization of the

5.3 Exercises

seed states improves forecasts of the M3 series by about 2.8%. It appears that full optimization is to be recommended where it is practicable.

When heuristics are used, they should be designed to apply to the model when the smoothing parameters are small. Fortunately, this is not difficult to achieve. For example, as α and β approach zero, model (5.4) reduces to the global linear trend model, so that our heuristic based on fitting the linear trend to the first ten observations can be expected to perform well in such circumstances. When α and β are large, the heuristic will perform less well, but the initial conditions are then discounted more rapidly, so the effect is reduced.

Heuristics originated in an era when computers were very much slower than those available today. They were devised to avoid what were then quite daunting computational loads. Nowadays, full optimization of the likeli-hood can be undertaken in a fraction of the time. For example, Hyndman et al. (2002) found that it took only 16 min to fit a collection of models, similar to those from Tables 2.2 and 2.3, to the 3,003 time series from the M3 fore-casting competition. Much faster times should now be possible with modern computers.

It may be argued that the above example with weekly seasonal effects suggests that there remain cases where full optimization may still not be practicable. However, in this type of situation, it is advisable to reduce the number of optimizable variables. For example, weeks with similar seasonal characteristics could be grouped to reduce the number of seasonal indexes. Alternatively, Fourier representations based on sine and cosine functions, might be used. Models with a reduced number of states are likely to yield more robust forecasts.

Heuristics, however, still have a useful place in forecasting. For example, they can be used to provide starting values for the optimizer. This way the seed values are still optimized, but the optimizer begins its search from a point that is likely to be closer to the optimal solution. An advantage of this approach is that it reduces the chances of sub-optimal solutions with multi-modal augmented sum of squared errors functions. Moreover, with series containing slow changes to states, the optimizer automatically prevents any deleterious effects of poor starting values, should they occur with the use of a heuristic. Finally, the time required for optimization is typically shortened by the use of heuristics.

5.3 Exercises

The following exercises should be completed using the quarterly US gross domestic product series available in the data set usgdp.

Exercise 5.1. Fit the local level model ETS(A,N,N) to the data. This willrequire calculating the sum of squared errors criterion and minimizing it

74 5 Estimation of Innovations State Space Models

with respect to the value of _ ₀ and α. Do the estimation using your own R code. Then compare the results with those obtained using the ets() function in the forecast package for R.

Exercise 5.2. Fit the local level model with drift (Sect. 3.5.2) to the log-transformed data.

Exercise 5.3. The multiplicative version of a local level model with driftmodel is

y_t = __{t− 1} b (1+ ε _t) __t = __{t− 1} b (1 + αε _t)

where b is a multiplicative drift term. Fit this model to the raw data using the augmented sum of squared errors criterion. Contrast the results with those from Exercise 5.2.

Prediction Distributions

and Intervals

Point forecasts for each of the state space models were given in Table 2.1 (p. 18). It is also useful to compute the associated prediction distributions and prediction intervals for each model. In this chapter, we discuss how to compute these distributions and intervals.

There are several sources of uncertainty when forecasting a future value of a time series (Chatfield 1993):

1. The uncertainty in model choice—maybe another model is correct, or maybe none of the candidate models is correct.

2. The uncertainty in the future innovations ε _{n +1},..., ε _{n + h}.

3. The uncertainty in the estimates of the parameters: α, β, γ, φ and x ₀.

Ideally, the prediction distribution and intervals should take all of these into account. However, this is a difficult problem, and in most time series analysis only the uncertainty in the future innovations is taken into account.

If we assume that the model and its parameters (including x ₀) are known, then we also know x_n, the state vector at the last period of observation, because the error in the transition equation can be calculated from the obser-vations up to time n. Consequently, we define the prediction distribution as the distribution of a future value of the series given the model, its estimated parameters, and x_n. A short-hand way of writing this is y_{n + h|n} ≡ y_{n + h} | x_n.

We briefly discuss how to allow for parameter estimation uncertainty in Sect. 6.1. We do not address how to allow for model uncertainty, although this is an important issue. Hyndman (2001) showed that model uncertainty is likely to be a much bigger source of error than parameter uncertainty.

The mean of the prediction distribution is called the forecast mean and is denoted by µ_{n + h|n} = E(y_{n + h} | x_n). The corresponding forecast variance is given by v_{n + h|n} = V(y_{n + h} | x_n). We will find expressions for these quantities for many of the models discussed in this book.

76 6 Prediction Distributions and Intervals

We are also interested in “lead-time demand” forecasting, where we pre-dict the aggregate of the next h observations rather than each of the next h observations individually. We discuss this briefly here and in more detail in Chap. 18.

The most direct method of obtaining prediction distributions is to simu-late many possible future sample paths from the fitted model, and to estimate the distributions from the simulated data. This approach will work for any time series model, including all of the models discussed in this book. We describe the simulation method in more detail in Sect. 6.1.

While the simulation approach is simple and can be applied to any well-specified time series model, the computations can be time-consuming. Furthermore, the resulting prediction intervals are only available numeri-cally rather than algebraically. Therefore, the approach does not allow for algebraic analysis of the prediction distributions.

An alternative approach is to derive the distributions analytically. Ana-lytical results on prediction distributions can provide additional insight and can be much quicker to compute. These results are relatively easy to derive for some models (particularly the linear models), but very difficult for others. In fact, there are analytical results on prediction distributions for only 15 of the 30 models in our exponential smoothing framework.

When discussing the analytical prediction distributions, it is helpful to divide the thirty state space models given in Tables 2.2 and 2.3 (pp. 21–22) into five classes; Classes 1–4 are shown in Table 6.1.

For each of Classes 1–3, we give expressions for the forecast means and variances. Class 1 consists of the linear models with homoscedastic errors; these are discussed in Sect. 6.2. In Sect. 6.3 we discuss Class 2, which contains the linear models with heteroscedastic errors. Class 3 models are discussed

Table 6.1. The models separated in the exponential smoothing framework split intoClasses 1–5.

Class 1 −→	A,N,N	A,N,A
A,A,N	A,A,A
	A,Ad,N A,Ad,A
Class 2 −→	M,N,N	M,N,A	M,N,M	←− Class 3
M,A,N	M,A,A	M,A,M
	M,Ad,N M,Ad,A M,Ad,M
Class 4 −→	M,M,N		M,M,M
M,Md,N		M,Md,M
Class 5 −→		M,M,A	A,N,M	A,M,N	A,Md,N
	M,Md,A	A,A,M	A,M,A	A,Md,A
			A,Ad,M	A,M,M	A,Md,M

6.1 Simulated Prediction Distributions and Intervals

in Sect. 6.4; these are the models with multiplicative errors and multiplicative seasonality but additive trend.

Class 4 consists of the models with multiplicative errors, multiplicative trend, and either no seasonality or multiplicative seasonality. For Class 4, there are no available analytical expressions for forecast means or variances, and so we recommend using simulation to find prediction intervals.

The remaining 11 models are in Class 5. For these models, we also rec-ommend using simulation to obtain prediction intervals. However, Class 5 models are those that can occasionally lead to numerical difficulties with very long forecast horizons. Specifically, the forecast variances are infinite, although this does not usually matter in practice for short- or medium-term forecasts. This issue is explored in Chap. 15.

Section 6.5 discusses the use of the forecast mean and variance formulae to construct prediction intervals even in cases where the prediction distribu-tions are not Gaussian. In Sect. 6.6, we discuss lead-time demand forecasting for Class 1 models.

Most of the results in this chapter are based on Hyndman et al. (2005) and Snyder et al. (2004), although we use a slightly different parameterization in this book, and we extend the results in some new directions.

To simplify some of the expressions, we introduce the following notation:

h = mh_m + h ⁺ _m,

where¹ h is the forecast horizon, m is the number of periods in each season, h_m =(h − 1)/ m_ and h ⁺ _m =^_(h − 1)mod m ^_+1. In other words, h _m is thenumber of complete years in the forecast period prior to time h, and h ⁺ _m is the number of remaining times in the forecast period up to and including time h. Thus, h ⁺ _m can take values 1, 2,..., m.

6.1 Simulated Prediction Distributions and Intervals

Recall from Chap. 4 that the general model with state vector

x_t = (__t, b_t, s_t, s_{t− 1},..., s_{t−m +1}) ^_

has the form

y_t = w (x_{t− 1}) + r (x_{t− 1}) ε _t, x_t = f (x_{t− 1}) + g (x_{t− 1}) ε _t,

where w (·) and r (·) are scalar functions, f (·) and g (·) are vector functions, and {ε _t } is a white noise process with variance σ ².

One simple approach to obtaining the prediction distribution is to sim-ulate sample paths from the models, conditional on the final state x_n. This

¹ The notation u_ means the integer part of u.

78 6 Prediction Distributions and Intervals

francs)
sales (thousands of
Quarterly	Historical data			Simulated future sample paths
			Year

Fig. 6.1. Quarterly French exports data with 20 simulated future sample paths gener-ated using the ETS(M,A,M) model assuming Gaussian innovations. The solid vertical line on the right shows a 90% prediction interval for the 16-step forecast horizon,calculated from the 0.05 and 0.95 quantiles of the 5,000 simulated values.

was the approach taken by Ord et al. (1997) and Hyndman et al. (2002). That

is, we generate observations {y ⁽ _t^{i )} }, for t = n + 1,..., n + h, starting with x_n from the fitted model. Each ε _t value is obtained from a random numbergenerator assuming a Gaussian or other appropriate distribution. This pro-cedure is repeated for i = 1,..., M, where M is a large integer. (In practice, we often use M = 5,000.)

Figure 6.1 shows a series of quarterly exports of a French company (in thousands of francs) taken from Makridakis et al. (1998, p. 162). We fit an ETS(M,A,M) model to the data. Then the model is used to simulate 5,000 future sample paths of the data. Twenty of these sample paths are shown in Fig. 6.1.

Characteristics of the prediction distribution of y_{n + h|n} can then be esti-mated from the simulated values at a specific forecast horizon: y_{n + h | n} =

{y ⁽ _n ¹₊⁾ _h,..., y ⁽ _n^M _{+ h} ⁾ }. For example, prediction intervals can be obtained usingquantiles of the simulated sample paths. An approximate 100(1 − α)% pre-diction interval for forecast horizon h is given by the α /2 and 1 − α /2 quantiles of y_{n + h|n}. The solid vertical line on the right of Fig. 6.1 is a 90% prediction interval computed in this way from the 0.05 and 0.95 quantiles of the simulated values at the 16-step horizon.

The full prediction density can be estimated using a kernel density esti-mator (Silverman 1986) applied to y_{n + h|n}. Figure 6.2 shows the prediction

6.1 Simulated Prediction Distributions and Intervals

Дата добавления: 2015-10-24