Nonlinear estimation algorithms are required for obtaining estimates of the parameters of a regression model with innovations that have an ARMA structure. Three estimation methods used by the ARIMA procedure in SAS/ETS software are implemented in the following SAS/IML program. The implemented algorithms are slightly different from those used by PROC ARIMA, but the results should be similar. This example uses the ARMALIK, PRODUCT, and RATIO functions to perform the estimation. Note the interactive nature of this example, illustrating how you can adjust the estimates when they venture outside the stationary or invertible regions.

 /*-------------------------------------------------------------*/
 /*----  Grunfeld's Investment Models Fit with ARMA Errors  ----*/
 /*-------------------------------------------------------------*/
 data grunfeld;
   input year gei gef gec wi wf wc;
   label gei='gross investment ge'
         gec='capital stock lagged ge'
         gef='value of outstanding shares ge lagged'
         wi ='gross investment w'
         wc ='capital stock lagged w'
         wf ='value of outstanding shares lagged w';
 /*---  GE STANDS FOR GENERAL ELECTRIC AND W FOR WESTINGHOUSE  ---*/
 datalines;
 1935     33.1      1170.6    97.8      12.93     191.5     1.8
 1936     45.0      2015.8    104.4     25.90     516.0     .8
 1937     77.2      2803.3    118.0     35.05     729.0     7.4
 1938     44.6      2039.7    156.2     22.89     560.4     18.1
 1939     48.1      2256.2    172.6     18.84     519.9     23.5
 1940     74.4      2132.2    186.6     28.57     628.5     26.5
 1941     113.0     1834.1    220.9     48.51     537.1     36.2
 1942     91.9      1588.0    287.8     43.34     561.2     60.8
 1943     61.3      1749.4    319.9     37.02     617.2     84.4
 1944     56.8      1687.2    321.3     37.81     626.7     91.2
 1945     93.6      2007.7    319.6     39.27     737.2     92.4
 1946     159.9     2208.3    346.0     53.46     760.5     86.0
 1947     147.2     1656.7    456.4     55.56     581.4     111.1
 1948     146.3     1604.4    543.4     49.56     662.3     130.6
 1949     98.3      1431.8    618.3     32.04     583.8     141.8
 1950     93.5      1610.5    647.4     32.24     635.2     136.7
 1951     135.2     1819.4    671.3     54.38     723.8     129.7
 1952     157.3     2079.7    726.1     71.78     864.1     145.5
 1953     179.5     2371.6    800.3     90.08     1193.5    174.8
 1954     189.6     2759.9    888.9     68.60     1188.9    213.5
 ;
 run;

proc iml;
/*-----------------------------------------------------------*/
/*  Estimation for regression model with ARMA errors         */
/*  The ARMAREG module uses the following global parameters: */
/*  x      - matrix of predictors.                           */
/*  y      - response vector.                                */
/*  iphi   - defines indices of nonzero AR parameters,       */
/*           omit the index 0 which corresponds to the zero  */
/*           order constant one.                             */
/*  itheta - defines indices of nonzero MA parameters,       */
/*           omit the index 0 which corresponds to the zero  */
/*           order constant one.                             */
/*  ml     - estimation option: -1 if Conditional Least      */
/*           Squares, 1 if Maximum Likelihood, otherwise     */
/*           Unconditional Least Squares.                    */
/*  delta  - step change in parameters (default 0.005).      */
/*  par    - initial values of parms. First ncol(iphi)       */
/*           values correspond to AR parms, next ncol(itheta)*/
/*           values correspond to MA parms, and remaining    */
/*           are regression coefficients.                    */
/*  init   - undefined or zero for first call to ARMAREG.    */
/*  maxit  - maximum number of iterations. No other          */
/*           convergence criterion is used. You can invoke   */
/*           ARMAREG without changing parameter values to    */
/*           continue iterations.                            */
/*  nopr   - undefined or zero implies no printing of        */
/*           intermediate results.                           */
/*                                                           */
/*  Notes: Optimization using Gauss-Newton iterations        */
/*                                                           */
/*  Invertibility and stationarity are not checked during    */
/*  theestimation process. The parameter array PAR can be    */
/*  modified after running ARMAREG to place estimates        */
/*  in the stationary and invertible regions, and then       */
/*  ARMAREG can be run again. If a nonstationary AR operator */
/*  is employed, a PAUSE will occur after calling ARMALIK    */
/*  because of a detected singularity. Using STOP will       */
/*  permit termination of ARMAREG so that the AR             */
/*  coefficients can be modified.                            */
/*                                                           */
/*  T-ratios are only approximate and can be undependable,   */
/*  especially for small series.                             */
/*                                                           */
/*  The notation is the same as for the ARMALIK function.    */
/*  The autoregressive and moving average coefficients have  */
/*  signs opposite those given by PROC ARIMA.                */

/* Begin ARMA estimation modules */

/* Generate residuals */
start gres;
   noise=y-x*beta;
   previous=noise[:];
   if ml=-1 then do;                       /* Conditional LS */
      noise=j(nrow(y),1,previous)//noise;
      resid=product(phi,noise`)[,nrow(y)+1:nrow(noise)];
      resid=ratio(theta,resid,ncol(resid));
      resid=resid[,1:ncol(resid)]`;
   end;
   else do;                            /* Maximum likelihood */
      free l;
      call armalik(l,resid,std,noise,phi,theta);

      /* Nonstationary condition produces PAUSE */
      if nrow(l)=0 then do;
         print
         'In GRES: Parameter estimates outside stationary region.';
      end;
      else do;
         temp=l[3,]/(2#nrow(resid));
         if ml=1 then resid=resid#exp(temp);
      end;
   end;
finish gres;                          /* finish module GRES  */

start getpar;                             /* get parameters  */
   if np=0 then phi=1;
   else do;
      temp=parm[,1:np];
      phi=1||j(1,p,0);
      phi[,iphi] =temp;
   end;
   if nq=0 then theta=1;
   else do;
      temp=parm[,np+1:np+nq];
      theta=1||j(1,q,0);
      theta[,itheta] =temp;
   end;
   beta=parm[,(np+nq+1):ncol(parm)]`;
finish getpar;   /* finish module GETPAR  */


/* Get SS Matrix - First Derivatives */
start getss;
   parm=par;
   run getpar;
   run gres;
   s=resid;
   oldsse=ssq(resid);
   do k=1 to ncol(par);
      parm=par;
      parm[,k]=parm[,k]+delta;
      run getpar;
      run gres;
      s=s || ((resid-s[,1])/delta);    /* append derivatives */
   end;
   ss=s`*s;
   if nopr^=0 then print ss[L='Gradient Matrix'];
   sssave=ss;
   do k=1 to 20;           /* Iterate if no reduction in SSE */
      do ii=2 to ncol(ss);
         ss[ii,ii]=(1+lambda)*ss[ii,ii];
      end;
      ss=sweep(ss,2:ncol(ss));       /* Gaussian elimination */
      delpar=ss[1,2:ncol(ss)];     /* update parm increments */
      parm=par+delpar;
      run getpar;
      run gres;
      sse=ssq(resid);
      ss=sssave;
      if sse<oldsse then do;      /* reduction, no iteration */
         lambda=max(lambda/10,1e-12);
         k=21;
      end;
      else do;                               /* no reduction */
                              /* increase lambda and iterate */
         if nopr^=0 then
            print lambda[L='Lambda='] sse oldsse,
                  ss[L='Gradient Matrix'];
         lambda=min(10*lambda,1e12);
         if k=20 then do;
            print ({'GETSS: No improvement in SSE after twenty iterations.',
                    'Possible Ridge Problem.'});
            return;
         end;
      end;
   end;
   if nopr^=0 then print ss[L='Gradient Matrix'];
finish getss;                     /* Finish module GETSS  */

start armareg;                        /* ARMAREG main module */
   /* Initialize options and parameters */
   if nrow(delta)=0 then   delta=0.005;
   if nrow(maxiter)=0 then maxiter=5;
   if nrow(nopr)=0 then    nopr=0;
   if nrow(ml)=0 then      ml=1;
   if nrow(init)=0 then    init=0;
   if init=0 then do;
      p=max(iphi);
      q=max(itheta);
      np=ncol(iphi);
      nq=ncol(itheta);

      /* Make indices one-based */
      do k=1 to np;
         iphi[,k]=iphi[,k]+1;
      end;
      do k=1 to nq;
         itheta[,k]=itheta[,k]+1;
      end;

      /* Create row labels for Parameter estimates */
      if p>0 then parmname = concat("AR",char(1:p,2));
      if q>0 then parmname = parmname||concat("MA",char(1:q,2));
      parmname = parmname||concat("B",char(1:ncol(x),2));

      /* Create column labels for Parameter estimates */
      pname = {"Estimate" "Std. Error" "T-Ratio"};
      init=1;
   end;

   /* Generate starting values */
   if nrow(par)=0 then do;
      beta=inv(x`*x)*x`*y;
      if np+nq>0 then par=j(1,np+nq,0)||beta`;
      else par=beta`;
   end;
   print par [colname=parmname L='Parameter Starting Values'];
   lambda=1e-6;                        /* Controls step size */
   do iter=1 to maxiter;            /* Do maxiter iterations */
      run getss;
      par=par+delpar;
      if nopr^=0 then do;
         print par[colname=parmname L='Parameter Update'];
         print lambda[L='Lambda='];
      end;
   end;

   sighat=sqrt(sse/(nrow(y)-ncol(par)));
   print sighat[L='Innovation Standard Deviation'];
   ss=sweep(ss,2:ncol(ss));          /* Gaussian elimination */
   estm=par`||(sqrt(diag(ss[2:ncol(ss),2:ncol(ss)]))
        *j(ncol(par),1,sighat));
   estm=estm||(estm[,1] /estm[,2]);
   if ml=1 then label='Maximum Likelihood Estimation Results';
   else if ml=-1 then label='Conditional Least Squares Estimation Results';
   else label='Unconditional Least Squares Estimation Results';
   print estm [rowname=parmname colname=pname L=label] ;
finish armareg;
/* End of ARMA Estimation modules */

/* Begin estimation for Grunfeld's investment models */
use grunfeld;
read all var {gei} into y;
read all var {gef gec} into x;
close grunfeld;

x=j(nrow(x),1,1)||x;
iphi=1;
itheta=1;
maxiter=10;
delta=0.0005;
ml=-1;
/*----  To prevent overflow, specify starting values  ----*/
par={-0.5  0.5  -9.956306  0.0265512  0.1516939};
run armareg;   /*----  Perform CLS estimation  ----*/

The results are shown in Output 9.13.1.

/*----  With CLS estimates as starting values,  ----*/
/*----  perform ML estimation.                  ----*/
ml=1;
maxiter=10;
run armareg;

The results are shown in Output 9.13.2.