| 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 written in IML in the following program. 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 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
;
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, */
/* 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 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 ----*/
The results are shown in Output 8.13.1.
Output 8.13.1: Conditional Least Squares Results
ml=1;
maxiter=10;
/*---- With CLS estimates as starting values, ----*/
/*---- perform ML estimation. ----*/
run armareg;
The results are shown in Output 8.13.2.
Output 8.13.2: Maximum Likelihood Results| AR 1 | MA 1 | B 1 | B 2 | B 3 |
| -0.071148 | 0.7850862 | -7.530983 | 0.0402554 | 0.0992474 |
| Innovation Standard Deviation: | 22.667286 |
| 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 |
Copyright © 2009 by SAS Institute Inc., Cary, NC, USA. All rights reserved.