# Part of the code is from Shumway and Stoffer (2011) page 100
# -----------------------------------------------------------------------------------
z = c(1,-1.5,.75)             # coefficients of the polynomial
(a = polyroot(z)[1])         # print one root: 1+0.57735i = 1 + i/sqrt(3)
arg = Arg(a)/(2*pi)          # arg in cycles/pt
1/arg                              # = 12, the pseudo period

phis = c(1.5,-0.75)
#phis = c(0.8,0.1)
#phis = c(-0.8,0.0)
#phis = c(0.449,0.449)


pdf(file="ar2graphs.pdf",width=10,height=6)
# Simulating AR(2) data and plotting the autocorrelation function
set.seed(90210)
ar2 = arima.sim(list(order=c(2,0,0), ar=phis), n = 144)
ACF = ARMAacf(ar=phis, ma=0, 50)

par(mfrow=c(1,2))
plot(1:144/12, ar2, type="l", xlab="Time (one unit = 12 points)")
abline(v=0:12, lty="dotted")
plot(ACF, type="h", xlab="lag")
abline(h=0)

# Studying the surface of the AR(2) likelihood function
like = function(phi1,phi2){
  sum(dnorm(y[3:n],phi1*y[2:(n-1)]+phi2*y[1:(n-2)],1,log=TRUE))
}

phi1 = seq(phis[1]-0.5,phis[1]+0.5,length=100)
phi2 = seq(phis[2]-0.5,phis[2]+0.5,length=100)
likes = matrix(0,100,100)
for (i in 1:100)
for (j in 1:100)
  likes[i,j] = like(phi1[i],phi2[j])

par(mfrow=c(1,1))
contour(phi1,phi2,exp(likes),xlim=c(-2,2),ylim=c(-1,1),drawlabels=FALSE)
abline(h=0,lty=3)
abline(v=0,lty=3)
segments(0,1,2,-1,lty=2)
segments(0,1,-2,-1,lty=2)
segments(-2,-1,2,-1,lty=2)
points(phis[1],phis[2],col=2,pch=16)
points(reg$coef[1],reg$coef[2],col=3,pch=16)
legend("topright",legend=c("True","MLE"),col=2:3,pch=16)

# Maximum likelihood estimation of the AR(2) model
y = ar2
n = length(y)
reg = lm(y[3:n]~y[2:(n-1)]+y[1:(n-2)]-1)
summary(reg)

plot(reg$res,type="l",main="Residual analysis",ylab="Residuals",xlab="Time")



dev.off()