Hyndman Ch 7 and 8 Exercises

Demetri Chokshi-Fox"

Chapter 7

library(forecast)

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

library(fpp2)

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

## Loading required package: ggplot2

## Warning: package 'ggplot2' was built under R version 3.5.3

## Loading required package: fma

## Warning: package 'fma' was built under R version 3.5.3

## Loading required package: expsmooth

## Warning: package 'expsmooth' was built under R version 3.5.3

library(seasonal)

## Warning: package 'seasonal' was built under R version 3.5.3

Question 1

autoplot(pigs)

sespigs = ses(pigs, h=4)
sespigs$model

## 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

autoplot(sespigs)

#Interval produced by R
sespigs$upper[1, "95%"]

##      95% 
## 119020.8

sespigs$lower[1, "95%"]

##      95% 
## 78611.97

#My intervals
s = sd(sespigs$residuals)
sespigs$mean[1] + 1.96*s

## [1] 118952.8

sespigs$mean[1] - 1.96*s

## [1] 78679.97

Question 5

5A

autoplot(books)

#The series has a definite upward trend, with daily seasonality.

5B

sespaper= ses(books[,"Paperback"], h=4)
seshard = ses(books[,"Hardcover"], h=4)
autoplot(sespaper)

autoplot(seshard)

#### 5C

spaper = sqrt(mean(sespaper$residuals^2))
shard = sqrt(mean(seshard$residuals^2))

Question 6

6A

holtpaper = holt(books[,"Paperback"], h=4)
holthard = holt(books[,"Hardcover"], h=4)

6B

hpaper = sqrt(mean(holtpaper$residuals^2))
hhard = sqrt(mean(holthard$residuals^2))
#Holt's method returned lower RMSE values

####6C

autoplot(holtpaper)

autoplot(sespaper)

autoplot(holthard)

autoplot(seshard)

#I prefer Holt's method, it captures the upward trend better

6D

#95% PI paperback sales Holt
holtpaper$upper[1, "95%"]

##      95% 
## 275.0205

holtpaper$lower[1, "95%"]

##     95% 
## 143.913

#95% PI paperback sales by formula
holtpaper$mean[1] +1.96*hpaper

## [1] 270.4951

holtpaper$mean[1] -1.96*hpaper

## [1] 148.4384

#95% PI hardcover sales by Holt
holthard$upper[1, "95%"]

##      95% 
## 307.4256

holthard$lower[1, "95%"]

##      95% 
## 192.9222

#95% PI hardcover sales by formula
holthard$mean[1] + 1.96*hhard

## [1] 303.4733

holthard$mean[1] - 1.96*hhard

## [1] 196.8745

Question 7

autoplot(eggs)

holtreg = holt(eggs, h=100)
holtdamped = holt(eggs, damped=TRUE, h=100)
holtdampexp = holt(eggs, damped=TRUE, exponential = TRUE, h=100)
holtexp = holt(eggs, damped=FALSE, exponential = TRUE, h=100)
autoplot(holtreg)

accuracy(holtreg)

##                      ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 0.04499087 26.58219 19.18491 -1.142201 9.653791 0.9463626
##                    ACF1
## Training set 0.01348202

autoplot(holtdamped)

accuracy(holtdamped)

##                     ME     RMSE     MAE       MPE     MAPE      MASE
## Training set -2.891496 26.54019 19.2795 -2.907633 10.01894 0.9510287
##                      ACF1
## Training set -0.003195358

autoplot(holtdampexp)

accuracy(holtdampexp)

##                      ME     RMSE      MAE       MPE     MAPE     MASE
## Training set -0.9089678 26.59113 19.54973 -2.125756 10.02328 0.964359
##                    ACF1
## Training set 0.01376115

autoplot(holtexp)

#Holt's method with exponential trend has the lowest RMSE

Question 10

10A

autoplot(ukcars)

#The series trends downward until around 1983, then trends upwards until 2000, where it takes a dip and then levels out.

10B

stlcars = stl(ukcars, s.window=4, robust=TRUE)
autoplot(stlcars)

seas.stlcars = seasadj(stlcars)
autoplot(seas.stlcars)

10C

add.damp.stl = stlf(seas.stlcars, etsmodel= "AAN", damped=TRUE, h=8)
autoplot(add.damp.stl)

10D

holtfc = stlf(seas.stlcars, etsmodel="AAN", damped=FALSE, h=8)
autoplot(holtfc)

10E

etsfc = ets(ukcars)
summary(etsfc)

## ETS(A,N,A) 
## 
## Call:
##  ets(y = ukcars) 
## 
##   Smoothing parameters:
##     alpha = 0.6199 
##     gamma = 1e-04 
## 
##   Initial states:
##     l = 314.2568 
##     s = -1.7579 -44.9601 21.1956 25.5223
## 
##   sigma:  25.9302
## 
##      AIC     AICc      BIC 
## 1277.752 1278.819 1296.844 
## 
## Training set error measures:
##                    ME     RMSE      MAE        MPE     MAPE      MASE
## Training set 1.313884 25.23244 20.17907 -0.1570979 6.629003 0.6576259
##                    ACF1
## Training set 0.02573334

#ETS used ANA method

10F

#Accuracy of additive damped trend method
accuracy(add.damp.stl)

##                   ME     RMSE      MAE      MPE     MAPE      MASE
## Training set 1.47663 21.29656 16.60585 0.167431 5.330547 0.5383996
##                     ACF1
## Training set 0.006062078

#Accuracy of holt method, not damped
accuracy(holtfc)

##                      ME     RMSE     MAE        MPE    MAPE      MASE
## Training set -0.2811218 21.31828 16.5201 -0.4191414 5.33054 0.5356195
##                     ACF1
## Training set 0.008478115

#Accuracy of ETS model auto-chosen
accuracy(etsfc)

##                    ME     RMSE      MAE        MPE     MAPE      MASE
## Training set 1.313884 25.23244 20.17907 -0.1570979 6.629003 0.6576259
##                    ACF1
## Training set 0.02573334

#Additive damped trend method gives the best RMSE

10G

autoplot(add.damp.stl)

autoplot(holtfc)

autoplot(forecast(etsfc, h=8))

#can't really compare the three, 2 are seasonally adjusted

10H

checkresiduals(add.damp.stl)

## Warning in checkresiduals(add.damp.stl): The fitted degrees of freedom is
## based on the model used for the seasonally adjusted data.

## 
##  Ljung-Box test
## 
## data:  Residuals from STL +  ETS(A,Ad,N)
## Q* = 46.556, df = 3, p-value = 4.32e-10
## 
## Model df: 5.   Total lags used: 8

checkresiduals(holtfc)

## Warning in checkresiduals(holtfc): The fitted degrees of freedom is based
## on the model used for the seasonally adjusted data.

## 
##  Ljung-Box test
## 
## data:  Residuals from STL +  ETS(A,A,N)
## Q* = 44.201, df = 4, p-value = 5.828e-09
## 
## Model df: 4.   Total lags used: 8

checkresiduals(forecast(etsfc, h=8))

## 
##  Ljung-Box test
## 
## data:  Residuals from ETS(A,N,A)
## Q* = 18.324, df = 3, p-value = 0.0003772
## 
## Model df: 6.   Total lags used: 9

Question 11

11A

autoplot(visitors)

#Upward trend line, looks like the seasonality increases propotional to the level

11B

test = window(visitors, start = c(2003,4))
train = window(visitors, end = c(2004, 3))
hwmult = hw(train, seasonal="multiplicative", h=24)
autoplot(hwmult)

#11C Multiplicative is needed because the seasonal variation in the data increases as the level of the series increases.

11D

etsfc <- train %>% ets() %>% forecast(h=24)
autoplot(etsfc)

etsadd.BCfc <- train %>% ets(lambda = BoxCox.lambda(train), additive.only=TRUE) %>% forecast(h=24)
autoplot(etsadd.BCfc)

snaive = snaive(train, h=24)
autoplot(snaive)

BCstl.ets <- train %>% stlm(
  lambda = BoxCox.lambda(train), 
  s.window = 13, 
  robust=TRUE, 
  method="ets") %>% forecast(h=24)
autoplot(BCstl.ets)

11E

accuracy(hwmult, test)

##                     ME     RMSE      MAE       MPE     MAPE      MASE
## Training set -1.557919 14.70228 11.17152 -1.019484 4.452697 0.4268183
## Test set     15.090935 20.85700 16.04326  3.122066 3.359691 0.6129476
##                    ACF1 Theil's U
## Training set 0.09876254        NA
## Test set     0.21025230 0.2869292

accuracy(etsfc, test)

##                     ME     RMSE      MAE       MPE     MAPE      MASE
## Training set  1.220635 14.90665 10.78433 0.2197614 4.018479 0.4120252
## Test set     31.260025 36.53961 31.26002 6.6962946 6.696295 1.1943180
##                     ACF1 Theil's U
## Training set -0.01726088        NA
## Test set      0.39874182 0.5090901

accuracy(etsadd.BCfc, test)

##                     ME     RMSE      MAE        MPE     MAPE      MASE
## Training set -0.373716 15.55474 11.35863 -0.1677915 4.070118 0.4339668
## Test set     12.013636 19.11730 15.31311  2.3827550 3.229648 0.5850514
##                     ACF1 Theil's U
## Training set -0.05004163        NA
## Test set      0.23875318 0.2521544

accuracy(snaive, test)

##                    ME     RMSE      MAE       MPE     MAPE     MASE
## Training set 16.78326 31.24472 26.17395  6.837756 10.14381 1.000000
## Test set     48.30000 55.23646 48.30000 11.417380 11.41738 1.845346
##                   ACF1 Theil's U
## Training set 0.6512811        NA
## Test set     0.6747630 0.7596856

accuracy(BCstl.ets, test)

##                      ME     RMSE       MAE        MPE     MAPE      MASE
## Training set -0.8379729 13.88885  9.778233 -0.2528935 3.536527 0.3735864
## Test set      9.4686802 17.54139 13.069371  1.8664731 2.747569 0.4993274
##                     ACF1 Theil's U
## Training set -0.04339136        NA
## Test set     -0.04401040 0.2405589

#The Box-Cox STL ETS model performed the best by RMSE

Question 12

fets <- function(y, h) {
  forecast(ets(y), h = h)
}
#fets(qcement, h=4) %>% tsCV()
#qcement %>% snaive() %>% forecast(h=4) %>% tsCV()
#couldn't figure this out, kept throwing error "invalid time series parameters specified"

Chapter 8

Question 2

autoplot(ibmclose)

plot(acf(ibmclose))

plot(pacf(ibmclose))

#Autocorrelation is clearly present as shown in the ACF plot as the CV is broken at every lag. The PACF plot only shows autocorrelation once at the beginning of the data.

Question 3

ndiffs(usnetelec)

## [1] 1

ndiffs(usgdp)

## [1] 2

ndiffs(mcopper)

## [1] 1

ndiffs(enplanements)

## [1] 1

ndiffs(visitors)

## [1] 1

Question 6

6A

y <- ts(numeric(100))
e <- rnorm(100)
for(i in 2:100)
  y[i] <- 0.6*y[i-1] + e[i]
plot(y[i])

6B

ar1 <- function(phi1){
  y <- ts(numeric(100))
  e <- rnorm(100)
  for(i in 2:100){
    y[i] <- phi1*y[i-1] + e[i]
  }
  return(y)
}

autoplot(ar1(0.3), series = "0.3") +
  geom_line(size = 1, colour = "red") +
  autolayer(y, series = "0.6", size = 1) +
  autolayer(ar1(0.9), size = 1, series = "0.9") +
  ylab("AR(1) models") +
  guides(colour = guide_legend(title = "Phi1"))

#The variation in y increases as phi increases.

6C

ma1 <- function(theta1){
  y <- ts(numeric(100))
  e <- rnorm(100)
  for(i in 2:100){
    y[i] <- theta1*e[i-1] + e[i]
  }
  return(y)
}

6D

autoplot(ma1(0.3), series = "0.3") +
  geom_line(size = 1, colour = "red") +
  autolayer(y, series = "0.6", size = 1) +
  autolayer(ma1(0.9), size = 1, series = "0.9") +
  ylab("MA(1) models") +
  guides(colour = guide_legend(title = "Theta1"))

#as theta increases, variation of y increases, just like phi.

6E

MYarima11 <- ts(numeric(50))
e <- rnorm(50)
for(i in 2:50){
  MYarima11[i] <- 0.6*MYarima11[i-1] + 0.6*e[i-1] + e[i]
}

6F

MYarima2 <- ts(numeric(50))
e <- rnorm(50)
for(i in 3:50){
  MYarima2[i] <- -0.8*MYarima2[i-1] + 0.3*MYarima2[i-2] + e[i]
}

6G

autoplot(MYarima11, series = "ARMA(1, 1)") +
  autolayer(MYarima2, series = "AR(2)") +
  ylab("y") +
  guides(colour = guide_legend(title = "ARIMA Method"))

#AR(2) is nonstationary data, while ARMA(1,1) is stationary data.

Question 7

7A

wmurderarima = auto.arima(wmurders)

7D

checkresiduals(wmurderarima)

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

#Residuals look good

7E

arimafc = forecast(wmurderarima, h=3)
arimafc$mean

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

7F

autoplot(arimafc)

Question 8

8A

arima1 = auto.arima(austa)
checkresiduals(arima1)

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

arima1fc = forecast(arima1, h=10)
autoplot(arima1fc)

8B

arima2 = Arima(austa, order= c(0,1,1), include.drift = F)
arima2fc = forecast(arima2, h=10)
autoplot(arima2fc)

autoplot(arima1fc)

#The model without drift did not capture the upward trend, similar to a naive model.

8C

arima3 = Arima(austa, order= c(2,1,3), include.drift=T)
arima3fc = forecast(arima3, h=10)
autoplot(arima3fc)

#can't remove constant, throws error "non-stationary AR part from CSS"

8D

arima4 = Arima(austa, order = c(0,0,1))
arima4fc = forecast(arima4, h=10)
autoplot(arima4fc)

8E

arima5 = Arima(austa, order=c(0,2,1))
arima5fc = forecast(arima5, h=10)
autoplot(arima5fc)

Question 9

usgdparima = auto.arima(usgdp)
usgdparima

## Series: usgdp 
## ARIMA(2,2,2) 
## 
## Coefficients:
##           ar1     ar2      ma1      ma2
##       -0.1228  0.3106  -0.5835  -0.3669
## s.e.   0.2869  0.0872   0.3004   0.2862
## 
## sigma^2 estimated as 1604:  log likelihood=-1199.57
## AIC=2409.13   AICc=2409.39   BIC=2426.43

usgdpets = ets(usgdp)
usgdpets

## ETS(A,A,N) 
## 
## Call:
##  ets(y = usgdp) 
## 
##   Smoothing parameters:
##     alpha = 0.9999 
##     beta  = 0.278 
## 
##   Initial states:
##     l = 1557.4589 
##     b = 18.6862 
## 
##   sigma:  41.8895
## 
##      AIC     AICc      BIC 
## 3072.303 3072.562 3089.643

checkresiduals(usgdparima)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(2,2,2)
## Q* = 8.6247, df = 4, p-value = 0.0712
## 
## Model df: 4.   Total lags used: 8

usgdparima %>% forecast(h=8) %>% autoplot()

usgdpets %>% forecast(h=8) %>% autoplot

Question 10

autoplot(austourists)

#A: increasing trend with seasonality
acf(austourists)

#B: Autocorrelation is frequent in the beginning of the data but decreases over time
pacf(austourists)

#C: 5 major spikes in the dataset
ggtsdisplay(diff(austourists, lag=4))

auto.arima(austourists)

## Series: austourists 
## ARIMA(1,0,0)(1,1,0)[4] with drift 
## 
## Coefficients:
##          ar1     sar1   drift
##       0.4705  -0.5305  0.5489
## s.e.  0.1154   0.1122  0.0864
## 
## sigma^2 estimated as 5.15:  log likelihood=-142.48
## AIC=292.97   AICc=293.65   BIC=301.6

Question 11

mausmelec = ma(usmelec, order=12, centre = TRUE)
autoplot(mausmelec)

ndiffs(usmelec)

## [1] 1

elecarima = auto.arima(usmelec)
elecarima

## Series: usmelec 
## ARIMA(1,0,2)(0,1,1)[12] with drift 
## 
## Coefficients:
##          ar1      ma1      ma2     sma1   drift
##       0.9717  -0.4374  -0.2774  -0.7061  0.3834
## s.e.  0.0163   0.0483   0.0493   0.0310  0.0868
## 
## sigma^2 estimated as 57.67:  log likelihood=-1635.13
## AIC=3282.26   AICc=3282.44   BIC=3307.22

checkresiduals(elecarima)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,0,2)(0,1,1)[12] with drift
## Q* = 42.725, df = 19, p-value = 0.001413
## 
## Model df: 5.   Total lags used: 24

#normal distribution of residuals and low autocorrelation indicate white noise
elecfc = forecast(elecarima, h= 12*15)
autoplot(elecfc)

Question 12

mcopper

##         Jan    Feb    Mar    Apr    May    Jun    Jul    Aug    Sep    Oct
## 1960  255.2  259.7  249.3  258.0  244.3  246.8  250.6  241.3  231.0  218.7
## 1961  216.6  220.1  221.7  225.6  238.6  232.8  226.1  227.2  225.8  225.0
## 1962  226.8  231.3  231.1  230.6  230.5  230.4  230.4  230.4  230.4  230.6
## 1963  230.6  230.6  230.6  230.6  230.6  230.6  230.6  230.6  230.6  230.6
## 1964  234.2  247.8  266.1  307.7  295.8  288.8  305.4  356.9  415.0  485.1
## 1965  356.7  419.2  440.6  480.5  490.8  466.1  404.0  431.6  473.5  500.1
## 1966  599.0  668.7  668.7  679.9  592.8  604.8  559.5  426.4  402.4  454.9
## 1967  443.7  435.4  391.8  354.9  369.5  362.3  355.9  372.7  378.3  406.3
## 1968  586.4  715.8  708.0  522.4  456.5  473.0  439.9  439.9  462.3  449.7
## 1969  523.8  535.0  532.4  578.1  579.6  617.0  605.4  669.0  657.6  647.4
## 1970  677.2  690.1  730.4  725.2  665.9  607.0  567.8  527.5  519.3  475.8
## 1971  420.7  424.9  476.5  521.2  464.1  447.2  464.4  451.0  427.5  417.7
## 1972  418.9  426.9  442.2  433.5  423.3  412.4  423.2  427.0  434.2  428.7
## 1973  474.7  511.8  610.2  638.8  613.0  678.4  795.4  843.5  799.6  849.6
## 1974  913.2 1006.5 1172.0 1267.7 1190.6 1020.1  802.8  767.9  630.1  599.2
## 1975  512.4  528.7  554.5  560.5  539.7  522.5  559.3  603.7  580.2  573.1
## 1976  587.8  601.6  683.7  817.9  836.0  878.0  922.1  862.2  844.9  784.8
## 1977  814.8  833.4  881.0  830.4  797.9  763.1  724.4  665.8  685.6  682.9
## 1978  651.2  627.0  657.6  694.4  716.0  725.0  705.4  734.6  735.8  750.1
## 1979  827.1  969.7 1005.5 1012.3  935.3  888.8  802.4  883.2  953.6  967.6
## 1980 1147.9 1220.2 1045.5  937.1  888.9  858.5  916.7  877.9  857.6  846.2
## 1981  777.4  784.9  814.3  837.0  834.0  861.0  897.3  981.4  942.2  904.9
## 1982  853.6  863.8  836.3  858.5  843.6  740.4  829.9  841.0  832.7  861.4
## 1983  997.3 1074.3 1071.9 1090.0 1122.8 1098.6 1115.6 1091.0 1041.0  958.4
## 1984  978.0  991.5 1030.9 1078.1 1022.3  991.2 1007.9 1018.2 1029.3 1043.5
## 1985 1205.5 1269.8 1234.7 1212.7 1225.3 1118.0 1067.5 1025.7 1001.4  973.8
## 1986  995.4  982.8  984.4  957.3  932.3  937.0  891.7  876.5  916.1  922.9
## 1987  893.8  902.7  920.1  909.3  911.9  964.6 1052.8 1097.5 1100.7 1182.7
## 1988 1476.8 1324.3 1286.3 1216.5 1307.1 1428.8 1297.5 1296.2 1445.7 1689.4
## 1989 1912.6 1765.0 1904.1 1832.3 1679.1 1638.8 1539.1 1731.5 1834.9 1801.8
## 1990 1432.7 1390.9 1615.7 1639.6 1633.7 1510.5 1529.6 1554.4 1611.6 1409.6
## 1991 1265.1 1246.5 1322.5 1412.3 1336.8 1344.8 1353.8 1325.6 1346.1 1371.2
## 1992 1182.2 1240.5 1291.7 1260.9 1224.6 1239.1 1313.8 1297.2 1307.0 1360.3
## 1993 1472.3 1536.7 1472.3 1261.9 1159.0 1228.5 1287.5 1305.0 1219.7 1108.7
## 1994 1208.4 1261.0 1284.0 1268.1 1430.5 1549.9 1589.5 1560.0 1601.0 1586.2
## 1995 1910.5 1830.3 1826.8 1805.7 1745.9 1876.8 1927.4 1936.6 1870.3 1782.7
## 1996 1708.7 1650.2 1676.5 1713.5 1757.3 1407.7 1276.6 1294.9 1244.3 1236.2
## 1997 1469.2 1480.3 1506.4 1466.8 1540.0 1588.0 1466.7 1403.6 1315.4 1256.4
## 1998 1032.1 1014.4 1051.4 1078.1 1057.8 1005.6 1004.2  991.8  979.1  935.6
## 1999  866.7  866.5  849.6  910.2  935.6  891.3 1041.2 1025.5 1077.3 1040.1
## 2000 1124.2 1124.7 1100.6 1059.8 1183.3 1161.6 1192.3 1245.1 1365.3 1308.1
## 2001 1209.8 1214.9 1202.8 1159.3 1178.9 1147.5 1078.5 1019.0  974.4  948.5
## 2002 1049.9 1097.4 1127.9 1101.4 1092.9 1110.0 1022.2  962.4  950.0  952.5
## 2003 1018.8 1046.7 1047.4 1007.9 1015.5 1015.4 1054.0 1103.4 1109.1 1143.9
## 2004 1329.0 1477.5 1647.3 1635.6 1529.5 1469.7 1523.4 1563.1 1614.9 1666.7
## 2005 1687.8 1723.8 1773.0 1790.0 1751.7 1938.0 2063.8 2115.8 2133.3 2301.1
## 2006 2677.7 2851.6 2927.3 3610.9 4306.0 3897.7 4170.5 4059.8 4032.9 3998.6
##         Nov    Dec
## 1960  222.8  227.5
## 1961  225.7  226.3
## 1962  230.4  230.5
## 1963  230.6  232.2
## 1964  500.1  453.3
## 1965  523.8  541.4
## 1966  464.4  433.4
## 1967  515.3  551.3
## 1968  458.0  493.4
## 1969  679.4  707.8
## 1970  451.9  435.4
## 1971  406.4  410.7
## 1972  427.8  435.7
## 1973  950.4  959.9
## 1974  608.2  553.2
## 1975  575.1  568.9
## 1976  780.3  767.3
## 1977  650.2  679.1
## 1978  749.3  771.7
## 1979  978.3 1005.6
## 1980  839.5  800.4
## 1981  867.7  869.3
## 1982  884.3  911.4
## 1983  940.4  986.7
## 1984 1085.1 1113.5
## 1985  950.9  962.4
## 1986  915.0  927.6
## 1987 1419.9 1566.6
## 1988 1826.2 1915.3
## 1989 1647.3 1514.1
## 1990 1315.9 1292.7
## 1991 1336.8 1216.4
## 1992 1413.2 1422.5
## 1993 1100.3 1156.4
## 1994 1763.6 1914.4
## 1995 1904.8 1898.7
## 1996 1343.7 1366.9
## 1997 1135.0 1061.5
## 1998  946.9  881.9
## 1999 1065.3 1093.7
## 2000 1258.9 1264.4
## 2001  994.2 1020.8
## 2002 1006.4 1005.7
## 2003 1216.0 1258.5
## 2004 1678.5 1631.1
## 2005 2461.6 2621.7
## 2006 3675.9 3435.8

mcopperarima = auto.arima(mcopper)
mcopperets = ets(mcopper)
mcopperarima

## Series: mcopper 
## ARIMA(0,1,1) 
## 
## Coefficients:
##          ma1
##       0.2900
## s.e.  0.0419
## 
## sigma^2 estimated as 6026:  log likelihood=-3248.53
## AIC=6501.07   AICc=6501.09   BIC=6509.73

mcopperets

## ETS(M,Ad,N) 
## 
## Call:
##  ets(y = mcopper) 
## 
##   Smoothing parameters:
##     alpha = 0.9999 
##     beta  = 0.2141 
##     phi   = 0.8 
## 
##   Initial states:
##     l = 265.0541 
##     b = -3.9142 
## 
##   sigma:  0.0632
## 
##      AIC     AICc      BIC 
## 8052.038 8052.189 8078.049

checkresiduals(mcopperarima)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,1)
## Q* = 30.716, df = 23, p-value = 0.1299
## 
## Model df: 1.   Total lags used: 24

copperfc = forecast(mcopperarima, h=24)
autoplot(copperfc)

mcopperets %>% forecast(h=24) %>% autoplot

### Question 13

ndiffs(auscafe)

## [1] 1

auto.arima(auscafe) -> ausarima
ausarima

## Series: auscafe 
## ARIMA(0,1,1)(0,1,1)[12] 
## 
## Coefficients:
##           ma1     sma1
##       -0.3673  -0.5991
## s.e.   0.0443   0.0348
## 
## sigma^2 estimated as 0.001673:  log likelihood=732.46
## AIC=-1458.92   AICc=-1458.87   BIC=-1446.85

checkresiduals(ausarima)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,1)(0,1,1)[12]
## Q* = 64.68, df = 22, p-value = 4.421e-06
## 
## Model df: 2.   Total lags used: 24

ausarima %>% forecast(h=24) %>% autoplot()

auscafe %>% ets %>% forecast(h=24) %>% autoplot()