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;