newton_fit_ARMA_01 = function(w_t,theta10=0.,verbose=F){
#
# Given a time series of differences in w_t and an initial guess at the MA(1) parameter theta0
#
# 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('uncond_sum_of_squares_ARMA_01.R')
source('numerically_evaluate_xt_ARMA_01.R')

theta11 = theta10
theta10 = Inf 

# Do each newton steps:
# 
while( abs(theta11 - theta10) > 1.e-6 ){
  if( verbose ){ print(theta11) }
  theta10 = theta11

  res = uncond_sum_of_squares_ARMA_01(theta10,w_t)
  a_t = res[[3]]
  x_t = numerically_evaluate_xt_ARMA_01(theta10,a_t,w_t)

  ### delta_theta = sum( a_t * x_t )/sum( x_t^2 )
  res = lsfit( x_t, a_t, intercept=FALSE )
  delta_theta = res$coeff

  theta11 = theta10 + delta_theta 
}
if( abs(theta11)>1.0 ){ theta11 = sign(theta11) } # respect invertable boundaries of a MA(1) process

# Get an approximation of the variance-covariance matrix of the betas:
# 
S_hat = res[[1]]
sigmaa2 = S_hat/length(w_t) # an approximation to \sigma_a^2
X = as.matrix(x_t)
Var_Cov_Beta = solve( t(X) %*% X ) * sigmaa2

return( list( theta11, Var_Cov_Beta ) ) 

}