Example 6.1: ETS(M,N,M) model

For the ETS(M,N,M) model, θ ₁ = _ ² _n, and for h ≥ 2,

h− 1

θ_h = _ ² _n + α ² σ ²∑ θ_h−j

j =1

= _ ² _n + α ² σ ²(θ ₁ + θ ₂ + · · · + θ_{h− 1}).

6.4 Class 3: Some Nonlinear Seasonal State Space Models

Then, by induction, we can show that θ_h = _ ² _n (1 + α ² σ ²) ^{h− 1}. Plugging this into (6.5) gives the following simpler expression for v_{n + h|n}:

_ _

v_{n + h|n} ≈ s ² _{n−m + h} + _m _ ² _n (1+ σ ²)(1+ α ² σ ²) ^{h− 1}(1+ γ ² σ ²) ^hm − 1.

The expression is exact for h ≤ m.

6.4.2 Exact Forecast Means and Variances

To obtain the exact formulae for h > m, we first write the models in Class 3 using the following nonlinear state space model:

y_t = w ₁ ^_x_{t− 1} w ₂ ^_z_{t− 1}(1+ ε _t),

x_t = (F ₁+ G ₁ ε _t) x_{t− 1}, z_t = (F ₂+ G ₂ ε _t) z_{t− 1},

where F ₁, F ₂, G ₁, G ₂, w ₁ ^_ and w ₂ ^_ are all matrix or vector coefficients, and x_t and z_t are unobserved state vectors at time t. As for Class 2, {ε _t } isNID(0, σ ²), where the lower tail of the distribution is truncated so that 1 + ε _t is positive.

Let k be the length of vector x_t and q be the length of vector z_t. Then the orders of the above matrices are as follows:

F

(k

×

k)

G

(k

×

k)

w^_

(1

×

k)

(q

q)

(q

q)

(1

q)

F

×

G

×

w^_

×

• For the ETS(M,N,M) model, x_t = __t, z_t = (s_t,..., s_{t−m +1}) ^_, and the matrix coefficients are w ₁ = 1, w ₂ ^_ = [0,..., 0, 1],

F 1=1,

F 2=

0 ^_

,

G 1

= α, and G 2

=

0 ^_

γ

.

_

m− 1

m− 1

0 m− 1^_

^Im− 1

0 m− 1^_

_

^Om− 1

•

For the ETS(M,A

,M) model, xt = (_t, bt) ^_, w^_

= [1, 1],

d

F 1=^_

1 φ

_ ,

G 1=^_

α α

0 φ

β β ^_,

and z ₂, w ₂, F ₂ and G ₂ are the same as for the ETS(M,N,M) model.

• The ETS(M,A,M) model is equivalent to the ETS(M,A_d,M) model with φ =1.

	6 Prediction Distributions and Intervals
	For models in this class,
	^µn + h|n = ^w 1 ^_^Mh− 1 ^w 2	(6.7)
and

v_{n + h|n} = (1+ σ ²)(w ₂ ^_ ⊗ w ₁ ^_) V_{n + h− 1 |n} (w ₂ ^_ ⊗ w ₁ ^_) ^_ + σ ² µ ² _{n + h|n},

where ⊗ denotes a Kronecker product (Schott 2005, Sect. 8.2), M ₀ V ₀= O _{2 m} , and for h ≥ 1,

M_h = F ₁ M_{h− 1} F ₂ ^_ + G ₁ M_{h− 1} G ₂ ^_σ ²

and

(6.8)

= x_nz_n^_,

(6.9)

^Vn + h|n = (^F 2 ^{⊗ F} 1) ^Vn + h− 1 |n (^F 2 ^{⊗ F} 1) ^_

+ σ ^2 ^_(F 2 ⊗ F 1) Vn + h− 1 |n (G 2 ⊗ G 1) ^_ + (G 2 ⊗ G 1) Vn + h− 1 |n (F 2 ⊗ F 1) ^_ ^_

+

σ

(G

2 ^⊗

F + F

G

) V

MM (G

F + F

⊗

G

)

_

2 ^⊗

_

n + h− 1 |n + ^−→h− 1 ^−→h^_− 1 ^_

2 ^⊗ 1

+

σ

(G

2 ^⊗

G

)

MM (G

2 ^⊗

G

)

_ ,

(6.10)

^_3 Vn + h− 1 |n +

2 ^−→h− 1 ^−→h_− 1^_

−→

where M _{h− 1} = vec(M_{h− 1}). (That is, the columns of M_{h− 1} are stacked to form a vector.) Note, in particular, that µ_{n +1 |n} = (w ₁ ^_x_n)(w ₂ ^_z_n) and v_{n +1 |n} = σ ² µ ² _{n +1 |n}. While these expressions look complicated and providelittle insight, it is relatively easy to compute them using computer matrix languages such as R and Matlab.

In Appendix “Derivation of results for Class 3,” we show that the approx-imations (6.4) and (6.5) follow from the exact expressions (6.7) and (6.8). Note that the usual point forecasts for these models are given by (6.4) rather than (6.7).

6.4.3 The Accuracy of the Approximations

In order to investigate the accuracy of the approximations (6.4) and (6.5) for the exact mean and variance given by (6.7) and (6.8), we provide some comparisons for the ETS(M,A,M) model in Class 3.

These comparisons are done for quarterly data, where the values for the components are assumed to be the following: __n = 100, b_n = 2, s_n = 0.80, s_{n− 1}=1.20, s_{n− 2}=0.90 and s _{n− 3}=1.10. We use the following base levelvalues for the parameters: α = 0.2, β = 0.06, γ = 0.1, and σ = 0.05. We vary these parameters one at a time as shown in Table 6.5.

The results in Table 6.5 show that the mean and approximate mean are always very close, and that the percentage difference in the standard

6.4 Class 3: Some Nonlinear Seasonal State Space Models

Table 6.5. Comparison of exact and approximate means and standard deviations forthe ETS(M,A,M) model in Class 3.

Period	Exact	Approximate	Exact	Approximate	SD percent
ahead	mean (6.7)	mean (6.4)	SD (6.8)	SD (6.5)	difference
h	^µn + h|n		^√vn + h|n
σ =0.05, α =0.2, β =0.06, γ =0.1
	121.01	121.00	7.53	7.33	2.69
	100.81	100.80	6.68	6.52	2.37
	136.81	136.80	9.70	9.50	2.07
	92.81	92.80	7.06	6.93	1.80
	129.83	129.80	10.85	10.45	3.68
	108.03	108.00	9.65	9.34	3.21
	146.44	146.40	13.99	13.60	2.81
	99.22	99.20	10.13	9.88	2.47
σ =0.1, α =0.2, β =0.06, γ =0.1
	121.05	121.00	15.09	14.68	2.73
	100.84	100.80	13.39	13.07	2.40
	136.86	136.80	19.45	19.04	2.11
	92.84	92.80	14.15	13.89	1.84
	129.93	129.80	21.77	20.96	3.75
	108.11	108.00	19.39	18.75	3.29
	146.55	146.40	28.11	27.30	2.89
	99.30	99.20	20.35	19.83	2.55
σ =0.05, α =0.6, β =0.06, γ =0.1
	121.02	121.00	10.87	10.60	2.47
	100.82	100.80	9.96	9.76	2.04
	136.83	136.80	14.76	14.51	1.72
	92.82	92.80	10.86	10.70	1.47
	129.86	129.80	16.64	16.19	2.71
	108.05	108.00	14.83	14.48	2.37
	146.46	146.40	21.45	21.00	2.09
	99.24	99.20	15.45	15.16	1.86
σ =0.05, α =0.2, β =0.18, γ =0.1
	121.03	121.00	10.19	9.87	3.08
	100.82	100.80	9.88	9.66	2.27
	136.83	136.80	15.55	15.29	1.69
	92.82	92.80	12.14	11.98	1.28
	129.87	129.80	19.67	19.16	2.56
	108.06	108.00	18.41	18.04	2.03
	146.48	146.40	27.86	27.41	1.64
	99.26	99.20	20.93	20.65	1.35
σ =0.05, α =0.2, β =0.06, γ =0.3
	121.04	121.00	8.10	7.53	7.12
	100.83	100.80	7.13	6.68	6.36
	136.84	136.80	10.28	9.70	5.64
	92.83	92.80	7.42	7.05	4.97
	129.90	129.80	11.89	10.77	9.46
	108.08	108.00	10.47	9.59	8.42
	146.51	146.40	15.04	13.91	7.49
	99.27	99.20	10.79	10.07	6.67

88 6 Prediction Distributions and Intervals

deviations only becomes substantial when we increase γ. This result for the standard deviation is not surprising because the approximation is exact if γ =0. In fact, we recommend that the approximation not be used if thesmoothing parameter for γ exceeds 0.10.

6.5 Prediction Intervals

The prediction distributions for Class 1 are clearly Gaussian, as the models are linear and the errors are Gaussian. Consequently, 100(1 − α)% prediction intervals can be calculated from the forecast means and variances in the usual

way, namely µ_{n + h|n} ± z_{α /2} ^√v_{n + h|n}, where z_q denotes the q th quantile of a standard Gaussian distribution.

In applying these formulae, the maximum likelihood estimator for σ ² (see

p. 68) is simply

n

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

t =1

where ε ^ˆ t = yt − µt|t− 1.

The prediction distributions for Classes 2 and 3 are non-Gaussian because of the nonlinearity of the state space equations. However, prediction inter-vals based on the above (Gaussian) formula will usually give reason-ably accurate results, as the following example shows. In cases where the Gaussian approximation may be unreasonable, it is necessary to use the simulation approach of Sect. 6.1.

6.5.1 Application: Quarterly French Exports

As a numerical example, we consider the quarterly French exports data given in Fig. 6.1, and use the ETS(M,A,M) model. We estimate the parameters to be α = 0.8185, β = 0.01, γ = 0.01 and σ = 0.0352, with the final states __n = 757.3, b_n = 15.7, and z_n = (0.873, 1.141, 1.022, 0.964) ^_.

Figure 6.4 shows the forecast standard deviations calculated exactly using (6.8) and approximately using (6.5). The approximate values are so close to the exact values in this case (because σ ² and γ are both very small) that it is almost impossible to distinguish the two lines.

The data with three years of forecasts are shown in Fig. 6.5. In this case, the conditional mean forecasts obtained from model ETS(M,A,M) are virtually indistinguishable from the usual forecasts because σ is so small (they are identical up to h = m). The solid lines show prediction intervals calculated as µ_{n + h |n} ± 1.96 ^√v_{n + h |n}, and the dotted lines show prediction intervals computed by generating 20,000 future sample paths from the fit-ted model and finding the 2.5 and 97.5% quantiles at each forecast horizon.

					6.5	Prediction Intervals
		Exact
	Small sigma approximation
deviation
standard
Forecast

Forecast horizon

Fig. 6.4. Forecast standard deviations calculated (a) exactly using (6.8); and (b)approximately using (6.5).

	Percentile−based interval
Variance−based interval
				Quarter

Fig. 6.5. Quarterly French exports data with 3 years of forecasts. The solid lines showprediction intervals calculated as µ_{n + h|n} ± 1.96 ^√v_{n + h|n}, and the dotted lines show pre-diction intervals computed by generating 20,000 future sample paths from the fitted model and finding the 2.5 and 97.5% quantiles at each forecast horizon.

90 6 Prediction Distributions and Intervals

Clearly, the variance-based intervals are a good approximation despite the non-Gaussianity of the prediction distributions.

6.6 Lead-Time Demand Forecasts for Linear Homoscedastic Models

For Class 1 models, it is also possible to obtain some analytical results on the distribution of lead-time demand, defined by

h
Yn (h) =∑ yn + j.	(6.11)

j =1

In particular, the variance of lead-time demand can be used when imple-menting an inventory strategy, although the basic exponential smoothing procedures originally provided only point forecasts, and rather ad hoc formulae were the vogue in inventory control software.

Harrison (1967) and Johnston and Harrison (1986) derived a variance for-mula for lead-time demand based on simple exponential smoothing using a state space model with independent error terms. They utilized the fact that simple exponential smoothing emerges as the steady state form of the model predictions in large samples. Adopting a different model, Snyder et al. (1999) were able to obtain the same formula without recourse to a restrictive large sample assumption. Around the same time, Graves (1999) obtained the formula using an ARIMA(0,1,1) model.

Harrison (1967) and Johnston and Harrison (1986) also obtained a variance formula for lead-time demand when trend-corrected exponential smoothing is employed. Yar and Chatfield (1990), however, suggested a slightly different formula. They also provide a formula that incorporates seasonal effects for use with the additive Holt-Winters method.

The approach we adopt here is based on Snyder et al. (2004), although the parameterization in this book is slightly different from that used in Snyder et al. (2004). The results obtained subsume those in Harrison (1967), Johnston and Harrison (1986), Yar and Chatfield (1990), Graves (1999) and Snyder et al. (1999). In addition, for ETS(A,A,A), the recursive variance formula in Yar and Chatfield (1990) has been replaced with a closed-form counterpart.

6.6.1 Means and Variances of Lead-Time Demand

In Appendix “Derivation of C_j values” we show that

j− 1

^yn + j ^{= µ}n + j|n ⁺∑ ^cj−i^ε n + i ^{+ εq}n + j ^, i =1

6.6 Lead-Time Demand Forecasts for Linear Homoscedastic Models

where µ

n + j|n

and ck are given in Table 6.2. Substitute this into (6.11) to give

h

_

j− 1

h

h

Yn (h) =∑

^µn + j

n ^+∑ ^cj

−

i ^ε n + i

^+ ^ε n + j

^= ∑

^µn + j

n ^+∑ ^Cj

−

1 ^ε n + h

−

j +1^,

j =1

|

i =1

_

j =1

|

j =1

(6.12)

where

j

j =1,..., h − 1.

C 0=1

and

Cj =1+∑ ci

for

(6.13)

i =1

Thus, lead-time demand can be resolved into a linear function of the uncorrelated level and error components.

From (6.12), it is easy to see that the point forecast (conditional mean) is simply

h
ˆ	(6.14)
Yn (h) =E(Yn (h) | ^xn) =∑ µn + j|n
j =1
and the conditional variance is given by
h− 1
V (Yn (h) | xn) = σ ^2 ∑ Cj ^2.	(6.15)

j =0

The value of C_j for each of the models is given in Table 6.6. These expressions are derived in Appendix “Derivation of C_j values.”

As with the equations for forecast variance at a specific forecast horizon, we can substitute these expressions into (6.15) to derive a specific formula for each model. This leads to a lot of tedious algebra that is of limited value. Therefore we only give the result for model ETS(A,N,N):

Table 6.6. Values of C_j to be used in computing the lead-time variance in (6.15).

Model	Cj
(A,N,N)		+ jα			_
		_
(A,A,N)		+ j	α +^1 β (j +1)
(A,Ad,N)		+ jα +	βφ	^_(j + 1)(1 − φ) − (1 − φ^j ^+^1)^_
(1 −φ)^2
(A,N,A)		+ jα + γ jm	_
		_
(A,A,A)		+ j	α +2^1 β (j +1) + γ jm
(A,Ad,A)		+ jα +	βφ	^_(j + 1)(1 − φ) − (1 − φ^j ^+^1)^_ + γ jm
(1 −φ)^2

Here m is the number of periods in each season and j_m = j / m_ is the number of complete seasonal cycles that occur within j time periods.

	6 Prediction Distributions and Intervals
	h− 1
	V(Yn (h) | xn) = ∑(1 + jα)^2
	j =0	1 + α (h − 1) + 6^1 α ^2 (h − 1)(2 h − 1)^_.
	= σ 2 h ^_	(6.16)

6.6.2 Matrix Calculation of Means and Variances

The mean and variance of the lead-time demand, and the forecast mean and variance for a single period, can also be computed recursively using matrix equations. From Chap. 3, we know that the form of the Class 1 models is

y_t = w^_x_{t− 1}+ ε _t,

x_t = F x_{t− 1}+ gε _t,

where w^_ is a row vector, g is a column vector, F is a matrix, x_t is the unobserved state vector at time t, and {ε _t } is NID(0, σ ²).

Observe that the lead-time demand can be determined recursively by

Yn (j) = Yn (j − 1) + yn + j,

(6.17)

where Y_n (0) = 0 and Y_n (j) = ∑ _i^j ₌₁ y_{n + i}. Consequently, (6.17) can be written as

Yn (j) = Yn (j − 1) + w^_xn + j− 1+ ε n + j.

(6.18)

So, if the state vector x_{n + j} is augmented with Y_n (j), the first-order recurrence relationship

^xn + j	=	F 0	^xn + j− 1	+ ^g		ε	n + j
^_ Yn (j)^_ _ w_ 1	^_ _ Yn (j − 1)^_ _1	_
is obtained. This has the general form zn + j = Azn + j	− 1	+ bε n + j. If the mean
			z				z	=
and variance of the zn + j are denoted by mn + j|n	= E(zn + j |	^xn)and ^Vn + j|n

V(z_{n + j} | x_n), then they can be computed recursively using the equations

m^z	= Am^z	,
n + j|n	n + j− 1 |n
Vn^z + j|n = AVn^z + j− 1 |nA^_ + σ ^2 bb^_.
The mean of the lead-time demand Yn (h) is the last element in m^z	, and
the variance of Yn (h) is the bottom right element of V ^z	n + h|n
.
		n + h|n

This same procedure of using an augmented matrix can also be applied to find the forecast mean and variance of y_{n + h} for any single future time period t = n + h. In this case, the state vector x_{n + j} is augmented with y_{n + j} in placeof Y_n (j), and

= ^_ F 0^_

^A _w_ ₀^.

6.6 Lead-Time Demand Forecasts for Linear Homoscedastic Models

Then, the mean and variance of y^{n + h} are the last elements in m^{z + |} and

n h n

^{V z + |} respectively. Of course, one can use A = [ F, w^_ ] ^_ and the general

n h n

form z_{n + j} = Ax_{n + j− 1} + bε _{n + j} to remove the unnecessary multiplications by 0 in an actual implementation.

6.6.3 Stochastic Lead-Times

In practice, lead-times are often stochastic, depending on various factors including demand in the previous time periods. We explore the effect of stochastic lead-times on forecast variances in the case of the ETS(A,N,N) model for simple exponential smoothing.

Let the lead-time, T, be stochastic with mean E(T) = h. The mean lead-time demand, given the level at time n, is

E(Y_n (T) | __n) = E _T [E(Y_n (T) | T, __n)] = h__n,

as in the case of a fixed lead-time. The variance of the lead-time demand reduces to

V(Y_n (T) | __n) = V _T [E(Y_n (T) | T, __n)] + E _T [V(Y_n (T) | T, __n)]

	T
= V T (_n T) + E T	σ ^2∑ Cj ^2, T
_	j =1	_
	T							_
= _ ^2 n V(T) + σ ^2E T	∑		+ 2 α (T	−	j) + α ^2(T	−	j)^2
^_ j =1	_			_
= _ ^2 n V(T) + σ ^2 h + σ ^2 α ^_(1 + 2^1 α) h [2] + 3^1 αh [3] ^_,
where h [ j ] = E[ T (T− 1)... (T − j + 1)], j	= 1, 2,..., is known as the j th
factorial moment of the distribution of T.

For example, when the lead-time is fixed, h _{[ j ]} = h (h − 1)... (h − j + 1).

When the lead-time is Poisson with mean h, then h _{[ j ]} = h^j. Therefore, the lead-time demand variance becomes

_ _

V(Y_n (T) | __n) = (_ ² _n + σ ²) h + σ ² α (1 + ¹₂ α) h ² + ¹₃ αh ³.

Compare this with the variance for a fixed lead-time as given in (6.16). The two variances are plotted in Fig. 6.6 for α = 0.1, σ = 1 and __n = 2, showing that a stochastic lead-time can substantially increase the lead-time demand variance.

94 6 Prediction Distributions and Intervals

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