#
# 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('est_theta1_MA1.R')
source('wt_std_error_MA1.R')

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

# For z: 
z_acf = c(0.95,0.90,0.86,0.80,0.76,0.71,0.66,0.61,0.57,0.52)
z_pacf = c(0.95,0.02,-0.02,-0.06,0.00,-0.02,-0.05,-0.01,-0.02,-0.04)

# For nabla z:
dz_acf = c(0.01,0.09,0.17,0.02,-0.04,-0.01,-0.24,-0.18,-0.03,-0.08)
dz_pacf = c(0.01,0.09,0.17,0.01,-0.08,-0.04,-0.25,-0.18,0.02,0.04)

n = 56

# For z:
# 
if( save_Plots ){
  fn=paste(paste(DN,"prob_4_z_acf",sep=""),".eps",sep="")
  postscript(fn, onefile=FALSE, horizontal=FALSE)
}
sdt = sqrt( ( 1 + 2 * sum( z_acf[1:3]^2 ) ) / n ) # an approximation
ylim = c( min( c(z_acf,-sdt) ), max( c(z_acf,+sdt) ) )
plot( z_acf, ylim = ylim, main="ACF(z)" )
abline(h= 2*sdt,col=2)
abline(h=-2*sdt,col=2)
if( save_Plots ){ dev.off() }

if( save_Plots ){
  fn=paste(paste(DN,"prob_4_z_pacf",sep=""),".eps",sep="")
  postscript(fn, onefile=FALSE, horizontal=FALSE)
}
plot( z_pacf, main="PACF(z)" )
abline(h= 2./sqrt(n),col=1)
abline(h=-2./sqrt(n),col=1)
if( save_Plots ){ dev.off() }

# For \nabla z:
# 
if( save_Plots ){
  fn=paste(paste(DN,"prob_4_delta_z_acf",sep=""),".eps",sep="")
  postscript(fn, onefile=FALSE, horizontal=FALSE)
}
sdt = sqrt( ( 1 + 2 * sum( dz_acf[1:3]^2 ) ) / n ) # an approximation
ylim = c( min( c(dz_acf,-sdt) ), max( c(dz_acf,+sdt) ) )
plot( dz_acf, ylim = ylim, main="ACF(Delta(z))" )
abline(h= 2*sdt,col=2)
abline(h=-2*sdt,col=2)
if( save_Plots ){ dev.off() }

if( save_Plots ){
  fn=paste(paste(DN,"prob_4_delta_z_pacf",sep=""),".eps",sep="")
  postscript(fn, onefile=FALSE, horizontal=FALSE)
}
plot( dz_pacf, main="PACF(z)" )
abline(h= 2./sqrt(n),col=1)
abline(h=-2./sqrt(n),col=1)
if( save_Plots ){ dev.off() }

# Based on these plots we consider a ARIMA(0,1,1) model ...
#
r1 = dz_acf[1]
r2 = dz_acf[2]
theta1 = est_theta1_MA1(r1)

zbar = 0.66 # for w_t 
c0 = sz2 = 0.7931 

# Lets see if the mean is significant (I would guess yes)
#
sigma_w_bar = wt_std_error_MA1(n,c0,r1,r2)

theta0 = zbar
sigmaa2 = c0 / ( 1 + theta1^2 )