EMM Estimation of a Stochastic Volatility Model
/*--------------------------------------------------------------
SAS Sample Library
Name: modex19.sas
Description: Example program from SAS/ETS User's Guide,
The MODEL Procedure
Title: EMM Estimation of a Stochastic Volatility Model
Product: SAS/ETS Software
Keys: nonlinear simultaneous equation models
PROC: MODEL
Notes:
--------------------------------------------------------------*/
%let nobs=1000;
data svdata;
a = -0.736; b=0.9; s=0.363;
ll=sqrt( exp(a/(1-b)));;
do i=-10 to &nobs;
u = rannor( 101 );
z = rannor( 101 );
lnssq = a+b*log(ll**2) +s*u;
st = sqrt(exp(lnssq));
ll = st;
y = st * z;
if i > 0 then output;
end;
run;
title1 'Efficient Method of Moments for Stochastic Volatility Model';
/* estimate GARCH(1,1) model */
proc autoreg data=svdata(keep=y)
outest=garchest
noprint covout;
model y = / noint garch=(q=1,p=1,type=nonneg);
run;
data aregparms;
keep mse w a b;
set garchest(where=(_type_="PARM"));
mse = _mse_;
w = _ah_0;
a = _ah_1;
b = _gh_1;
run;
data parmsout;
if _n_=1 then set aregparms;
set svdata(keep=y);
run;
data scores;
keep y var dlldv dlldw dllda dlldb;
retain lagy2 0
lagvar 0
currdvdb 0
lagdvdw 0
lagdvda 0
lagdvdb 0
laglagdvdb 0;
set parmsout;
if _n_=1 then lagy2=mse;
/* derivative of loglik wrt sigma-sq */
lagvar = w + a*lagy2 + lagvar*b;
var = lagvar + mse*b**_n_;
dlldv = (-1 + y**2/var)/var/2;
/* arch 0 */
dvdw = b*lagdvdw + 1;
dlldw = dlldv*dvdw;
/* arch 1 */
dvda = b*lagdvda + lagy2;
dllda = dlldv*dvda;
/* garch 1 */
dvdb = - b*b*laglagdvdb + 2*b*lagdvdb + currdvdb;
dlldb = dlldv*(dvdb + _n_*b**(_n_-1)*mse);
/* reset lagged values */
currdvdb = w + a*lagy2;
lagy2 = y*y;
lagdvdw = dvdw;
lagdvda = dvda;
laglagdvdb = lagdvdb;
lagdvdb = dvdb;
run;
/* compute the V matrix */
data vvalues;
set scores;
array score{*} dlldw dllda dlldb;
array v_t{*} v_t_1-v_t_6;
array v{*} v_1-v_6;
/* compute external product of score vector */
do i=1 to 3;
do j=i to 3;
v_t{j*(j-1)/2 + i} = score{i}*score{j};
end;
end;
/* average them over t */
do s=1 to 6;
v{s}+ v_t{s}/&nobs;
end;
run;
/* Create a VDATA data set acceptable to PROC MODEL */
/* Transpose the last obs in the data set */
proc transpose data=vvalues(firstobs=&nobs keep=v_1-v_6)
out=tempv;
run;
/* Add eq and inst labels */
data vhat;
set tempv(drop=_name_);
value = col1;
drop col1;
input _type_ $ eq_row $ eq_col $ inst_row $ inst_col $; *$;
datalines;
gmm m1 m1 1 1 /* intcpt is the only inst we use */
gmm m1 m2 1 1
gmm m2 m2 1 1
gmm m1 m3 1 1
gmm m2 m3 1 1
gmm m3 m3 1 1
;
/* USE SMM TO FIND EMM ESTIMATES */
/* Generate data set of length T */
data emm;
set garchest(where=(_type_="PARM") rename=(_ah_0=w _ah_1=a _gh_1=b _mse_=mse)
keep=_type_ _ah_0 _ah_1 _gh_1 _mse_);
do i=1 to 20000;
output;
end;
drop i;
run;
title2 'EMM estimates';
/* Find the EMM estimates */
proc model data=emm maxiter=1000 plot=none;
parms aa -0.736 bb 0.9 ss 0.363;
instruments _exog_ / intonly;
/* Describe the structural model */
u = rannor( 8801 );
z = rannor( 9701 );
lsigmasq = xlag(sigmasq,exp(aa));
lnsigmasq = aa + bb * log(lsigmasq) + ss * u;
sigmasq = exp( lnsigmasq );
ysim = sqrt(sigmasq) * z;
/* Compute scores of GARCH(1,1) */
/* derivative of loglik wrt sigma-sq */
ysim2 = ysim*ysim;
lagvar = w + a*xlag(ysim2,mse) + xlag(lagvar,0)*b;
var = lagvar + mse*b**_n_;
dlldv = (-1 + ysim2/var)/var/2;
/* arch 0 */
dvdw = b*xlag(dvdw,0) + 1;
dlldw = dlldv*dvdw;
/* arch 1 */
dvda = b*xlag(dvda,0) + xlag(ysim2,mse);
dllda = dlldv*dvda;
/* garch 1 */
currdvdb = w + a*xlag(ysim2,mse);
dvdb = - b*b*xlag2(dvdb,0) + 2*b*xlag(dvdb,0) + xlag(currdvdb,0);
dlldb = dlldv*(dvdb + _n_*b**(_n_-1)*mse);
/* Use scores of the GARCH model as moment conditions */
eq.m1 = dlldw;
eq.m2 = dllda;
eq.m3 = dlldb;
/* Fit scores using SMM and estimated Vhat */
fit m1 m2 m3 / gmm npreobs=10 ndraw=1 /* smm options */
vdata=vhat /* use estimated Vhat */
kernel=(bart,0,) /* turn smoothing off */;
bounds ss > 0, 0 < bb < 1;
quit;