#aar2 <- function(x, f=-1)
  
 x <- fig89.dat
 f <- -1   #0.6 

   #Fitting an AAR(2) model with p = 2 below.  The program can be
   #adapted to fit other AAR(p) model
   # X is a time series  and f is the
   #parameter that is related to the bandwidth.  f < 0 means automatic
   #selection.  For f in (0,1), the fraction of number of
   #observations in a neighborhood. This is the smoothing parameter
   #for a loess fit.  For p other than 2, 
   #the following program needs to be revised.
   #
########**************************************################
########DUE TO SOME WEIRD PROBLEMS in Splus   ###############
########The following needed to be excuted by batch##########
########The result is "res.term"              ###############
  if(f >= 1) 
   {
    stop("Warning:  f is out of range (0,1). Automatic bandwidth is used")
   }
  n <- length(x)

  Q = 3; m = floor(0.1*n); p = 2
  #### creating regression structure
  Y <- x[(p+1):n]
  X1 <- x[p:(n-1)]
  X2 <- x[(p-1):(n-2)]
				#add more terms if p > 2 is large
				

   ################## Selecting optimal bandwidth ###########
 if(f < 0)  
 {
 MSE = rep(0,4)                 #initial values

 for(i in 1:4)
   {
    f = i*0.2
    for(q in 1:Q)
    {
     print(c(i,q))
     j = n-q*m-p;               #length of the series
     g = (n/j)^0.2*f            #adjusted bandwidth
     series <- data.frame(y= Y[1:j],x1=X1[1:j],x2=X2[1:j])
     res.gam   <- gam(y ~ lo(x1,g) + lo(x2,g),data=series)
     res.term  <- predict(res.gam, type="term")
     x1=X1[(j+1):(j+m)]; 
     x2=X2[(j+1):(j+m)];
     pred = approx(X1[1:j],res.term[,1],x1,rule=3)$y +
	    approx(X2[1:j],res.term[,2],x2,rule=3)$y + mean(Y[1:j])
     MSE[i] = MSE[i] + sum((Y[(j+1):(j+m)] - pred)^2)
    }
   }

  f=0.2* (1:4)[MSE==min(MSE)]                     #optimal bandiwdth 
 }

 print(paste("bandwidth = ", f))
 series <- data.frame(y= Y,x1=X1,x2=X2)
 res.gam   <- gam(y ~ lo(x1,f) + lo(x2,f),data=series)
 resid     <- res.gam$resid
 res.term  <- predict(res.gam, type="term")
 res.term  <- cbind(res.term,resid)              #last column is residuals



 rm(MSE,f, res.gam, resid, x2, Q, m, X1, n, series, X2, g, p, x, Y, j, pred, x1)
 print("result is stored at `res.term'") 
 print("First two-column are estimated functions")
 print("The last column are residuals")
 fig89.res <- res.term