# Exercise 6.7

airpass.vec=scan("http://www.uio.no/studier/emner/matnat/math/STK4060/v06/airpass.tsm")

airpass<-ts(log(airpass.vec[1:(length(airpass.vec)-12)]),start=1949,frequency=12)

# Plot the original data:
plot.ts(exp(airpass),xlab="time", ylab="airpass")
# Plot the log-tranformed data:
plot.ts(airpass,xlab="time", ylab="log(airpass)")

acf(airpass)


# I did not find any code for the tests described ins ection 6.3
# so I'll use visual inspection to choose the number of 
# differenciations
ts.plot(diff(airpass))
acf(diff(airpass),lag=40)
# Seems to have a seasonal term, but else, ok.


ts.plot(diff(airpass, lag=12))
acf(diff(airpass, lag=12),lag=40)
# Slow but steady decline of the acf; possible AR-model?


# See which of these two candidates for differenciation
# that gives the smallest model from the BIC criterion

# We'll now try both seasonal and normal ARMA models
# for I1=1 and I12=0:
zero.I1.model=arima(airpass,order=c(0,1,0))
best.I1.model=zero.I1.model
bicbest.I1.model=zero.I1.model
n=length(airpass)
best.I1.bic=-2*zero.I1.model$loglik+2*log(n)
bicbest.I1.p=0
bicbest.I1.q=0
bicbest.I1.P=0
bicbest.I1.Q=0
for(Q in 0:3)
{
  for(P in 0:3)
    {
      for(q in 0:3)
		 {
		   for(p in 0:3)
		     {
		       model<-arima(airpass,order=c(p,1,q),seasonal=c(P,0,Q),method="ML")
		       curr.bic=-2*model$loglik+log(n)*(P+Q+p+q+1+1)
		       if(curr.bic<best.I1.bic)
		 		 {
		 		   best.I1.bic=curr.bic
		 		   bicbest.I1.model=model
		 		   bicbest.I1.p=p
		 		   bicbest.I1.q=q
		 		   bicbest.I1.P=P
		 		   bicbest.I1.Q=Q
		 		   
		 		   show(c("BIC:",P,Q,p,q,best.I1.bic))
		 		 }
		     }
		 }
    }
}

bicbest.I1.model
#Coefficients:
#          ar1    sar1    sar2
#      -0.3068  0.5314  0.4333
#s.e.   0.0852  0.0815  0.0842
# 
#sigma^2 estimated as 0.001405:  log likelihood = 229.08,  aic = -450.15
# AR(1) + seasonal AR(2)




# Try seasonal differenciation, too...
zero.I12.model=arima(airpass,seasonal=c(0,1,0))
best.I12.model=zero.I12.model
bicbest.I12.model=zero.I12.model
n=length(I12)
best.I12.bic=-2*zero.I12.model$loglik+2*log(n)
bicbest.I12.p=0
bicbest.I12.q=0
bicbest.I12.P=0
bicbest.I12.Q=0
for(Q in 0:3)
{
  for(P in 0:3)
    {
      for(q in 0:3)
		 {
		   for(p in 0:3)
		     {
		       model<-arima(airpass,order=c(p,0,q),seasonal=c(P,1,Q),method="ML")
		       curr.bic=-2*model$loglik+log(n)*(P+Q+p+q+1+1)
		       if(curr.bic<best.I12.bic)
		 		 {
		 		   best.I12.bic=curr.bic
		 		   bicbest.I12.model=model
		 		   bicbest.I12.p=p
		 		   bicbest.I12.q=q
		 		   bicbest.I12.P=P
		 		   bicbest.I12.Q=Q
		 		   
		 		   # Watch the progress:
		 		   show(c("BIC:",P,Q,p,q,best.I12.bic))
		 		 }
		     }
		 }
    }
}

bicbest.I12.model
#Coefficients:
#         ar1     ar2     ma1     ma2      ma3     sar1     sar2     sma1
#      0.4096  0.5838  0.2236  -0.182  -0.0332  -0.1008  -0.0294  -0.4766
#s.e.     NaN     NaN     NaN     NaN      NaN   0.2407   0.1603   0.2277
#sigma^2 estimated as 0.001297:  log likelihood = 225.06,  aic = -432.11



# Best model from I1 seems better than that from I12
# We go for this model:
airpass.mod=bicbest.I1.model
airpass.mod
#Coefficients:
#          ar1    sar1    sar2
#      -0.3068  0.5314  0.4333
#s.e.   0.0852  0.0815  0.0842
#sigma^2 estimated as 0.001405:  log likelihood = 229.08,  aic = -450.15




# a)
# Explanation: Visual inspection of lag=1 differenciation seems to
# favour that model. SARMA coefficients chosen from BIC criterion
# gives a seasonal ARIMA (1,1,0) x (2,0,0)_12



# b) 95% bounds:
# AR-component:
c(airpass.mod$coef[1]-1.96*sqrt(airpass.mod$var[1,1]),airpass.mod$coef[1],
  airpass.mod$coef[1]+1.96*sqrt(airpass.mod$var[1,1]))
#        ar1        ar1        ar1
# -0.4737862 -0.3067691 -0.1397520

# SAR-component 1:
c(airpass.mod$coef[2]-1.96*sqrt(airpass.mod$var[2,2]),airpass.mod$coef[2],
  airpass.mod$coef[2]+1.96*sqrt(airpass.mod$var[2,2]))
#      sar1      sar1      sar1
# 0.3716280 0.5314117 0.6911955

# SAR-component 2:
c(airpass.mod$coef[3]-1.96*sqrt(airpass.mod$var[3,3]),airpass.mod$coef[3],
  airpass.mod$coef[3]+1.96*sqrt(airpass.mod$var[3,3]))
#      sar2      sar2      sar2
# 0.2681824 0.4332741 0.5983658



# c) Plot residuals plus extra diagnostic plots:
tsdiag(airpass.mod)
# Looks ok



# d) Plot predictions

airpass.full=ts(log(airpass.vec),start=1949,frequency=12)
pred=predict(airpass.mod, n.ahead=12)$pred
s.pred=predict(airpass.mod, 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.




# e) Show predictions and prediction errors

cbind(pred, pred-1.96*s.pred, pred+1.96*s.pred)
#Jan 1960 6.036146             5.962680             6.109613
#Feb 1960 5.982345             5.892951             6.071738
#Mar 1960 6.128905             6.022430             6.235381
#Apr 1960 6.098793             5.978618             6.218968
#May 1960 6.148275             6.015535             6.281015
#Jun 1960 6.288723             6.144587             6.432860
#Jul 1960 6.420523             6.265805             6.575241
#Aug 1960 6.443268             6.278653             6.607882
#Sep 1960 6.246455             6.072505             6.420406
#Oct 1960 6.126783             5.943973             6.309592
#Nov 1960 6.000934             5.809675             6.192194
#Dec 1960 6.096765             5.897413             6.296116



# f) Show forecast errors 
len=length(airpass.full)
index=(len-11):len
cbind(airpass.full[index],airpass.full[index]-pred)
#Jan 1960            6.033086               -0.003060221
#Feb 1960            5.968708               -0.013637089
#Mar 1960            6.037871               -0.091034208
#Apr 1960            6.133398                0.034605262
#May 1960            6.156979                0.008703780
#Jun 1960            6.282267               -0.006456715
#Jul 1960            6.432940                0.012416782
#Aug 1960            6.406880               -0.036387932
#Sep 1960            6.230481               -0.015973639
#Oct 1960            6.133398                0.006615526
#Nov 1960            5.966147               -0.034787737
#Dec 1960            6.068426               -0.028339072

# The last real value is inside the forecast confidence interval


plot(airpass.full[index]-pred, type="o")
# Doesn't seem to bad

# Take a look at the RMSE:
sqrt(mean((airpass.full[index]-pred)^2))
#[1] 0.03365016





# Note that the roots for the SAR-coefficients of the model
# lies dangerously close to the unit circle ((0.53+sqrt(0.53^2+4*0.43))/2=0.93). 
# Could try lag=12 differenciation in addition to lag=1 differenciation:

acf(diff(diff(airpass),lag=12),lag=40)
# Seems to have smaller acf than the other differenciations

zero.I1_12.model=arima(airpass,order=c(0,1,0),seasonal=c(0,1,0))
best.I1_12.model=zero.I1_12.model
bicbest.I1_12.model=zero.I1_12.model
n=length(I1_12)
best.I1_12.bic=-2*zero.I1_12.model$loglik+2*log(n)
bicbest.I1_12.p=0
bicbest.I1_12.q=0
bicbest.I1_12.P=0
bicbest.I1_12.Q=0
for(Q in 0:3)
{
  for(P in 0:3)
    {
      for(q in 0:3)
		 {
		   for(p in 0:3)
		     {
		       model<-arima(airpass,order=c(p,1,q),seasonal=c(P,1,Q),method="ML")
		       curr.bic=-2*model$loglik+log(n)*(P+Q+p+q+1+1)
		       if(curr.bic<best.I1_12.bic)
		 		 {
		 		   best.I1_12.bic=curr.bic
		 		   bicbest.I1_12.model=model
		 		   bicbest.I1_12.p=p
		 		   bicbest.I1_12.q=q
		 		   bicbest.I1_12.P=P
		 		   bicbest.I1_12.Q=Q
		 		   
		 		   show(c("BIC:",P,Q,p,q,best.I1_12.bic))
		 		 }
		     }
		 }
    }
}

airpass.mod=bicbest.I1_12.model
airpass.mod
#Coefficients:
#         ar1      ma1     sar1     sar2
#      0.1169  -0.4734  -0.5510  -0.2180
#s.e.  0.3183   0.2919   0.0908   0.0998 

# We get a model with one AR-component, one MA-component
# and two seasonal components: (1,1,1)x(2,1,0)_12. 
# Slighly more complex than what we got earlier.


# c) Plot residuals plus extra diagnostic plots:
tsdiag(airpass.mod)
# Looks quite good



# d) Plot predictions

airpass.full=ts(log(airpass.vec),start=1949,frequency=12)
pred=predict(airpass.mod, n.ahead=12)$pred
s.pred=predict(airpass.mod, 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)
# The predictions looks good

index=(len-11):len
plot(airpass.full[index]-pred, type="o")
# Seems about the same as for the (1,1,0)x(2,0,0)_12 model

# Take a look at the RMSE:
sqrt(mean((airpass.full[index]-pred)^2))
#[1] 0.03538331
# Slightly (but maybe not significantly) worse than the 
# (1,1,0)x(2,0,0)_12 model







# 6.8

# Fit a seasonal + quadratic trend 

t=1:12
ses.mean=vector(length=12)
for(i in t)
  ses.mean[i]=mean(airpass[i+12*(0:10)])
plot(t,ses.mean, type="l")

t=time(airpass)
t1=rep(c(1,rep(0,11)),11)
t2=rep(c(0,1,rep(0,10)),11)
t3=rep(c(rep(0,2),1,rep(0,9)),11)
t4=rep(c(rep(0,3),1,rep(0,8)),11)
t5=rep(c(rep(0,4),1,rep(0,7)),11)
t6=rep(c(rep(0,5),1,rep(0,6)),11)
t7=rep(c(rep(0,6),1,rep(0,5)),11)
t8=rep(c(rep(0,7),1,rep(0,4)),11)
t9=rep(c(rep(0,8),1,rep(0,3)),11)
t10=rep(c(rep(0,9),1,rep(0,2)),11)
t11=rep(c(rep(0,10),1,rep(0,1)),11)
t12=rep(c(rep(0,11),1),11)

airpass.lin<-lm(airpass~t1+t2+t3+t4+t5+t6+t7+t8+t9+t10+t11+t12+I(t-1960)+I((t-1960)*(t-1960))-1,x=T)
remainder<-airpass.lin$res
acf(remainder)
pacf(remainder)

# Seems like AR 1 or AR 2

# Find best ARMA model, using BIC:
zero.rem.model=arima(remainder)
best.rem.model=zero.rem.model
bicbest.rem.model=zero.rem.model
n=length(remainder)
best.rem.bic=-2*zero.rem.model$loglik+log(n)
bicbest.rem.p=0
bicbest.rem.q=0
for(q in 0:3)
  {
    for(p in 0:3)
      {
		 model<-arima(remainder,order=c(p,0,q),method="ML")
		 curr.bic=-2*model$loglik+log(n)*(p+q+1)
		 if(curr.bic<best.rem.bic)
		   {
		     best.rem.bic=curr.bic
		     bicbest.rem.model=model
		     bicbest.rem.p=p
		     bicbest.rem.q=q
		     
		     show(c("BIC:",p,q,best.rem.bic))
		   }
      }
  }

bicbest.rem.model
#Coefficients:
#         ar1  intercept
#      0.6995     0.0012
#s.e.  0.0618     0.0093
# 
#sigma^2 estimated as 0.001069:  log likelihood = 263.85,  aic = -521.71

# Best model seems to be AR(1) fro the remainder


# a) Chose to let each month have it's own regression coefficient
#    The remaining signal seems to be best described by a 
#    AR(1) model.




# b) Bounds for ARMA-coefficients:

# For the given remainder, it's:
c(bicbest.rem.model$coef[1]-1.96*sqrt(bicbest.rem.model$var[1,1]),
  bicbest.rem.model$coef[1],
  bicbest.rem.model$coef[1]+1.96*sqrt(bicbest.rem.model$var[1,1]))
#       ar1       ar1       ar1
# 0.5783635 0.6994532 0.8205430

# This does however not take into account the uncertainty involved
# in making the remainder


# Could try a parametric bootstrap algorithm for that:
# but this is much work!




# c) Plot residuals
tsdiag(bicbest.rem.model)
# Looks ok


# d) Forecast graph

sigma=summary(airpass.lin)$sigma
x=airpass.lin$x
v=solve((t(x)%*%x))*sigma^2

airpass.full=ts(log(airpass.vec),start=1949,frequency=12)
pred=predict(bicbest.rem.model, n.ahead=12)$pred
s.pred=predict(bicbest.rem.model, n.ahead=12, se.fit=T)$se
for(t in 1:12)
{
  x0=c(t==1,t==2,t==3,t==4,t==5,t==6,
      t==7,t==8,t==9,t==10,t==11,t==12, 
      (t-1)/12, ((t-1)/12)^2)

  pred[t]=pred[t]+t(airpass.lin$coef)%*%x0
  s.pred[t]=s.pred[t]+t(x0)%*%v%*%x0
}

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:12]+1.96*s.pred[1:12]), type="l", lty=3)
lines(t, c(airpass[length(airpass)], pred[1:12]-1.96*s.pred[1:12]), type="l", lty=4)

# The predictions seems to be okay here to. The confidence limits
# seems to be about the same width.


# e) Show forecast uncertainty
cbind(pred, pred[1:12]-1.96*s.pred[1:12], pred[1:12]+1.96*s.pred[1:12])
#133 6.027210                         5.953019                         6.101401
#134 6.014607                         5.924471                         6.104742
#135 6.156623                         6.049387                         6.263860
#136 6.120509                         5.999553                         6.241464
#137 6.122768                         5.989226                         6.256309
#138 6.251836                         6.106877                         6.396796
#139 6.358890                         6.203327                         6.514453
#140 6.358585                         6.193102                         6.524067
#141 6.224378                         6.049535                         6.399220
#142 6.090077                         5.906350                         6.273804
#143 5.956127                         5.763925                         6.148329
#144 6.078699                         5.878378                         6.279019


# f) Show forecast errors:
len=length(airpass.full)
index=(len-11):len
cbind(airpass.full[index],airpass.full[index]-pred)
#133            6.033086                0.005876460
#134            5.968708               -0.045899328
#135            6.037871               -0.118752288
#136            6.133398                0.012889464
#137            6.156979                0.034211343
#138            6.282267                0.030430556
#139            6.432940                0.074050249
#140            6.406880                0.048295188
#141            6.230481                0.006103747
#142            6.133398                0.043320889
#143            5.966147                0.010019808
#144            6.068426               -0.010273084

# The actual values sems to be inside the confidence intervals again.
# and the conficence intervals seems to have the same width


# Take a look at the RMSE:
sqrt(mean((airpass.full[index]-pred)^2))
# 0.04868411



# The root mean square error seems to be somewhat larger here than
# for the seasonal ARIMA model.