#
# Written by:
# -- 
# John L. Weatherwax                2009-04-21
# 
# email: wax@alum.mit.edu
# 
# Please send comments and especially bug reports to the
# above email address.
# 
#-----

source('wt_std_error_AR2.R')

save_Plots = F
DN="../../WriteUp/Graphics/Chapter6/"

n = 60 

r = c(0.93,0.80,0.65,0.49,0.32,0.16,0.03,-0.09,-0.16,-0.22,-0.25,-0.25,-0.21,-0.12,-0.01,0.10)
if( save_Plots ){
  fn=paste(paste(DN,"prob_3_acf",sep=""),".eps",sep="")
  postscript(fn, onefile=FALSE, horizontal=FALSE)
}
sdt = sqrt( ( 1 + 2 * sum( r[1:3]^2 ) ) / n ) # an approximation
ylim = c( min( c(r,-sdt) ), max( c(r,+sdt) ) )
plot( r, ylim = ylim, main="ACF" )
abline(h= 2*sdt,col=2)
abline(h=-2*sdt,col=2)
if( save_Plots ){ dev.off() }


p = c(0.93,-0.41,-0.14,-0.11,-0.07,-0.10,0.05,-0.07,0.12,-0.14,0.03,0.09,0.19,0.20,0.03,-0.11)
if( save_Plots ){
  fn=paste(paste(DN,"prob_3_pacf",sep=""),".eps",sep="")
  postscript(fn, onefile=FALSE, horizontal=FALSE)
}
plot( p, main="PACF" )
abline(h= 2./sqrt(n),col=1)
abline(h=-2./sqrt(n),col=1)
if( save_Plots ){ dev.off() }


# Assume an AR(2) model and obtain preliminary estimates of the coefficients:
#
r1 = r[1]
r2 = r[2]

phi1 = r1 * ( 1 - r2 )/( 1 - r1^2 ) 
phi2 = ( r2 - r1^2 )/( 1 - r1^2 )

# Test if the mean of z is significant:
#
zbar = 2.56
c0 = 0.01681

sigma_w_bar = wt_std_error_AR2(n,c0,r1,r2)

# Yes it is significant
theta0 = ( 1. - phi1 - phi2 ) * zbar

sigmaa2 = ( 1 - phi1*r1 - phi2*r2 ) * c0