/*-------------------------------------------------------------*/
 /*----  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
;

proc iml;
reset noname;
/*-----------------------------------------------------------*/
/*  name: ARMAREG Modules                                    */
/*  purpose: Perform Estimation for regression model with    */
/*  ARMA errors                                              */
/*  usage: Before invoking the command                       */
/*                                                           */
/*  run armareg;                                             */
/*                                                           */
/*  define the 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        */
/*                                                           */
/*  No checking for invertibility or stationarity during     */
/*  estimation 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 follows that of the IML function ARMALIK;   */
/*  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`) [, 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 ,'Gradient Matrix', ss;
   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);
      if sse<oldsse then
      do;                         /* reduction, no iteration */
         lambda=lambda/10;
         k=21;
      end;

      else do;                               /* no reduction */
                              /* increase lambda and iterate */
         if nopr^=0 then print ,
            'Lambda=' lambda 'SSE=' sse 'OLDSSE=' oldsse,
            'Gradient Matrix', ss ;
         lambda=10*lambda;
         ss=sssave;
         if k=20 then
         do;
            print 'In module GETSS:
                   No improvement in SSE after twenty iterations.';
            print ' Possible Ridge Problem. ';
            return;
         end;
      end;
   end;
   if nopr^=0 then print ,'Gradient Matrix', ss;
   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:p,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 ,'Parameter Starting Values',;
   print par [colname=parmname];            /* stderr tratio */
   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 ,'Parameter Update',;
         print par [colname=parmname];      /* stderr tratio */
         print ,'Lambda=' lambda,;
      end;
      if abs(par[,1] )>1 then par[,1] =-.8;
   end;

   sighat=sqrt(sse/(nrow(y)-ncol(par)));
   print ,'Innovation Standard Deviation:' sighat;
   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 print ,'Maximum Likelihood Estimation Results',;
   else if ml=-1 then print ,
      'Conditional Least Squares Estimation  Results',;
   else print ,'Unconditional Least Squares Estimation Results',;
   print estm [rowname=parmname colname=pname] ;
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;
 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  ----*/