# forecasting:  predict(arima) uses Kalman filter 

ts.plot(AirPassengers)
airpass.log<-log(AirPassengers)



airpass.arima3<-arima(airpass.log, order=c(0,1,1), seasonal=list(order=c(0,1,1),
period=12 )  ) 
tsdiag(airpass.arima3) 

years<-seq(1949, 1960+11/12, by=1/12) 

# take the last year out for the forecast 

airpass.short<-airpass.log[1:132]
years.short<-years[1:132] 

plot.ts(airpass.short) 




par(mfrow=c(1,2))

acf(airpass.short)
pacf(airpass.short) 


# differencing

airpass.short.diff<-diff(airpass.short)
ts.plot(airpass.short.diff) 

airpass.short.diff2<-diff(airpass.short.diff, lag=12)
ts.plot(airpass.short.diff2) 

# fit SARIMA(0,1,1) x (0,1,1) 

airpass.short.arima<-arima(airpass.short, order=c(0,1,1), seasonal=list(order=c(0,1,1),
period=12 )  ) 
tsdiag(airpass.short.arima) 

fit.res<-airpass.short.arima$residuals 

# are consistent with white noise 

cpgram(fit.res) 


# forecasting 

fyears<-years[133:144]
pred.1<-predict(airpass.short.arima, 12)
pred.1

forecast.1<-pred.1$pred
flimits.1<-1.96* round(pred.1$se, 3) 



plot(years.short, airpass.short, type="l", xlab="year", ylab="log passenger numbers", xlim=range(years))
lines(fyears, forecast.1, lwd=2) 

lines(fyears, forecast.1 + flimits.1, lty=2)
lines(fyears, forecast.1 - flimits.1, lty=2) 

# how good is the forecast? 


plot(airpass.log, type="l", xlab="year", ylab="log passenger numbers", xlim=range(years))
lines(fyears, forecast.1, col="red") 

lines(fyears, forecast.1 + flimits.1, lty=2, col="blue")
lines(fyears, forecast.1 - flimits.1, lty=2, col="blue") 




# calculate the mean squared prediction error 

mse<- mean(( 10^forecast.1 - 10^airpass.log[133:144])^2) 

mse 

# we could now compare this for different models.