Homework 2

Exercise 7.8

1

A)

??pigs

library(fma)

## Loading required package: forecast

## Warning: package 'forecast' was built under R version 3.4.4

pigs

##         Jan    Feb    Mar    Apr    May    Jun    Jul    Aug    Sep    Oct
## 1980  76378  71947  33873  96428 105084  95741 110647 100331  94133 103055
## 1981  76889  81291  91643  96228 102736 100264 103491  97027  95240  91680
## 1982  76892  85773  95210  93771  98202  97906 100306  94089 102680  77919
## 1983  81225  88357 106175  91922 104114 109959  97880 105386  96479  97580
## 1984  90974  98981 107188  94177 115097 113696 114532 120110  93607 110925
## 1985 103069 103351 111331 106161 111590  99447 101987  85333  86970 100561
## 1986  82719  79498  74846  73819  77029  78446  86978  75878  69571  75722
## 1987  63292  59380  78332  72381  55971  69750  85472  70133  79125  85805
## 1988  69069  79556  88174  66698  72258  73445  76131  86082  75443  73969
## 1989  66269  73776  80034  70694  81823  75640  75540  82229  75345  77034
## 1990  75982  78074  77588  84100  97966  89051  93503  84747  74531  91900
## 1991  81022  78265  77271  85043  95418  79568 103283  95770  91297 101244
## 1992  93866  95171 100183 103926 102643 108387  97077  90901  90336  88732
## 1993  73292  78943  94399  92937  90130  91055 106062 103560 104075 101783
## 1994  82413  83534 109011  96499 102430 103002  91815  99067 110067 101599
## 1995  88905  89936 106723  84307 114896 106749  87892 100506              
##         Nov    Dec
## 1980  90595 101457
## 1981 101259 109564
## 1982  93561 117062
## 1983 109490 110191
## 1984 103312 120184
## 1985  89543  89265
## 1986  64182  77357
## 1987  81778  86852
## 1988  78139  78646
## 1989  78589  79769
## 1990  81635  89797
## 1991 114525 101139
## 1992  83759  99267
## 1993  93791 102313
## 1994  97646 104930
## 1995

ses_pigs <-ses(pigs, h = 4)
summary(ses_pigs)

## 
## Forecast method: Simple exponential smoothing
## 
## Model Information:
## Simple exponential smoothing 
## 
## Call:
##  ses(y = pigs, h = 4) 
## 
##   Smoothing parameters:
##     alpha = 0.2971 
## 
##   Initial states:
##     l = 77260.0561 
## 
##   sigma:  10308.58
## 
##      AIC     AICc      BIC 
## 4462.955 4463.086 4472.665 
## 
## Error measures:
##                    ME    RMSE      MAE       MPE     MAPE      MASE
## Training set 385.8721 10253.6 7961.383 -0.922652 9.274016 0.7966249
##                    ACF1
## Training set 0.01282239
## 
## Forecasts:
##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Sep 1995       98816.41 85605.43 112027.4 78611.97 119020.8
## Oct 1995       98816.41 85034.52 112598.3 77738.83 119894.0
## Nov 1995       98816.41 84486.34 113146.5 76900.46 120732.4
## Dec 1995       98816.41 83958.37 113674.4 76092.99 121539.8

B)

ses_pigs$upper[1, "95%"]

##      95% 
## 119020.8

ses_pigs$lower[1, "95%"]

##      95% 
## 78611.97

Prediction interval for first forecast

s <- sd(ses_pigs$residuals)
ses_pigs$mean[1] + 1.96*s

## [1] 118952.8

ses_pigs$mean[1] - 1.96*s

## [1] 78679.97

Larger confidence interval than the one produced by R

autoplot(ses_pigs) +
  autolayer(ses_pigs$fitted)

Plotted forecasts

12

A)

  fets <- function(y, h) {
  forecast(ets(y), h = h)
  }

library(fpp2)

## Warning: package 'fpp2' was built under R version 3.4.4

## Loading required package: ggplot2

## Loading required package: expsmooth

naive.qcement <- naive(qcement)

fets.qcement <- fets(qcement, 4)

snaive.qcement <- snaive(qcement, h = 4)

mean(qcement, na.rm = TRUE)

## [1] 1.482369

B)

mean(tsCV(qcement, fets.qcement, h = 4)^2, na.rm = TRUE)

## [1] NaN

mean(tsCV(qcement, snaive(qcement), h = 4)^2, na.rm = TRUE)

## [1] NaN

NaNs produced? Cant figure it out.

13

Ausbeer

Training and Test Sets

ausbeer_train <- subset(
  ausbeer, end = length(ausbeer) - 12
  )
ausbeer_test <- subset(
  ausbeer, start = length(ausbeer) - 11
  )

Models

ets_ausbeer_train <- forecast(
  ets(ausbeer_train), h = 12
)
snaive_ausbeer_train <- snaive(ausbeer_train,  h = 12)
stlf_ausbeer_train <- stlf(
  ausbeer_train, 
  h = 12,
  s.window = 5,
  robust = TRUE,
  lambda = BoxCox.lambda(ausbeer_train))

Model Selection

accuracy(ets_ausbeer_train, ausbeer_test)

##                      ME      RMSE       MAE          MPE     MAPE
## Training set -0.3466095 15.781973 11.979955 -0.064073435 2.864311
## Test set      0.1272571  9.620828  8.919347  0.009981598 2.132836
##                   MASE       ACF1 Theil's U
## Training set 0.7567076 -0.1942276        NA
## Test set     0.5633859  0.3763918 0.1783972

accuracy(snaive_ausbeer_train, ausbeer_test)

##                     ME     RMSE      MAE        MPE     MAPE      MASE
## Training set  3.336634 19.67005 15.83168  0.9044009 3.771370 1.0000000
## Test set     -2.916667 10.80509  9.75000 -0.6505927 2.338012 0.6158537
##                      ACF1 Theil's U
## Training set 0.0009690785        NA
## Test set     0.3254581810 0.1981463

accuracy(stlf_ausbeer_train, ausbeer_test)

##                      ME     RMSE      MAE        MPE     MAPE      MASE
## Training set  0.5704833 13.61609 9.816091  0.1279204 2.334194 0.6200283
## Test set     -4.0433841 10.11436 8.429747 -0.9670308 2.051915 0.5324605
##                    ACF1 Theil's U
## Training set -0.1131945        NA
## Test set      0.2865992 0.1689574

Best Model is STL + ETS

autoplot(stlf_ausbeer_train) +
  autolayer(ausbeer_test)

Forecasting Function (GitHub)

forecast.models <- function(y, h){
  m <- frequency(y)
  
  y_train <- subset(
    y, end = length(y) - m*h
    )
  
  y_test <- subset(
    y, start = length(y) - m*h + 1
    )
  ets_y_train <- forecast(
    ets(y_train), h = m*h
  )
  
  snaive_y_train <- snaive(y_train,  h = m*h)
  
  stlf_y_train <- stlf(
    y_train, 
    h = m*h,
    s.window = m + 1,
    robust = TRUE
    )
  
  stlf_y_train_with_BoxCox <- stlf(
    y_train, 
    h = m*h,
    s.window = m + 1,
    robust = TRUE,
    lambda = BoxCox.lambda(y_train))
  models <- list(ets_y_train, 
                 snaive_y_train, 
                 stlf_y_train,
                 stlf_y_train_with_BoxCox)
  
  names(models) <- c("ets", "snaive", "stlf", "stlf_with_BoxCox")
  acc_ets <- accuracy(ets_y_train, y_test)
  acc_snaive <- accuracy(snaive_y_train, y_test)
  acc_stlf <- accuracy(stlf_y_train, y_test)
  acc_stlf_with_BoxCox <- accuracy(stlf_y_train_with_BoxCox, y_test)
  accuracies <- list(acc_ets, 
                     acc_snaive, 
                     acc_stlf,
                     acc_stlf_with_BoxCox)
  
  names(accuracies) <- c("acc_ets", "acc_snaive", "acc_stlf", "acc_stlf_with_BoxCox")
  output <- list(y_train, y_test, m, models, accuracies)
  names(output) <- c("train", "test", "m", "models", "accuracies")
  
  return(output)
}

Bricksq

fc_bricksq <- forecast.models(bricksq, 3)
fc_bricksq$accuracies

## $acc_ets
##                     ME     RMSE      MAE       MPE     MAPE      MASE
## Training set  1.228708 21.96146 15.81586 0.2606171 3.912882 0.4299638
## Test set     37.682916 42.36696 37.68292 8.4157549 8.415755 1.0244329
##                   ACF1 Theil's U
## Training set 0.2038074        NA
## Test set     0.6190546  1.116608
## 
## $acc_snaive
##                     ME     RMSE      MAE      MPE     MAPE      MASE
## Training set  6.438849 50.25482 36.78417 1.430237 8.999949 1.0000000
## Test set     15.333333 37.15284 33.66667 3.231487 7.549326 0.9152487
##                     ACF1 Theil's U
## Training set  0.81169827        NA
## Test set     -0.09314867 0.9538702
## 
## $acc_stlf
##                     ME     RMSE      MAE        MPE      MAPE      MASE
## Training set  1.318736 21.49082 14.41766  0.2855898  3.501287 0.3919528
## Test set     44.266931 50.20039 45.89301 10.0933739 10.487099 1.2476294
##                   ACF1 Theil's U
## Training set 0.1502759        NA
## Test set     0.1053316   1.31732
## 
## $acc_stlf_with_BoxCox
##                    ME     RMSE      MAE       MPE     MAPE      MASE
## Training set  1.42724 21.40466 14.30375 0.3185472 3.500221 0.3888560
## Test set     34.69481 40.13618 36.34649 7.7322091 8.132133 0.9881014
##                   ACF1 Theil's U
## Training set 0.1625098        NA
## Test set     0.4746242   1.04561

autoplot(fc_bricksq$models$snaive) +
  autolayer(fc_bricksq$test)

Seasonal Naive is best

Dole

fc_dole <- forecast.models(dole, 3)
fc_dole$accuracies

## $acc_ets
##                        ME      RMSE        MAE        MPE     MAPE
## Training set     98.00253  16130.58   9427.497  0.5518769  6.20867
## Test set     221647.42595 279668.65 221647.426 32.5447637 32.54476
##                   MASE      ACF1 Theil's U
## Training set 0.2995409 0.5139536        NA
## Test set     7.0424271 0.9394267  11.29943
## 
## $acc_snaive
##                     ME      RMSE       MAE       MPE     MAPE     MASE
## Training set  12606.58  56384.06  31473.16  3.350334 27.73651 1.000000
## Test set     160370.33 240830.92 186201.00 20.496197 27.38678 5.916184
##                   ACF1 Theil's U
## Training set 0.9781058        NA
## Test set     0.9208325  9.479785
## 
## $acc_stlf
##                       ME       RMSE       MAE        MPE      MAPE
## Training set    -96.4811   7801.958   4530.06 -0.2719573  5.116149
## Test set     328005.9874 407547.190 328005.99 48.6435815 48.643581
##                    MASE       ACF1 Theil's U
## Training set  0.1439341 0.08037493        NA
## Test set     10.4217690 0.93373958  16.57033
## 
## $acc_stlf_with_BoxCox
##                       ME       RMSE        MAE        MPE     MAPE
## Training set    145.1468   6688.083   3659.404  0.1464869  3.82385
## Test set     205385.6547 268135.992 207317.238 29.3540384 29.85051
##                   MASE        ACF1 Theil's U
## Training set 0.1162706 -0.09446256        NA
## Test set     6.5871126  0.93778748  10.68587

autoplot(fc_dole$models$snaive) +
  autolayer(fc_dole$test)

Seasonal Naive is Best

A10

fc_a10 <- forecast.models(a10, 3)
fc_a10$accuracies

## $acc_ets
##                      ME     RMSE       MAE       MPE     MAPE      MASE
## Training set 0.04378049 0.488989 0.3542285 0.2011915 3.993267 0.3611299
## Test set     1.62559331 2.562046 2.0101563 7.0011207 9.410027 2.0493197
##                     ACF1 Theil's U
## Training set -0.05238173        NA
## Test set      0.23062972 0.6929693
## 
## $acc_snaive
##                     ME     RMSE       MAE      MPE     MAPE     MASE
## Training set 0.9577207 1.177528 0.9808895 10.86283 11.15767 1.000000
## Test set     4.3181585 5.180738 4.3306974 19.80542 19.89352 4.415071
##                   ACF1 Theil's U
## Training set 0.3779746        NA
## Test set     0.6384822  1.383765
## 
## $acc_stlf
##                      ME      RMSE       MAE       MPE      MAPE      MASE
## Training set 0.06704027 0.6809938 0.4297557 0.3200743  4.732989 0.4381286
## Test set     1.83324369 3.1118106 2.4570607 7.0963340 11.223584 2.5049311
##                     ACF1 Theil's U
## Training set -0.01970509        NA
## Test set      0.23918322 0.8089697
## 
## $acc_stlf_with_BoxCox
##                      ME      RMSE       MAE         MPE     MAPE      MASE
## Training set 0.01686093 0.4244591 0.3098511 0.003515845 3.594686 0.3158879
## Test set     1.25869582 2.2242431 1.8428354 5.032385108 8.701050 1.8787390
##                    ACF1 Theil's U
## Training set -0.1780966        NA
## Test set      0.1502463 0.5964321

autoplot(fc_a10$models$stlf_with_BoxCox) +
  autolayer(fc_a10$test)

STL + ETS with Box Cox is best

H20

fc_h02 <- forecast.models(h02, 3)
fc_h02$accuracies

## $acc_ets
##                       ME       RMSE        MAE          MPE     MAPE
## Training set 0.001532209 0.04653434 0.03463451 0.0008075938 4.575471
## Test set     0.023667588 0.07667744 0.06442193 1.7152319315 7.030603
##                   MASE        ACF1 Theil's U
## Training set 0.5850192  0.07982687        NA
## Test set     1.0881653 -0.12074883  0.450176
## 
## $acc_snaive
##                       ME       RMSE        MAE       MPE     MAPE    MASE
## Training set  0.03880677 0.07122149 0.05920234  5.220203 8.156509 1.00000
## Test set     -0.01479486 0.08551752 0.07161055 -1.308298 7.884556 1.20959
##                    ACF1 Theil's U
## Training set 0.40465289        NA
## Test set     0.02264221 0.5009677
## 
## $acc_stlf
##                       ME       RMSE        MAE        MPE     MAPE
## Training set 0.001987465 0.04609626 0.03396307 -0.1760405 4.751272
## Test set     0.046520958 0.09174622 0.07544541  3.5302072 8.001263
##                   MASE       ACF1 Theil's U
## Training set 0.5736778 0.01953927        NA
## Test set     1.2743653 0.05661406 0.5085785
## 
## $acc_stlf_with_BoxCox
##                       ME       RMSE        MAE        MPE     MAPE
## Training set 0.003628017 0.04230875 0.03068561 0.02682146 4.159581
## Test set     0.031706735 0.07574975 0.06381807 2.51742657 6.903027
##                   MASE        ACF1 Theil's U
## Training set 0.5183175 -0.09600393        NA
## Test set     1.0779653 -0.25297574 0.4385025

autoplot(fc_h02$models$stlf_with_BoxCox) +
  autolayer(fc_h02$test)

STL + ETS is best

USMELEC

fc_usmelec <- forecast.models(usmelec, 3)
fc_usmelec$accuracies

## $acc_ets
##                        ME      RMSE       MAE        MPE     MAPE
## Training set  -0.07834851  7.447167  5.601963 -0.1168709 2.163495
## Test set     -10.20080652 15.291359 13.123441 -3.1908161 3.926338
##                   MASE      ACF1 Theil's U
## Training set 0.6225637 0.2719249        NA
## Test set     1.4584491 0.4679496 0.4859495
## 
## $acc_snaive
##                    ME     RMSE       MAE      MPE     MAPE     MASE
## Training set 4.921564 11.58553  8.998217 2.000667 3.511648 1.000000
## Test set     5.475972 17.20037 13.494750 1.391437 3.767514 1.499714
##                   ACF1 Theil's U
## Training set 0.4841860        NA
## Test set     0.2692544 0.4765145
## 
## $acc_stlf
##                        ME      RMSE       MAE         MPE     MAPE
## Training set  -0.07064178  6.659696  4.907696 -0.07729202 1.910234
## Test set     -11.41131201 17.776488 15.918880 -3.66455560 4.779167
##                   MASE      ACF1 Theil's U
## Training set 0.5454076 0.1126383        NA
## Test set     1.7691150 0.5099758 0.5639709
## 
## $acc_stlf_with_BoxCox
##                       ME      RMSE       MAE         MPE     MAPE
## Training set  -0.1351692  6.474779  4.761143 -0.05026415 1.848092
## Test set     -19.9969503 24.406877 20.951753 -6.06977735 6.315723
##                   MASE       ACF1 Theil's U
## Training set 0.5291207 0.08255205        NA
## Test set     2.3284338 0.60877348 0.7777035

autoplot(fc_usmelec$models$ets) +
  autolayer(fc_usmelec$test)

ETS is best

Exercise 8.11

7

A)

autoplot(wmurders)

autoplot(diff(wmurders))

ndiffs(wmurders)

## [1] 2

2 differencing needed

autoplot(diff(wmurders, differences = 2))

Stationary

diff(wmurders, differences = 2) %>% ggtsdisplay()

Lag 1 and 2 show spikes..

Use ARIMA(0, 2, 2)

B)

No costant should be used becasue of twice differencing

C)

???

D)

wmurders_arima <- Arima(wmurders, 
                              order = c(0, 2, 2))
checkresiduals(wmurders_arima)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,2,2)
## Q* = 11.764, df = 8, p-value = 0.1621
## 
## Model df: 2.   Total lags used: 10

Residuals look good.

E)

fc_wmurders_arima <- forecast(
  wmurders_arima, h = 3
)

fc_wmurders_arima$mean

## Time Series:
## Start = 2005 
## End = 2007 
## Frequency = 1 
## [1] 2.480525 2.374890 2.269256

F)

autoplot(fc_wmurders_arima)

####G)

fc_wmurders_autoarima <- forecast(
  auto.arima(wmurders), h = 3
)
fc_wmurders_autoarima$mean

## Time Series:
## Start = 2005 
## End = 2007 
## Frequency = 1 
## [1] 2.470660 2.363106 2.252833

Slightly different

accuracy(fc_wmurders_arima)

##                      ME      RMSE       MAE        MPE     MAPE      MASE
## Training set -0.0113461 0.2088162 0.1525773 -0.2403396 4.331729 0.9382785
##                     ACF1
## Training set -0.05094066

accuracy(fc_wmurders_autoarima)

##                       ME      RMSE       MAE        MPE     MAPE      MASE
## Training set -0.01065956 0.2072523 0.1528734 -0.2149476 4.335214 0.9400996
##                    ACF1
## Training set 0.02176343

Very close, second model slightly better

8

A)

autoplot(austa)

fc_austa_autoarima <- forecast(
  auto.arima(austa), h = 10
)
fc_austa_autoarima$model

## Series: austa 
## ARIMA(0,1,1) with drift 
## 
## Coefficients:
##          ma1   drift
##       0.3006  0.1735
## s.e.  0.1647  0.0390
## 
## sigma^2 estimated as 0.03376:  log likelihood=10.62
## AIC=-15.24   AICc=-14.46   BIC=-10.57

checkresiduals(fc_austa_autoarima)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,1) with drift
## Q* = 2.297, df = 5.2, p-value = 0.8266
## 
## Model df: 2.   Total lags used: 7.2

autoplot(fc_austa_autoarima)

ARIMA (0,1,1) was chosen

The residuals look like white noise

B)

fc_austa_arima.0.1.1 <- forecast(
  Arima(austa, order = c(0, 1, 1)), h = 10
)
autoplot(fc_austa_arima.0.1.1)

fc_austa_arima.0.1.0 <- forecast(
  Arima(austa, order = c(0, 1, 0)), h = 10
)
autoplot(fc_austa_arima.0.1.0)

The increasing trend is not reflected in the forecast

C)

fc_austa_arima.2.1.3.drift <- forecast(
  Arima(austa, order = c(2, 1, 3), include.drift = TRUE),
  h = 10
)
autoplot(fc_austa_arima.2.1.3.drift)

The forecast now refelcts the increasing trend

drift_austa <- fc_austa_arima.2.1.3.drift$model$coef[6]
fc_austa_arima.2.1.3.nodrift <- fc_austa_arima.2.1.3.drift$mean - drift_austa*seq_len(10)
autoplot(fc_austa_arima.2.1.3.drift) +
  autolayer(fc_austa_arima.2.1.3.nodrift)

Forecast looks worse without drift constant

D)

fc_austa_arima.0.0.1.const <- forecast(
  Arima(
    austa, order = c(0, 0, 1), include.constant = TRUE
    ),
  h = 10
)
autoplot(fc_austa_arima.0.0.1.const)

Forecast immediately drops to the overall mean of the data’s history

E)

fc_austa_arima.0.2.1 <- forecast(
  Arima(austa, order = c(0, 2, 1)),
  h = 10
)
autoplot(fc_austa_arima.0.2.1)

Interval grows as forecast goes further out, shows increasing trend.

9

A)

autoplot(usgdp)

autoplot(BoxCox(usgdp, BoxCox.lambda(usgdp)))

Box Cox is good

lambda_usgdp <- BoxCox.lambda(usgdp)

B)

usgdp_autoarima <- auto.arima(usgdp, 
                              lambda = lambda_usgdp)
autoplot(usgdp, series = "Data") +
  autolayer(usgdp_autoarima$fitted, series = "Fitted")

Model fits well

usgdp_autoarima

## Series: usgdp 
## ARIMA(2,1,0) with drift 
## Box Cox transformation: lambda= 0.366352 
## 
## Coefficients:
##          ar1     ar2   drift
##       0.2795  0.1208  0.1829
## s.e.  0.0647  0.0648  0.0202
## 
## sigma^2 estimated as 0.03518:  log likelihood=61.56
## AIC=-115.11   AICc=-114.94   BIC=-101.26

C)

ndiffs(BoxCox(usgdp, lambda_usgdp))

## [1] 1

1 differencing will make it stationary

usgdp_arima.1.1.0 <- Arima(
  usgdp, lambda = lambda_usgdp, order = c(1, 1, 0)
)
usgdp_arima.1.1.0

## Series: usgdp 
## ARIMA(1,1,0) 
## Box Cox transformation: lambda= 0.366352 
## 
## Coefficients:
##          ar1
##       0.6326
## s.e.  0.0504
## 
## sigma^2 estimated as 0.04384:  log likelihood=34.39
## AIC=-64.78   AICc=-64.73   BIC=-57.85

autoplot(usgdp, series = "Data") +
  autolayer(usgdp_arima.1.1.0$fitted, series = "Fitted")

Looks good too

usgdp_arima.1.1.0.drift <- Arima(
  usgdp, lambda = lambda_usgdp, order = c(1, 1, 0),
  include.drift = TRUE
)
usgdp_arima.1.1.0.drift

## Series: usgdp 
## ARIMA(1,1,0) with drift 
## Box Cox transformation: lambda= 0.366352 
## 
## Coefficients:
##          ar1   drift
##       0.3180  0.1831
## s.e.  0.0619  0.0179
## 
## sigma^2 estimated as 0.03555:  log likelihood=59.83
## AIC=-113.66   AICc=-113.56   BIC=-103.27

autoplot(usgdp, series = "Data") +
  autolayer(usgdp_arima.1.1.0.drift$fitted, series = "Fitted")

Both fit really well

D)

accuracy(usgdp_autoarima)

##                    ME    RMSE      MAE         MPE      MAPE      MASE
## Training set 1.195275 39.2224 29.29521 -0.01363259 0.6863491 0.1655687
##                     ACF1
## Training set -0.03824844

accuracy(usgdp_arima.1.1.0)

##                    ME     RMSE      MAE       MPE      MAPE     MASE
## Training set 15.45449 45.49569 35.08393 0.3101283 0.7815664 0.198285
##                    ACF1
## Training set -0.3381619

accuracy(usgdp_arima.1.1.0.drift)

##                    ME     RMSE     MAE         MPE      MAPE      MASE
## Training set 1.315796 39.90012 29.5802 -0.01678591 0.6834509 0.1671794
##                     ACF1
## Training set -0.08544569

Choose between first and third model

checkresiduals(usgdp_autoarima)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(2,1,0) with drift
## Q* = 6.5772, df = 5, p-value = 0.254
## 
## Model df: 3.   Total lags used: 8

checkresiduals(usgdp_arima.1.1.0.drift)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,1,0) with drift
## Q* = 10.274, df = 6, p-value = 0.1136
## 
## Model df: 2.   Total lags used: 8

Third model looks to be the best

E)

fc_usgdp_autoarima <- forecast(
  usgdp_autoarima
)
autoplot(fc_usgdp_autoarima)

Forecasts look good

F)

fc_usgdp_ets <- forecast(
  ets(usgdp)
)
autoplot(fc_usgdp_ets)

These look slightly better than those of the ARIMA model

Exercise 9.7

4

A)

autoplot(gasoline)

3 notable parts

t <- time(gasoline)
t.knot1 <- 2007.5
t.knot2 <- 2013
t.pw1 <- ts(pmax(0, t - t.knot1), start = t[1],
            frequency = 365.25/7)
t.pw2 <- ts(pmax(0, t - t.knot2), start = t[1],
            frequency = 365.25/7)

Default AIC infinity and pairs 0

AICc <- Inf
K_min.Aicc <- 0

for(num in c(1:26)){
  gasoline_tslm <- tslm(
    gasoline ~ trend + t.pw1 + t.pw2 + fourier(
      gasoline, K = num
    )
  )
  AICc_value <- CV(gasoline_tslm)["AICc"]
  
  if(AICc > AICc_value){
    AICc <- AICc_value
  }else{
    K_min.Aicc <- num
    break 
  }
}
K_min.Aicc

## [1] 11

11 pairs chosen

autoplot(gasoline) +
  geom_line(color = "gray") +
  autolayer(gasoline_tslm$fitted.values)

Fitted values follow trend of data and are very similar

B)

gasoline_autoarima <- Arima( gasoline, xreg = cbind(t=t, t.pw1=t.pw1, t.pw2=t.pw2, Fourier = fourier(gasoline, K = 11)), order = c(4, 0, 2), seasonal = c(1, 0, 0) ) gasoline_autoarima

Would not run code

C)

checkresiduals(gasoline_autoarima)

???

10.8

4

str(visnights)

##  Time-Series [1:76, 1:20] from 1998 to 2017: 9.05 6.96 6.87 7.15 7.96 ...
##  - attr(*, "dimnames")=List of 2
##   ..$ : NULL
##   ..$ : chr [1:20] "NSWMetro" "NSWNthCo" "NSWSthCo" "NSWSthIn" ...

require("hts")

## Loading required package: hts

## Warning: package 'hts' was built under R version 3.4.4

library(hts)
tourism.hts <- hts(visnights, characters = c(3, 5))
str(smatrix(tourism.hts))

##  num [1:27, 1:20] 1 1 0 0 0 0 0 1 0 0 ...
##  - attr(*, "dimnames")=List of 2
##   ..$ : NULL
##   ..$ : NULL

fc_visnights_arima_opt <- forecast(
  tourism.hts, h = 8,
  method = "comb", weights = "mint", covariance = "sam",
  fmethod = "arima"
)

fc_visnights_arima <- forecast(
  tourism.hts, h = 8, method = "bu", fmethod = "arima"
)
plot(fc_visnights_arima_opt, levels = 0, col = "red")
par(new = TRUE, xpd = TRUE)
plot(fc_visnights_arima, levels = 0, col = "blue")
title(main = "Total number of visitors")
legend("bottomright", legend = c("Optimal", "Bottom-up"), title = "Coherent forecast", col = c("red", "blue"), lty = c(1, 1), bty = "n", cex = 0.5)

This method yielded higher forecasts than the bottom up method but smaller than the ARIMA model

plot(fc_visnights_arima_opt, levels = 1, color_lab = TRUE)

plot(fc_visnights_arima, levels = 1, color_lab = TRUE)

plot(fc_visnights_arima_opt, levels = 2, color_lab = TRUE)

plot(fc_visnights_arima, levels = 2, color_lab = TRUE)

The optimally reconciled method yielded higher forecast values especially for peaks than the bottom-up method.

Some forecasts also became bigger when optimally reconciled method was used and for some the forecasts became smaller.

11.5

3

Cant load data