A Regression Model with ARMA Errors
/****************************************************************/
/* S A S S A M P L E L I B R A R Y */
/* */
/* NAME: ARMAEST */
/* TITLE: A Regression Model with ARMA Errors */
/* PRODUCT: IML */
/* SYSTEM: ALL */
/* KEYS: MATRIX TIME SUGI6 */
/* PROCS: IML */
/* DATA: */
/* */
/* SUPPORT: TJW UPDATE: 22APR90 */
/* REF: Sep 2013 */
/* */
/****************************************************************/
/*-------------------------------------------------------------*/
/*---- 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 ----*/
/*---- With CLS estimates as starting values, ----*/
/*---- perform ML estimation. ----*/
ml=1;
maxiter=10;
run armareg;
