#
# 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('cond_sum_of_squares_ARMA_01.R')
source('uncond_sum_of_squares_ARMA_01.R')

thetas = seq(-0.5,+0.5,by=0.1)

z_t = scan("../../Data/series_b.dat",strip.white=T)
w_t = diff(z_t)
#w_t = w_t[1:9] # for debugging the example worked in the book 

# The *conditional* sum of squares function S_*(theta)
# 
all_s_stars = c()
for( ti in 1:length(thetas) ){
  ss = cond_sum_of_squares_ARMA_01(thetas[ti],w_t)
  all_s_stars = c(all_s_stars,ss)
}
print("The conditional sum of squares function")
print(all_s_stars)

# The *unconditional* sum of squares function S_(theta)
#
all_s = c()
for( ti in 1:length(thetas) ){
  ss = uncond_sum_of_squares_ARMA_01(thetas[ti],w_t)
  all_s = c(all_s,ss[[1]])
}
print("The unconditional sum of squares function")
print(all_s)


# Sampled over a much finer grid:
#
thetas = seq( -0.9, +0.9, by=0.01 )
all_s = c()
for( ti in 1:length(thetas) ){
  ss = uncond_sum_of_squares_ARMA_01(thetas[ti],w_t)
  all_s = c(all_s,ss[[1]])
}

spot = which.min( all_s )
print(sprintf("spot= %d; ss_value= %10.6f; theta= %10.6f",spot,all_s[spot],thetas[spot]))

#postscript("../../WriteUp/Graphics/Chapter7/dup_fig_7_1.eps", onefile=FALSE, horizontal=FALSE)
plot( thetas, all_s, type="l" )
abline( v = thetas[spot], col="red" )
#dev.off()

# Compute the standard error in theta using the techniques from the book:
#
delta_s  = diff( all_s )
delta2_s = diff( all_s, differences=2 ) 

print( all_s[spot] )
print(delta_s[spot])
print(delta2_s[spot])

dtheta = thetas[2] - thetas[1]
Sij = delta2_s[spot]/(dtheta^2) # this is the approximation to Sij

n = length(w_t)

sigmaa2 = all_s[spot]/n # this is the approximation to \sigma_a^2
k = 1 # we are estimating only one parameter theta
q = qchisq( 0.95, k ) # parameter determining confidence interval

# the right-hand-side of (theta - \hat{theta})^2 = 2 sigmaa^2 chisq() / Sij is 
# 
(2 * sigmaa2 * q)/Sij