uncond_sum_of_squares_ARMA_20 = function(phi1,phi2,w_t,number_of_sweeps=2){
#
# Given a time series of differences in w_t this function computes the unconditional sum of squares
# for an AR(2) MA(0) model.
#
# 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.
# 
#-----

# For our model we only need to pad with this many zeros
Q = 10

# expand w_t beyond its original limits and initialized e_t and a_t
w_t = c( rep(0,Q), w_t, rep(0,Q) )
e_t = rep(0,length(w_t)) 
a_t = rep(0,length(w_t)) 

np = length(w_t)

for( si in seq(1,number_of_sweeps) ){
  
 # Perform a backwards sweep to compute [e_t]
 #
 for( ii in seq(np-Q-1,Q+1,by=-1) ){
   e_t[ii] = w_t[ii] - phi1 * w_t[ii+1] - phi2 * w_t[ii+2] 
 }

 # Perform a backwards sweep to compute [w_t] for negative t
 #
 for( ii in seq(Q,1,by=-1) ){
   w_t[ii] = e_t[ii] + phi1 * w_t[ii+1] + phi2 * w_t[ii+2] 
 }

 # Perform a forwards sweep to compute [a_t] for all t
 #
 for( ii in seq(3,np-Q,by=+1) ){
   a_t[ii] = w_t[ii] - phi1 * w_t[ii-1] - phi2 * w_t[ii-2] 
 }

 # Perform a forwards sweep to compute [w_t] for t > n
 #
 for( ii in seq(np-Q,np-1,by=+1) ){
   w_t[ii] = a_t[ii] + phi1 * w_t[ii-1] + phi2 * w_t[ii-2] 
 }
}

return( list( sum(a_t^2), e_t, a_t ) ) 

}