numerically_evaluate_xt_ARMA_02 = function(theta10,theta20,a_t,w_t){
#
# Given a time series of differences in w_t and the initial time 
# series of white noise in a_t and the initial MA(2) parameters theta1,theta2
# this routine numerically computes an estimate of the (negative) derivative of 
# a_t with resepect to theta1 and
# a_t with resepect to theta2
#
# 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')    

param_step = 0.02

if( theta10 != 0 ){ 
  delta_theta1 = param_step * theta10 # take a step in theta1 ... call S() using arguments in lambda
}else{
  delta_theta1 = param_step 
}
new_theta1 = theta10 + delta_theta1
tmp = theta_to_lambda_ARMA_02(new_theta1,theta20); lambda00 = tmp[1]; lambda10 = tmp[2] 
res = uncond_sum_of_squares_ARMA_02(lambda00,lambda10,w_t)
dx_dTheta1 = ( a_t - res[[3]] )/delta_theta1

if( theta20 != 0 ){ 
  delta_theta2 = param_step * theta20 # take a step in theta2 ... call S() using arguments in lambda
}else{
  delta_theta2 = param_step 
}
new_theta2 = theta20 + delta_theta2
tmp = theta_to_lambda_ARMA_02(theta10,new_theta2); lambda00 = tmp[1]; lambda10 = tmp[2] 
res = uncond_sum_of_squares_ARMA_02(lambda00,lambda10,w_t)
dx_dTheta2 = ( a_t - res[[3]] )/delta_theta2

out = as.matrix( cbind( dx_dTheta1, dx_dTheta2 ) ) # returns a matrix with derivatives along each column
colnames(out) = NULL
return(out) 
    
}