newton_fit_ARMA_02 = function(w_t,theta10=0,theta20=1/3,verbose=F){
#
# Given a time series of differences in w_t and an initial guess at the MA(2)
# parameters theta0, theta1 this routine iterates using Newton's method to
# find convergent estimates. 
#
# The initial guess theta10,theta20 is taken to be the centroid of the invertable region.
#
# 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('theta_to_lambda_ARMA_02.R')
source('uncond_sum_of_squares_ARMA_02.R')
source('numerically_evaluate_xt_ARMA_02.R')

thetas0 = c( Inf, Inf ) # 0 => old guess 
thetas1 = as.matrix(c( theta10, theta20 )) # 1 => new guess 

# Do each newton step:
# 
while( sqrt( sum( (thetas1 - thetas0)^2 ) ) > 1.e-3 ){
  if( verbose ){ print(t(thetas1)) }
  thetas0 = thetas1 # get the most recent solution 

  # Get the two values of lambda from the values of theta
  tmp = theta_to_lambda_ARMA_02(thetas0[1],thetas0[2]); lambda00 = tmp[1]; lambda10 = tmp[2]; 

  res = uncond_sum_of_squares_ARMA_02(lambda00,lambda10,w_t)
  a_t = res[[3]]
  x_t = numerically_evaluate_xt_ARMA_02(thetas0[1],thetas0[2],a_t,w_t)
  
  # Solve the least square problem: 
  ### delta_theta = solve( t(x_t) %*% x_t, t(x_t) %*% as.matrix(a_t) )
  res = lsfit( x_t, a_t, intercept=FALSE )
  delta_theta = res$coeff

  thetas1 = thetas0 + delta_theta
}
# respect invertable boundaries of a MA(2) process:
#
if( abs(thetas1[1])>2 ){ thetas1[1] = 2*sign(thetas1[1]) }
if( thetas1[2] < 0 ){ thetas1[2] = 0. }
if( thetas1[1] > 0 & thetas1[2] > 1 - 0.5*thetas1[1] ){ thetas1[2] = 1 - 0.5*thetas1[1] }
if( thetas1[1] > 0 & thetas1[2] > 1 + 0.5*thetas1[1] ){ thetas1[2] = 1 + 0.5*thetas1[1] }

# 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( thetas1, Var_Cov_Beta ) ) 

}