Chapter Contents Previous Next
 General Statistics Examples

## Example 8.13: Parameter Estimation for a Regression Model with ARMA Errors

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

```   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,       */
/*           omitting index 0 corresponding to the zero      */
/*           order constant one.                             */
/*  itheta - defines indices of nonzero MA parameters,       */
/*           omitting index 0 corresponding 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 may be       */
/*  modified after running armareg to place estimates        */
/*  in the stationary and invertible regions, and then       */
/*  armareg may 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 may be modified.                            */
/*                                                           */
/*  T-ratios are only approximate and may 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,
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 */
/* See SAS/ETS User's Guide, Version 5 Edition       */

/* Access Grunfeld data as produced by DATA step in  */
/* SAS/ETS User's Guide Example                      */
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;
run armareg; /* Perform CLS estimation */
```

The output is shown below.

 AR 1 MA 1 B 1 B 2 B 3 0 0 -9.956306 0.0265512 0.1516939

 Innovation Standard Deviation: 18.6391

 Estimate Std. Error T-Ratio AR 1 -0.071148 0.3248431 -0.219022 MA 1 1.2737455 0.2319205 5.4921656 B 1 -7.530983 20.447977 -0.3683 B 2 0.0402554 0.0170277 2.3641054 B 3 0.0992474 0.0354776 2.7974682

```      /* Transform MA coefficient to invertible region */
par[,2]=1/par[,2];
maxiter=2;
run armareg; /* Continue CLS estimation */
```

The output is shown below.

 AR 1 MA 1 B 1 B 2 B 3 -0.071148 0.7850862 -7.530983 0.0402554 0.0992474

 Innovation Standard Deviation: 22.6673

 Estimate Std. Error T-Ratio AR 1 -0.191675 0.3360688 -0.570345 MA 1 0.7367182 0.2101849 3.5050966 B 1 -19.45436 31.327362 -0.621002 B 2 0.038099 0.0168731 2.2579666 B 3 0.121766 0.0433174 2.8110191

```   ml=1;
maxiter=10;

/* With CLS estimates as starting values, */
/* perform ML estimation.                 */
run armareg;
```

The output is shown below.

 AR 1 MA 1 B 1 B 2 B 3 -0.191675 0.7367182 -19.45436 0.038099 0.121766

 Innovation Standard Deviation: 23.0393

 Estimate Std. Error T-Ratio AR 1 -0.196224 0.3510867 -0.558904 MA 1 0.6816035 0.2712038 2.5132518 B 1 -26.47514 33.752825 -0.784383 B 2 0.0392213 0.0165545 2.3692243 B 3 0.1310306 0.0425996 3.0758634

 Chapter Contents Previous Next Top