# Order selection and estimation of coefficients for ARMA models
# Main functions: ar.yw, ar.mle,ar.burg,ar.ols,spec.ar
#                 arima, arima.sim, tsdiag
                  

# 1. Lynx: Comparing the 4 methods Yule-Walker, MLE, Burg and Least Squares
plot(log(lynx))

par(mfrow=c(1,2))
acf(log(lynx))
acf(log(lynx),type="partial")

fit.1 <- ar.yw(log(lynx))
fit.1
str(fit.1)

fit.2 <- ar.mle(log(lynx))
fit.2

fit.3 <- ar.burg(log(lynx))
fit.3

fit.4 <- ar.burg(log(lynx),order=11)
fit.4

fit.5 <- ar.ols(log(lynx))
fit.5
fit.5$var.pred <- drop(fit.5$var.pred) ## work around R buglet

fit.6 <- ar.ols(log(lynx),order=11)
fit.6

fit.7 <- ar.ols(log(lynx),order=12)


plot(fit.1$ar,ylab="",xlab="lag",main="estimated coefficients")
points(fit.2$ar,col=2)
points(fit.4$ar,col=3)
points(fit.6$ar,col=4)
legend(x="topright",c("Yule-Walker","MLE","Burg","OLS"),col=1:4,pch=rep(1,4))

plot(0:20,fit.1$aic,type="b",ylab="",xlab="order",main="AIC")
lines(0:12,fit.2$aic,type="b",col=2)
lines(0:20,fit.3$aic,type="b",col=3)
lines(0:20,fit.5$aic,type="b",col=4)
legend(x="topright",c("Yule-Walker","MLE","Burg","OLS"),col=1:4,pch=rep(1,4))
abline(h=0)

# standard errors
round(sqrt(diag(fit.1$asy.var.coef)),3)
round(sqrt(diag(fit.2$asy.var.coef)),3)
round(sqrt(diag(fit.4$asy.var.coef)),3)
round(fit.6$asy.se.coef$ar,3)

# spectra of fitted models
spec.ar(fit.1)
spec.ar(fit.2,add=TRUE,col=2)
spec.ar(fit.3,add=TRUE,col=3)
spec.ar(fit.5,add=TRUE,col=4) 
legend(x="topright",c("Yule-Walker","MLE","Burg","OLS"),col=1:4,pch=rep(1,4))


# for comparison: Spectrum of fitted AR(2)
fit.7 <- ar.mle(log(lynx),order.max=2)
spec.ar(fit.7)

# analysis of residuals

resid <- fit.1$resid[-(1:11)] # acf does not take NA's
par(mfrow=c(1,2))
plot(resid,type="l")
plot(resid^2,type="l")
acf(resid)
acf(resid,type="partial")
acf(resid^2)
acf(resid^2,type="partial") #lag 2: borderline case
qqnorm(resid)
plot(log(lynx))
plot(arima.sim(n = 114, list(ar=fit.3$ar),sd=sqrt(fit.3$var.pred))+
     mean(log(lynx)),ylab="simulated series")
# using what is provided in R
fit.8 <- arima(log(lynx),order=c(11,0,0)) #need to fit with a different function
tsdiag(fit.8)

# 2. Simulations: Normal and t-distributed innovations
ar.coef1 <- matrix(0,nrow=50,ncol=2)
sd.1 <- c(0,0)
for (i in (1:50)) {
  x <- arima.sim(n = 100, list(ar=c(1.38, -0.74)),
               rand.gen = function(n, ...) rt(n, df = 5))
  ar.fit <- ar.burg(x,order.max=2)
  ar.coef1[i,] <- (ar.fit$ar -c(1.38, -0.74))
  sd.1 <- sd.1 + (sqrt(diag(ar.fit$asy.var.coef))-sd.1)/i
}
sd.1 #mean of standard errors

ar.coef2 <- matrix(0,nrow=50,ncol=2)
sd.2 <- c(0,0)
for (i in (1:50)) {
  x <- arima.sim(n = 100, list(ar=c(1.38, -0.74)),
               rand.gen = function(n, ...) rnorm(n))
  ar.fit <- ar.burg(x,order.max=2)
  ar.coef2[i,] <- (ar.fit$ar -c(1.38, -0.74))#/sqrt(diag(ar.fit$asy.var.coef))
  sd.2 <- sd.2 + (sqrt(diag(ar.fit$asy.var.coef))-sd.2)/i
}
sd.2 #mean of standard errors

boxplot(ar.coef1,main=("estimation error, t innovations"))
abline(h=2*sd.1[1],col=2)
abline(h=-2*sd.1[1],col=2)
boxplot(ar.coef2,main=("estimation error, normal innovations"))
abline(h=2*sd.2[1],col=2)
abline(h=-2*sd.2[1],col=2)

# Do the same thing with Yule Walker instead of Burg. Notice slightly
# larger bias and standard deviation for Yule Walker


# 3. Airline data
plot(AirPassengers)
plot(log(AirPassengers))
acf(log(AirPassengers))
acf(log(AirPassengers),type="partial")
acf(diff(log(AirPassengers)))
acf(diff(log(AirPassengers)),type="partial")
acf(diff(diff(log(AirPassengers)),lag=12))
acf(diff(diff(log(AirPassengers)),lag=12),type="partial")

fit.1 <- arima(log(AirPassengers),order=c(0,1,1),seasonal=list(order=c(0,1,1),
               period=12))
tsdiag(fit.1)
plot(fit.1$resid)
qqnorm(fit.1$resid)
acf(fit.1$resid)
acf(fit.1$resid,type="partial")

pred.1 <- predict(fit.1,n.ahead=12)
pred.1
airp <- ts(c(log(AirPassengers),pred.1$pred),start=c(1949,1),freq=12)
par(mfrow=c(1,1))
plot(airp)
lines(pred.1$pred,col=2)
lines(pred.1$pred+2*pred.1$se,col=3)
lines(pred.1$pred-2*pred.1$se,col=3)