# K. Ensor (ensor@rice.edu)
# Rice University
#
# Example of different estimating procedures for
# an AR(2) model - easily generalized
#


# Plotting function used later

prettyplot<-function(nrow,ncol){
	par(mfrow=c(nrow,ncol),omi=c(.5,1,.5,1))
}


x<-arima.sim(model=list(ar=c(0,.9)),n=200)
kdescripb(x,"Simulated AR(2) Realization")

#Parameter estimation

ywe<-ar(x,order=2,method="yule-walker")
burg<-ar(x,order=2,method="burg")
leastsq<-lm(x[3:200]~x[2:199]+x[1:198]-1)
mle<-arima.mle((x-mean(x)),model=list(ar=c(.2,.6)))
mlediag<-arima.diag(mle, resid=T)

prettyplot(2,2)
tsplot(ywe$resid)
tsplot(burg$resid)
tsplot(leastsq$residuals)
tsplot(mlediag$resid)

# Compare variance estimates - true value is 1
ywe$var.pred
burg$var.pred
var(resid(leastsq),SumSquares=T)/(198-2)
mle$sigma2

leastsq$resid<-c(NA,NA,residuals(leastsq))
resid<-data.frame(ywe$resid,burg$resid,leastsq$resid,mlediag$resid)
prettyplot(1,1)
tsplot(resid)

prettyplot(2,3)
for (i in 1:3) {
	for (j in (i+1):4)
		plot(resid[,i],resid[,j])
}

prettyplot(1,1)
plot(leastsq$resid,ywe$resid)