2.
(i) The auto.arima() function in R uses a variation of the Hyndman and Khandakar algorithm which combines unit root tests, minimization of the AICc and MLE to obtain an ARIMA model. The algorithm follows these steps.
1. The number of differences d is determined using repeated ADF tests.
2. The values of p and q are then chosen by minimizing the AICc after differencing the data d times. Rather than considering every possible combination of p and q, the algorithm uses a stepwise search to traverse the model space.
(a) The best model (with smallest AICc) is selected from the following four:
ARIMA(2,d,2),
ARIMA(0,d,0),
ARIMA(1,d,0),
ARIMA(0,d,1).
If d=0 then the constant cc is included; if d≥1 then the constant c is set to zero. This is called the “current model”.
(b) Variations on the current model are considered: vary p and/or q from the current model by ±1; include/exclude c from the current model. The best model considered so far (either the current model, or one of these variations) becomes the new current model.
(c) Repeat Step 2(b) until no lower AICc can be found.
(iii)
3.
library(tseries)
#Fit AR(1) and AR(2) models
fit1<-arima(myts,order=c(1,0,0),method='ML')
fit2<-arima(myts,order=c(2,0,0),method='ML')
#run the test
adf.test(fit1$residuals)
## Warning in adf.test(fit1$residuals): p-value smaller than printed p-value
##
## Augmented Dickey-Fuller Test
##
## data: fit1$residuals
## Dickey-Fuller = -4.7067, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary
adf.test(fit2$residuals)
## Warning in adf.test(fit2$residuals): p-value smaller than printed p-value
##
## Augmented Dickey-Fuller Test
##
## data: fit2$residuals
## Dickey-Fuller = -4.2581, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary
##both p-values are less than 0.01 so accept alternative hypothesis of stationary.
4.
#(a)
#same as above that using auto.arima
airpass.v<-scan("airpass.tsm")
airpass<-ts(log(airpass.v[1:(length(airpass.v)-12)]),start=1949,frequency=12)
fit<-auto.arima(airpass,d=1, ic="aicc", test="adf")
#(b)
tsdiag(fit)
## the result looks fine
#(c)
airpass.full<-ts(log(airpass.v),start=1949,frequency=12)
pred<-predict(fit, n.ahead=12)$pred
s.pred<-predict(fit, n.ahead=12, se.fit=T)$se
ts.plot(airpass.full, type="l", xlim=c(1955,1961), ylim=c(5.5,7.0))
t<-seq(1960-1/12, 1961-1/12,1/12)
lines(t, c(airpass[length(airpass)], pred), type="l", lty=2)
lines(t, c(airpass[length(airpass)], pred+1.96*s.pred), type="l", lty=3)
lines(t, c(airpass[length(airpass)], pred-1.96*s.pred), type="l", lty=4)
##Predictions are quite near the actual data for the next year.
##The confidence limit gets wider, but is still smaller than the annual variation at the end.
#(d)
cbind(pred, pred-1.96*s.pred, pred+1.96*s.pred)
## pred pred - 1.96 * s.pred pred + 1.96 * s.pred
## Jan 1960 6.038647 5.966964 6.110330
## Feb 1960 5.988763 5.903207 6.074319
## Mar 1960 6.145428 6.047954 6.242902
## Apr 1960 6.118993 6.010907 6.227079
## May 1960 6.159652 6.041907 6.277397
## Jun 1960 6.304666 6.177996 6.431336
## Jul 1960 6.433288 6.298282 6.568294
## Aug 1960 6.445958 6.303101 6.588815
## Sep 1960 6.266719 6.116422 6.417017
## Oct 1960 6.136192 5.978805 6.293579
## Nov 1960 6.007899 5.843727 6.172070
## Dec 1960 6.114338 5.943652 6.285023
#(f)
len<-length(airpass.full)
index<-(len-11):len
cbind(airpass.full[index],airpass.full[index]-pred)
## airpass.full[index] airpass.full[index] - pred
## Jan 1960 6.033086 -0.0055606120
## Feb 1960 5.968708 -0.0200554916
## Mar 1960 6.037871 -0.1075567720
## Apr 1960 6.133398 0.0144046681
## May 1960 6.156979 -0.0026728622
## Jun 1960 6.282267 -0.0223994376
## Jul 1960 6.432940 -0.0003481843
## Aug 1960 6.406880 -0.0390781075
## Sep 1960 6.230481 -0.0362380037
## Oct 1960 6.133398 -0.0027939616
## Nov 1960 5.966147 -0.0417517659
## Dec 1960 6.068426 -0.0459122744
##The last real value is inside the forecast confidence interval.