Processing math: 100%
  • Section 8.8 Problem#6
  • Section 8.8 Problem#7
  • Section 8.8 Problem#8

Section 8.8 Problem#6

6. Forecast the Chinese GDP from the global_economy data set using an ETS model. Experiment with the various options in the ETS() function to see how much the forecasts change with damped trend, or with a Box-Cox transformation. Try to develop an intuition of what each is doing to the forecasts.

[Hint: use a relatively large value of h when forecasting, so you can clearly see the differences between the various options when plotting the forecasts.]

china <- global_economy |>
  filter(Country == "China")
china |> autoplot(GDP)

It clearly needs a relatively strong transformation due to the increasing variance.

china |> autoplot(box_cox(GDP, 0.2))

china |> features(GDP, guerrero)
ABCDEFGHIJ0123456789
Country
<fctr>
lambda_guerrero
<dbl>
China-0.03446284

Making λ=0.2 looks ok.

The Guerrero method suggests an even stronger transformation. Let’s also try a log.

fit <- china |>
  model(
    ets = ETS(GDP),
    ets_damped = ETS(GDP ~ trend("Ad")),
    ets_bc = ETS(box_cox(GDP, 0.2)),
    ets_log = ETS(log(GDP))
  )

fit
ABCDEFGHIJ0123456789
Country
<fctr>
ets
<S3: lst_mdl>
ets_damped
<S3: lst_mdl>
ets_bc
<S3: lst_mdl>
ets_log
<S3: lst_mdl>
China<ETS(M,A,N)><ETS(M,Ad,N)><ETS(A,A,N)><ETS(A,A,N)>
augment(fit)
ABCDEFGHIJ0123456789
Country
<fctr>
.model
<chr>
Year
<dbl>
GDP
<dbl>
.fitted
<dbl>
.resid
<dbl>
.innov
<dbl>
Chinaets19605.971647e+104.900169e+101.071478e+102.186614e-01
Chinaets19615.005687e+106.634664e+10-1.628977e+10-2.455252e-01
Chinaets19624.720936e+105.160737e+10-4.398009e+09-8.522057e-02
Chinaets19635.070680e+104.738649e+103.320305e+097.006860e-02
Chinaets19645.970834e+105.191909e+107.789252e+091.500267e-01
Chinaets19657.043627e+106.335042e+107.085845e+091.118516e-01
Chinaets19667.672029e+107.628919e+104.310994e+085.650858e-03
Chinaets19677.288163e+108.270838e+10-9.826744e+09-1.188120e-01
Chinaets19687.084654e+107.580482e+10-4.958286e+09-6.540858e-02
Chinaets19697.970591e+107.222226e+107.483647e+091.036197e-01
fit |>
  forecast(h = "20 years") |>
  autoplot(china, level = NULL)

The transformations have a big effect, with small lambda values creating big increases in the forecasts. The damping has relatively a small effect.

Section 8.8 Problem#7

7.Find an ETS model for the Gas data from aus_production and forecast the next few years. Why is multiplicative seasonality necessary here? Experiment with making the trend damped. Does it improve the forecasts?

aus_production |> autoplot(Gas)

There is a huge increase in variance as the series increases in level which makes it necessary to use multiplicative seasonality.

fit <- aus_production |>
  model(
    hw = ETS(Gas ~ error("M") + trend("A") + season("M")),
    hwdamped = ETS(Gas ~ error("M") + trend("Ad") + season("M")),
  )

fit |> glance()
ABCDEFGHIJ0123456789
.model
<chr>
sigma2
<dbl>
log_lik
<dbl>
AIC
<dbl>
AICc
<dbl>
BIC
<dbl>
MSE
<dbl>
AMSE
<dbl>
hw0.003244894-831.46431680.9291681.7941711.38921.1150732.16296
hwdamped0.003285280-832.01391684.0281685.0911717.87321.0850032.04319

The non-damped model seems to be doing slightly better here, probably because the trend is very strong over most of the historical data.

fit |>
  select(hw) |>
  gg_tsresiduals()



fit |> tidy()
ABCDEFGHIJ0123456789
.model
<chr>
term
<chr>
estimate
<dbl>
hwalpha0.65285449
hwbeta0.14416745
hwgamma0.09784922
hwl[0]5.94559214
hwb[0]0.07062881
hws[0]0.93092356
hws[-1]1.17788274
hws[-2]1.07485100
hws[-3]0.81634270
hwdampedalpha0.64890445
# A tibble: 19 × 3

fit |>
  augment() |>
  filter(.model == "hw") |>
  features(.innov, ljung_box, lag = 24)
ABCDEFGHIJ0123456789
.model
<chr>
lb_stat
<dbl>
lb_pvalue
<dbl>
hw57.11360.0001613395

There is still some small correlations left in the residuals, showing the model has not fully captured the available information. There also appears to be some heteroskedasticity in the residuals with larger variance in the first half the series.

fit |>
  forecast(h = 36) |>
  filter(.model == "hw") |>
  autoplot(aus_production)

While the point forecasts look ok, the intervals are excessively wide.

Section 8.8 Problem#8

8.Recall your retail time series data (from Exercise 7 in Section 2.10).

a.Why is multiplicative seasonality necessary for this series?

set.seed(12345678)
myseries <- aus_retail |>
  filter(`Series ID` == sample(aus_retail$`Series ID`, 1))
myseries |> autoplot(Turnover)

The variation in the seasonal pattern increases as the level of the series rises. (This may not be true for every series, but is true for almost all of them.)

b. Apply Holt-Winters’ multiplicative method to the data. Experiment with making the trend damped.

fit <- myseries |>
  model(
    hw = ETS(Turnover ~ error("M") + trend("A") + season("M")),
    hwdamped = ETS(Turnover ~ error("M") + trend("Ad") + season("M"))
  )
fc <- fit |> forecast(h = 36)
fc |> autoplot(myseries)

c.Compare the RMSE of the one-step forecasts from the two methods. Which do you prefer?

accuracy(fit)
ABCDEFGHIJ0123456789
State
<chr>
Industry
<chr>
.model
<chr>
Northern TerritoryClothing, footwear and personal accessory retailinghw
Northern TerritoryClothing, footwear and personal accessory retailinghwdamped

The non-damped method is doing slightly better (on RMSE), but the damped method is doing better on most other scores. I’d be inclined to use the damped method here as the trend at the end of the series seems to be flattening.

d.Check that the residuals from the best method look like white noise.

fit |>
  select("hwdamped") |>
  gg_tsresiduals()

There are significant spikes at lags 8 and 18 in the ACF, but they are relatively small and probably of no consequence.

augment(fit) |>
  filter(.model == "hwdamped") |>
  features(.innov, ljung_box, lag = 36)
ABCDEFGHIJ0123456789
State
<chr>
Industry
<chr>
.model
<chr>
Northern TerritoryClothing, footwear and personal accessory retailinghwdamped

Overall, the autocorrelation in the residuals is not significant.

e.Now find the test set RMSE, while training the model to the end of 2010. Can you beat the seasonal naïve approach from Exercise 7 in Section 5.11?

myseries |>
  filter(year(Month) < 2011) |>
  model(
    snaive = SNAIVE(Turnover),
    hw = ETS(Turnover ~ error("M") + trend("A") + season("M")),
    hwdamped = ETS(Turnover ~ error("M") + trend("Ad") + season("M"))
  ) |>
  forecast(h = "7 years") |>
  accuracy(myseries)
ABCDEFGHIJ0123456789
.model
<chr>
State
<chr>
Industry
<chr>
hwNorthern TerritoryClothing, footwear and personal accessory retailing
hwdampedNorthern TerritoryClothing, footwear and personal accessory retailing
snaiveNorthern TerritoryClothing, footwear and personal accessory retailing

The SNAIVE model is doing much better than the HW model for this data set.

LS0tDQp0aXRsZTogIkV4cG9uZW50aWFsIFNtb290aGluZyBNb2RlbHMiDQpvdXRwdXQ6IA0KICBodG1sX25vdGVib29rOg0KICAgIHRvYzogdHJ1ZQ0KICAgIHRvY19mbG9hdDogdHJ1ZQ0KLS0tDQoNCiMjIFNlY3Rpb24gOC44IFByb2JsZW1bIzZdKGh0dHBzOi8vb3RleHRzLmNvbS9mcHAzL2V4cHNtb290aC1leGVyY2lzZXMuaHRtbCkNCg0KX182LiBGb3JlY2FzdCB0aGUgQ2hpbmVzZSBHRFAgZnJvbSB0aGUgYGdsb2JhbF9lY29ub215YCBkYXRhIHNldCB1c2luZyBhbiBFVFMgbW9kZWwuIEV4cGVyaW1lbnQgd2l0aCB0aGUgdmFyaW91cyBvcHRpb25zIGluIHRoZSBgRVRTKClgIGZ1bmN0aW9uIHRvIHNlZSBob3cgbXVjaCB0aGUgZm9yZWNhc3RzIGNoYW5nZSB3aXRoIGRhbXBlZCB0cmVuZCwgb3Igd2l0aCBhIEJveC1Db3ggdHJhbnNmb3JtYXRpb24uIFRyeSB0byBkZXZlbG9wIGFuIGludHVpdGlvbiBvZiB3aGF0IGVhY2ggaXMgZG9pbmcgdG8gdGhlIGZvcmVjYXN0cy5fXw0KDQpfX1tIaW50OiB1c2UgYSByZWxhdGl2ZWx5IGxhcmdlIHZhbHVlIG9mIGggd2hlbiBmb3JlY2FzdGluZywgc28geW91IGNhbiBjbGVhcmx5IHNlZSB0aGUgZGlmZmVyZW5jZXMgYmV0d2VlbiB0aGUgdmFyaW91cyBvcHRpb25zIHdoZW4gcGxvdHRpbmcgdGhlIGZvcmVjYXN0cy5dX18NCg0KYGBge3J9DQpjaGluYSA8LSBnbG9iYWxfZWNvbm9teSB8Pg0KICBmaWx0ZXIoQ291bnRyeSA9PSAiQ2hpbmEiKQ0KY2hpbmEgfD4gYXV0b3Bsb3QoR0RQKQ0KYGBgDQoNCkl0IGNsZWFybHkgbmVlZHMgYSByZWxhdGl2ZWx5IHN0cm9uZyB0cmFuc2Zvcm1hdGlvbiBkdWUgdG8gdGhlIGluY3JlYXNpbmcgdmFyaWFuY2UuDQoNCmBgYHtyfQ0KY2hpbmEgfD4gYXV0b3Bsb3QoYm94X2NveChHRFAsIDAuMikpDQpgYGANCg0KYGBge3J9DQpjaGluYSB8PiBmZWF0dXJlcyhHRFAsIGd1ZXJyZXJvKQ0KYGBgDQpNYWtpbmcgJFxsYW1iZGEgPSAwLjIkIGxvb2tzIG9rLg0KDQpUaGUgR3VlcnJlcm8gbWV0aG9kIHN1Z2dlc3RzIGFuIGV2ZW4gc3Ryb25nZXIgdHJhbnNmb3JtYXRpb24uIExldOKAmXMgYWxzbyB0cnkgYSBsb2cuDQoNCmBgYHtyfQ0KZml0IDwtIGNoaW5hIHw+DQogIG1vZGVsKA0KICAgIGV0cyA9IEVUUyhHRFApLA0KICAgIGV0c19kYW1wZWQgPSBFVFMoR0RQIH4gdHJlbmQoIkFkIikpLA0KICAgIGV0c19iYyA9IEVUUyhib3hfY294KEdEUCwgMC4yKSksDQogICAgZXRzX2xvZyA9IEVUUyhsb2coR0RQKSkNCiAgKQ0KDQpmaXQNCmBgYA0KDQpgYGB7cn0NCmF1Z21lbnQoZml0KQ0KYGBgDQoNCmBgYHtyfQ0KZml0IHw+DQogIGZvcmVjYXN0KGggPSAiMjAgeWVhcnMiKSB8Pg0KICBhdXRvcGxvdChjaGluYSwgbGV2ZWwgPSBOVUxMKQ0KYGBgDQoNClRoZSB0cmFuc2Zvcm1hdGlvbnMgaGF2ZSBhIGJpZyBlZmZlY3QsIHdpdGggc21hbGwgbGFtYmRhIHZhbHVlcyBjcmVhdGluZyBiaWcgaW5jcmVhc2VzIGluIHRoZSBmb3JlY2FzdHMuIFRoZSBkYW1waW5nIGhhcyByZWxhdGl2ZWx5IGEgc21hbGwgZWZmZWN0Lg0KDQoNCg0KIyMgU2VjdGlvbiA4LjggUHJvYmxlbVsjN10oaHR0cHM6Ly9vdGV4dHMuY29tL2ZwcDMvZXhwc21vb3RoLWV4ZXJjaXNlcy5odG1sKQ0KX183LkZpbmQgYW4gRVRTIG1vZGVsIGZvciB0aGUgR2FzIGRhdGEgZnJvbSBgYXVzX3Byb2R1Y3Rpb25gIGFuZCBmb3JlY2FzdCB0aGUgbmV4dCBmZXcgeWVhcnMuIFdoeSBpcyBtdWx0aXBsaWNhdGl2ZSBzZWFzb25hbGl0eSBuZWNlc3NhcnkgaGVyZT8gRXhwZXJpbWVudCB3aXRoIG1ha2luZyB0aGUgdHJlbmQgZGFtcGVkLiBEb2VzIGl0IGltcHJvdmUgdGhlIGZvcmVjYXN0cz9fXw0KDQpgYGB7cn0NCmF1c19wcm9kdWN0aW9uIHw+IGF1dG9wbG90KEdhcykNCmBgYA0KDQpUaGVyZSBpcyBhIGh1Z2UgaW5jcmVhc2UgaW4gdmFyaWFuY2UgYXMgdGhlIHNlcmllcyBpbmNyZWFzZXMgaW4gbGV2ZWwgd2hpY2ggbWFrZXMgaXQgbmVjZXNzYXJ5IHRvIHVzZSBtdWx0aXBsaWNhdGl2ZSBzZWFzb25hbGl0eS4NCg0KYGBge3J9DQpmaXQgPC0gYXVzX3Byb2R1Y3Rpb24gfD4NCiAgbW9kZWwoDQogICAgaHcgPSBFVFMoR2FzIH4gZXJyb3IoIk0iKSArIHRyZW5kKCJBIikgKyBzZWFzb24oIk0iKSksDQogICAgaHdkYW1wZWQgPSBFVFMoR2FzIH4gZXJyb3IoIk0iKSArIHRyZW5kKCJBZCIpICsgc2Vhc29uKCJNIikpLA0KICApDQoNCmZpdCB8PiBnbGFuY2UoKQ0KYGBgDQoNClRoZSBub24tZGFtcGVkIG1vZGVsIHNlZW1zIHRvIGJlIGRvaW5nIHNsaWdodGx5IGJldHRlciBoZXJlLCBwcm9iYWJseSBiZWNhdXNlIHRoZSB0cmVuZCBpcyB2ZXJ5IHN0cm9uZyBvdmVyIG1vc3Qgb2YgdGhlIGhpc3RvcmljYWwgZGF0YS4NCg0KYGBge3J9DQpmaXQgfD4NCiAgc2VsZWN0KGh3KSB8Pg0KICBnZ190c3Jlc2lkdWFscygpDQoNCg0KZml0IHw+IHRpZHkoKQ0KIyBBIHRpYmJsZTogMTkgw5cgMw0KDQpmaXQgfD4NCiAgYXVnbWVudCgpIHw+DQogIGZpbHRlcigubW9kZWwgPT0gImh3IikgfD4NCiAgZmVhdHVyZXMoLmlubm92LCBsanVuZ19ib3gsIGxhZyA9IDI0KQ0KYGBgDQoNClRoZXJlIGlzIHN0aWxsIHNvbWUgc21hbGwgY29ycmVsYXRpb25zIGxlZnQgaW4gdGhlIHJlc2lkdWFscywgc2hvd2luZyB0aGUgbW9kZWwgaGFzIG5vdCBmdWxseSBjYXB0dXJlZCB0aGUgYXZhaWxhYmxlIGluZm9ybWF0aW9uLg0KVGhlcmUgYWxzbyBhcHBlYXJzIHRvIGJlIHNvbWUgaGV0ZXJvc2tlZGFzdGljaXR5IGluIHRoZSByZXNpZHVhbHMgd2l0aCBsYXJnZXIgdmFyaWFuY2UgaW4gdGhlIGZpcnN0IGhhbGYgdGhlIHNlcmllcy4NCg0KYGBge3J9DQpmaXQgfD4NCiAgZm9yZWNhc3QoaCA9IDM2KSB8Pg0KICBmaWx0ZXIoLm1vZGVsID09ICJodyIpIHw+DQogIGF1dG9wbG90KGF1c19wcm9kdWN0aW9uKQ0KYGBgDQoNCldoaWxlIHRoZSBwb2ludCBmb3JlY2FzdHMgbG9vayBvaywgdGhlIGludGVydmFscyBhcmUgZXhjZXNzaXZlbHkgd2lkZS4NCg0KIyMgU2VjdGlvbiA4LjggUHJvYmxlbVsjOF0oaHR0cHM6Ly9vdGV4dHMuY29tL2ZwcDMvZXhwc21vb3RoLWV4ZXJjaXNlcy5odG1sKQ0KX184LlJlY2FsbCB5b3VyIHJldGFpbCB0aW1lIHNlcmllcyBkYXRhIChmcm9tIFtFeGVyY2lzZSA3XShodHRwczovL290ZXh0cy5jb20vZnBwMy9ncmFwaGljcy1leGVyY2lzZXMuaHRtbCNncmFwaGljcy1leGVyY2lzZXMpIGluIFNlY3Rpb24gMi4xMCkuX18NCg0KX19hLldoeSBpcyBtdWx0aXBsaWNhdGl2ZSBzZWFzb25hbGl0eSBuZWNlc3NhcnkgZm9yIHRoaXMgc2VyaWVzP19fDQoNCmBgYHtyfQ0Kc2V0LnNlZWQoMTIzNDU2NzgpDQpteXNlcmllcyA8LSBhdXNfcmV0YWlsIHw+DQogIGZpbHRlcihgU2VyaWVzIElEYCA9PSBzYW1wbGUoYXVzX3JldGFpbCRgU2VyaWVzIElEYCwgMSkpDQpteXNlcmllcyB8PiBhdXRvcGxvdChUdXJub3ZlcikNCmBgYA0KDQpUaGUgdmFyaWF0aW9uIGluIHRoZSBzZWFzb25hbCBwYXR0ZXJuIGluY3JlYXNlcyBhcyB0aGUgbGV2ZWwgb2YgdGhlIHNlcmllcyByaXNlcy4gKFRoaXMgbWF5IG5vdCBiZSB0cnVlIGZvciBldmVyeSBzZXJpZXMsIGJ1dCBpcyB0cnVlIGZvciBhbG1vc3QgYWxsIG9mIHRoZW0uKQ0KDQpfX2IuIEFwcGx5IEhvbHQtV2ludGVyc+KAmSBtdWx0aXBsaWNhdGl2ZSBtZXRob2QgdG8gdGhlIGRhdGEuIEV4cGVyaW1lbnQgd2l0aCBtYWtpbmcgdGhlIHRyZW5kIGRhbXBlZC5fXw0KDQpgYGB7cn0NCmZpdCA8LSBteXNlcmllcyB8Pg0KICBtb2RlbCgNCiAgICBodyA9IEVUUyhUdXJub3ZlciB+IGVycm9yKCJNIikgKyB0cmVuZCgiQSIpICsgc2Vhc29uKCJNIikpLA0KICAgIGh3ZGFtcGVkID0gRVRTKFR1cm5vdmVyIH4gZXJyb3IoIk0iKSArIHRyZW5kKCJBZCIpICsgc2Vhc29uKCJNIikpDQogICkNCmZjIDwtIGZpdCB8PiBmb3JlY2FzdChoID0gMzYpDQpmYyB8PiBhdXRvcGxvdChteXNlcmllcykNCmBgYA0KX19jLkNvbXBhcmUgdGhlIFJNU0Ugb2YgdGhlIG9uZS1zdGVwIGZvcmVjYXN0cyBmcm9tIHRoZSB0d28gbWV0aG9kcy4gV2hpY2ggZG8geW91IHByZWZlcj9fXw0KDQpgYGB7cn0NCmFjY3VyYWN5KGZpdCkNCmBgYA0KDQpUaGUgbm9uLWRhbXBlZCBtZXRob2QgaXMgZG9pbmcgc2xpZ2h0bHkgYmV0dGVyIChvbiBSTVNFKSwgYnV0IHRoZSBkYW1wZWQgbWV0aG9kIGlzIGRvaW5nIGJldHRlciBvbiBtb3N0IG90aGVyIHNjb3Jlcy4gSeKAmWQgYmUgaW5jbGluZWQgdG8gdXNlIHRoZSBkYW1wZWQgbWV0aG9kIGhlcmUgYXMgdGhlIHRyZW5kIGF0IHRoZSBlbmQgb2YgdGhlIHNlcmllcyBzZWVtcyB0byBiZSBmbGF0dGVuaW5nLg0KDQpfX2QuQ2hlY2sgdGhhdCB0aGUgcmVzaWR1YWxzIGZyb20gdGhlIGJlc3QgbWV0aG9kIGxvb2sgbGlrZSB3aGl0ZSBub2lzZS5fXw0KDQpgYGB7cn0NCmZpdCB8Pg0KICBzZWxlY3QoImh3ZGFtcGVkIikgfD4NCiAgZ2dfdHNyZXNpZHVhbHMoKQ0KYGBgDQoNClRoZXJlIGFyZSBzaWduaWZpY2FudCBzcGlrZXMgYXQgbGFncyA4IGFuZCAxOCBpbiB0aGUgQUNGLCBidXQgdGhleSBhcmUgcmVsYXRpdmVseSBzbWFsbCBhbmQgcHJvYmFibHkgb2Ygbm8gY29uc2VxdWVuY2UuDQoNCmBgYHtyfQ0KYXVnbWVudChmaXQpIHw+DQogIGZpbHRlcigubW9kZWwgPT0gImh3ZGFtcGVkIikgfD4NCiAgZmVhdHVyZXMoLmlubm92LCBsanVuZ19ib3gsIGxhZyA9IDM2KQ0KYGBgDQpPdmVyYWxsLCB0aGUgYXV0b2NvcnJlbGF0aW9uIGluIHRoZSByZXNpZHVhbHMgaXMgbm90IHNpZ25pZmljYW50Lg0KDQoNCl9fZS5Ob3cgZmluZCB0aGUgdGVzdCBzZXQgUk1TRSwgd2hpbGUgdHJhaW5pbmcgdGhlIG1vZGVsIHRvIHRoZSBlbmQgb2YgMjAxMC4gQ2FuIHlvdSBiZWF0IHRoZSBzZWFzb25hbCBuYcOvdmUgYXBwcm9hY2ggZnJvbSBbRXhlcmNpc2UgN10oaHR0cHM6Ly9vdGV4dHMuY29tL2ZwcDMvdG9vbGJveC1leGVyY2lzZXMuaHRtbCN0b29sYm94LWV4ZXJjaXNlcykgaW4gU2VjdGlvbiA1LjExP19fDQoNCmBgYHtyfQ0KbXlzZXJpZXMgfD4NCiAgZmlsdGVyKHllYXIoTW9udGgpIDwgMjAxMSkgfD4NCiAgbW9kZWwoDQogICAgc25haXZlID0gU05BSVZFKFR1cm5vdmVyKSwNCiAgICBodyA9IEVUUyhUdXJub3ZlciB+IGVycm9yKCJNIikgKyB0cmVuZCgiQSIpICsgc2Vhc29uKCJNIikpLA0KICAgIGh3ZGFtcGVkID0gRVRTKFR1cm5vdmVyIH4gZXJyb3IoIk0iKSArIHRyZW5kKCJBZCIpICsgc2Vhc29uKCJNIikpDQogICkgfD4NCiAgZm9yZWNhc3QoaCA9ICI3IHllYXJzIikgfD4NCiAgYWNjdXJhY3kobXlzZXJpZXMpDQpgYGANCg0KVGhlIFNOQUlWRSBtb2RlbCBpcyBkb2luZyBtdWNoIGJldHRlciB0aGFuIHRoZSBIVyBtb2RlbCBmb3IgdGhpcyBkYXRhIHNldC4=