#
# 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.
# 
#-----

savePlots = F

zt = c( 166,172,172,169,164,168,171,167,168,172 ) # z_91, z_92, to z_{100} 

plot( 91:100, zt, xlab="index", ylab="z_t", main="z_t", xlim=c(91,112) )

# past a_t computed from the given z_t using:
# 
zt_hat = c() # will hold \hat{z}_{91}(1), \hat{z}_{92}(1), ... \hat{z}_{99}(1), \hat{z}_{100}(1)
at_hat = c(0,0) # will hold a_hat(90), a_hat(91), ... a_hat(99), a_hat(100)
for( ti in 91:100 ){
  ztiprime = ti-90 # maps to the observed measurment z_t
  atiprime = ti-90+1 # maps to the estimated at_hat 

  z_value = zt[ztiprime] - 1.1 * at_hat[atiprime] + 0.28 * at_hat[atiprime-1]
  zt_hat = c(zt_hat,z_value)
  
  a_value = zt[ztiprime+1] - zt_hat[ztiprime]
  at_hat = c(at_hat,a_value)
  
}

# print the results thus far: 
zt_hat
at_hat

#lines( 92:101, zt_hat ); # plot the onestep ahead predictions

# append to zt_hat values for zhat_{100}(2), and zhat_{100}(l)
#
z_100_1 = zt_hat[length(zt_hat)] # the value for \hat{z}_{100}(1) 
z_100_2 = 168.8347 # the value for \hat{z}_{100}(2) 
zhat_l = c(z_100_1,z_100_2)
for( l in 1:10){ # append \hat{z}_{100}(l) for l = 3, 4, 5, ... 12 (or 10 more numbers)
  zhat_l = c(zhat_l,z_value)
}
lines( 101:112, zhat_l, col="red" ); # plot the predictions 

# compute the confidence intervals:
#
psi_1 = -0.1
psi_j = 0.18

sigma_a = sqrt(1.103)

sigma2_l = c( sigma_a^2 )
for( l in 2:12 ){
  v = sigma_a^2 * ( 1 + psi_1^2 + (l-2) * psi_j^2 )
  sigma2_l = c(sigma2_l, v)
}

p_value = 0.8
c = qnorm( 1 - (1-p_value)/2 )
lower_limits = zhat_l - c * sqrt( sigma2_l )
upper_limits = zhat_l + c * sqrt( sigma2_l )

# plot the confidence intervals: 
#
lines( 101:112, lower_limits, col="green" )
lines( 101:112, upper_limits, col="green" ) 

if( savePlots ){
  postscript("../../WriteUp/Graphics/Chapter5/prob_2.eps", onefile=FALSE, horizontal=FALSE)
  plot( 91:100, zt, xlab="index", ylab="z_t", main="z_t", xlim=c(91,112) )
  lines( 101:112, zhat_l, col="red" ); # plot the predictions 
  lines( 101:112, lower_limits, col="green" )
  lines( 101:112, upper_limits, col="green" )   
  dev.off() 
}