install.packages("forecast") 
library(forecast) 

par(mfrow=c(2,4))
plot(AirPassengers)
plot(diff(AirPassengers),ylab="Diff 1",main=expression(y[t]-y[t-1]))
plot(diff(AirPassengers,12),ylab="Diff 12",main=expression(y[t]-y[t-12]))
plot(diff(diff(AirPassengers,1),12),ylab="Diff 1,12",main=expression(y[t]-y[t-1]-y[t-12]+y[t-13]))

plot(log(AirPassengers))
plot(diff(log(AirPassengers)),ylab="Diff 1",main=expression(y[t]-y[t-1]))
plot(diff(log(AirPassengers),12),ylab="Diff 12",main=expression(y[t]-y[t-12]))
plot(diff(diff(log(AirPassengers),1),12),ylab="Diff 1,12",main=expression(y[t]-y[t-1]-y[t-12]+y[t-13]))

# Training sample: 1949.1 - 1957.12
# Testing sample: 1958.1 - 1960.12
n1 = 108 
n = length(AirPassengers)
trend = 1:n
d1 = rep(c(1,rep(0,11)),12)
d2 = c(0,d1)[1:n]
d3 = c(0,d2)[1:n]
d4 = c(0,d3)[1:n]
d5 = c(0,d4)[1:n]
d6 = c(0,d5)[1:n]
d7 = c(0,d6)[1:n]
d8 = c(0,d7)[1:n]
d9 = c(0,d8)[1:n]
d10 = c(0,d9)[1:n]
d11 = c(0,d10)[1:n]
X = cbind(1,trend,d1,d2,d3,d4,d5,d6,d7,d8,d9,d10,d11)


par(mfrow=c(2,3))
# Deterministic seasonality - AirPassengers
###################################
reg = lm(AirPassengers[1:n1]~X[1:n1,]-1)
fit = ts(reg$fit,start=c(1949,1),end=c(1957,12),frequency=12)
plot(AirPassengers,ylim=range(AirPassengers,fit))
lines(fit,col=2)
legend("topleft",legend=c("Data","Fit"),col=1:2,lty=1)
title("Deterministic cycle")

# Deterministic seasonality - log-AirPassengers
######################################
reg = lm(log(AirPassengers[1:n1])~X[1:n1,]-1)
fit = ts(reg$fit,start=c(1949,1),end=c(1957,12),frequency=12)
plot(log(AirPassengers),ylim=range(log(AirPassengers),fit))
lines(fit,col=2)
legend("topleft",legend=c("Data","Fit"),col=1:2,lty=1)
title("Deterministic cycle")

# Seasonal ARIMA modeling of AirPassengers
######################################
bics = matrix(0,16,5)
s = 0
for (i in 0:1)
for (j in 0:1)
for (k in 0:1)
for (l in 0:1){
 s = s + 1
bics[s,] = c(i,j,k,l,Arima(window(AirPassengers,end=1957+11/12), 
                      order=c(i,1,j),seasonal=list(order=c(k,1,l),period=12), lambda=0)$bic)
}
bics[bics[,5]==min(bics[,5]),]

# Fitting the best model
air.model = Arima(window(AirPassengers,end=1957+11/12), order=c(0,1,1),
                              seasonal=list(order=c(0,1,1),period=12), lambda=0) 
air.model
plot(AirPassengers,ylim=range(AirPassengers,air.model$fit))
lines(air.model$fit,col=2)
legend("topleft",legend=c("Data","Fit"),col=1:2,lty=1)
title("ARIMA(0,1,1)(0,1,1)[12]")
sarima1 = air.model


# Forecasting exercise
##################
reg = lm(AirPassengers[1:n1]~X[1:n1,]-1)
fit = ts(reg$fit,start=c(1949,1),end=c(1957,12),frequency=12)
plot(AirPassengers,ylim=range(AirPassengers,fit))
lines(fit,col=2)
fit1 = ts(X[(n1+1):n,]%*%reg$coef,start=c(1958,1),end=c(1960,12),frequency=12)
lines(fit1,col=4)
title("Deterministic cycle")


reg = lm(log(AirPassengers[1:n1])~X[1:n1,]-1)
fit = ts(reg$fit,start=c(1949,1),end=c(1957,12),frequency=12)
fit1 = ts(X[(n1+1):n,]%*%reg$coef,start=c(1958,1),end=c(1960,12),frequency=12)
plot(log(AirPassengers),ylim=range(log(AirPassengers),fit,fit1))
lines(fit,col=2)
lines(fit1,col=4)
title("Deterministic cycle")


air.model = Arima(window(AirPassengers,end=1957+11/12), order=c(0,1,1),
                              seasonal=list(order=c(0,1,1),period=12), lambda=0) 
plot(forecast(air.model,h=36),ylab="AirPassengers")
lines(AirPassengers)
sarima2 = air.model