HW9

library(TSA)

## 
## Attaching package: 'TSA'

## The following objects are masked from 'package:stats':
## 
##     acf, arima

## The following object is masked from 'package:utils':
## 
##     tar

library(forecast)

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

## Registered S3 methods overwritten by 'forecast':
##   method       from
##   fitted.Arima TSA 
##   plot.Arima   TSA

library(tseries)

data("airpass")
plot(airpass, type='o',ylab='Air Passengers')

plot(log(airpass), type='o',ylab='Log(Air Passengers)')

plot(diff(log(airpass)),type='o',ylab='Difference of Log(Air Passengers)')

plot(diff(log(airpass)),type='l',ylab='Difference of Log(Air Passengers)')
points(diff(log(airpass)),x=time(diff(log(airpass))),pch=as.vector(season(diff(log(airpass)))))

plot(diff(diff(log(airpass)),lag=12),type='l',
     ylab='First & Seasonal Differences of Log(AirPass)')
points(diff(diff(log(airpass)),lag=12),x=time(diff(diff(log(airpass)),lag=12)),
       pch=as.vector(season(diff(diff(log(airpass)),lag=12))))

acf(as.vector(diff(diff(log(airpass)),lag=12)),ci.type='ma',
    main='First & Seasonal Differences of Log(AirPass)')

model=arima(log(airpass),order=c(0,1,1),seasonal=list(order=c(0,1,1),period=12))
model

## 
## Call:
## arima(x = log(airpass), order = c(0, 1, 1), seasonal = list(order = c(0, 1, 
##     1), period = 12))
## 
## Coefficients:
##           ma1     sma1
##       -0.4018  -0.5569
## s.e.   0.0896   0.0731
## 
## sigma^2 estimated as 0.001348:  log likelihood = 244.7,  aic = -485.4

par(mar = c(1, 1, 1, 1))
tsdiag(model)

hist(residuals(model),xlab='Residuals')
#qqplot(residuals(model))
qqline(residuals(model))

shapiro.test(residuals(model))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(model)
## W = 0.98637, p-value = 0.1674

forecasts <- predict(model,n.ahead = 24)
prediction <- predict(model, n.ahead = 24)$pred
upper_lim <- prediction + 2*forecasts$se
lower_lim <- prediction - 2*forecasts$se
ts.plot(log(airpass),prediction,upper_lim, lower_lim, col=c("black","blue",'red','red'),lty = c(1,1,3,3) ,lwd=2)
legend("topleft",c("Actual","Forecast","95% CI"),col=c("black","blue",'red'),lty = c(1,1,3))

ts.plot(airpass,exp(prediction),exp(upper_lim), exp(lower_lim), col=c("black","blue",'red','red'),lty = c(1,1,3,3
) ,lwd=2)
legend("topleft",c("Actual","Forecast","95% CI"),col=c("black","blue",'red'),lty = c(1,1,3))