/*-----------------------------------------------------------------
     Example: Estimating GARCH Models
     Requires: SAS/ETS
     Version: 9.0                                                                       
     ------------------------------------------------------------------*/

   %let df = 7.5;
   %let sig1 = 1;
   %let sig2 = 0.1 ;
   %let var2 = 2.5;
   %let nobs = 1000 ;
   %let nobs2 = 2000 ;
   %let arch0 = 0.1 ;
   %let arch1 = 0.2 ;
   %let garch1 = 0.75 ;
   %let intercept = 0.5 ;
   
   data normal;
      lu = &var2;
      lh = &var2;

      do i= -500  to &nobs ;
              /* GARCH(1,1) with normally distributed residuals */
         h = &arch0 + &arch1*lu**2 + &garch1*lh;
         u = sqrt(h) * rannor(12345) ;
         y = &intercept + u;
         lu = u;
         lh = h;
         if i > 0 then output;
      end;
   run;




   data t;
      lu = &var2 ;
      lh = &var2 ;
      do i = -500 to &nobs2;
              /* GARCH(1,1) with t-distributed residuals */
         t = &sig1*tinv(ranuni(1234),&df) ;
         h = &arch0 + &arch1*lu**2 + &garch1*lh;
         u = sqrt( h) * t ;
         y = &intercept + u ;
         lu = u;
         lh = h;
         if i > 0 then output;
      end;
   run;

      

   data cauchy;
      lu = &var2 ;
      lh = &var2 ;
      do i = -500  to &nobs ;
          /* GARCH(1,1) with Cauchy distributed residuals */   
         cauchy = &sig2*rancau(1234);
         h = &arch0 + &arch1*lu**2 + &garch1*lh;
         u = sqrt( h) * cauchy ;
         y = &intercept + u ;
         lu = u;
         lh = h;
         if i > 0 then output;
      end;
   run;



   data garchm;
      lu = &var2;
      lh = &var2;
      gamma = 0.5;
      do i= -500  to &nobs ;
              /* GARCH-M  */
         h = &arch0 + &arch1*lu**2 + &garch1*lh  ;
         u = sqrt(h) * rannor(1234) ;
         y = &intercept  + gamma*sqrt(h) + u;
         lu = u;
         lh = h;
         if i > 0 then output;
      end;
   run;



   data egarch ;
      lu = &var2 ;
      lh = &var2 ;
      theta = .65 ;
      lz = lu/sqrt(lh) ;
      lg = theta*lz + abs(lz)-sqrt(2/3.14159) ;
      do i = -500 to &nobs ;
             /* EGARCH */
         h = exp( &arch0 + &arch1*lg + &garch1*log(lh)) ;
         u = sqrt(h)*rannor(12346) ;
         z = u/sqrt(h) ;
         g = theta*z + abs(z) -sqrt(2/3.14159) ;
         y = &intercept + u ;
         lu = u ;
         lh = h ;
         lg = g ;
         if i > 0 then output ;
      end ;
   run ;    



   data qgarch;
      lu = &var2;
      lh = &var2;
      phi = .2;
      do i= -500  to &nobs ;
              /* Quadratic GARCH  */
         h = &arch0 + &arch1*lu**2 + &garch1*lh + phi*lu ;
         u = sqrt(h) * rannor(1234) ;
         y = &intercept  + u;
         lu = u;
         lh = h;
         if i > 0 then output;
      end;
   run;  

    

   data gjrgarch;
      lu = &var2;
      lh = &var2;
      phi = 0.1;
      do i= -500  to &nobs ;
              /* GJR-GARCH  */
         if lu >= 0 then
            h = &arch0 + &arch1*lu**2 + &garch1*lh + phi*lu**2 ;
         else
            h = &arch0 + &arch1*lu**2 + &garch1*lh ;
         u = sqrt(h) * rannor(1234) ;
         y = &intercept  + u;
         lu = u;
         lh = h;
         if i > 0 then output;
      end;
   run;


   data tgarch;
      lu = &var2;
      lh = &var2;
      arch1_plus = 0.1;
      arch1_minus = 0.1;
      do i= -500 to &nobs ;
              /* TGARCH  */
         if lu > 0 then
            h = (&arch0 + arch1_plus*lu + &garch1*sqrt(lh))**2 ;
         else
            h = (&arch0 + arch1_minus*lu + &garch1*sqrt(lh))**2 ;
         u = sqrt(h) * rannor(1234) ;
         y = &intercept  + u;
         lu = u;
         lh = h;
         if i > 0 then output;
      end;
   run;



   proc autoreg data = normal ;
    /* Estimate GARCH(1,1) with normally distributed residuals with AUTOREG*/  
      model y = / garch = ( q=1,p=1 ) ;
   run ;
   quit ;

    

    /* Estimate GARCH(1,1) with normally distributed residuals with MODEL*/    
   proc model data = normal ;
       parms arch0 .1 arch1 .2 garch1 .75 ;
       /* mean model */
       y = intercept ;
       /* variance model */
       h.y = arch0 + arch1*xlag(resid.y**2,mse.y) +
             garch1*xlag(h.y,mse.y) ;
       /* fit the model */
       fit y / method = marquardt fiml ;
   run ;
   quit ;



    /* Estimate GARCH(1,1) with t-distributed residuals with AUTOREG*/
   proc autoreg data = t ;
      model y = / garch=( q=1, p=1 ) dist = t ;
   run;
   quit;



    /* Estimate GARCH(1,1) with t-distributed residuals with MODEL*/    
   proc model data = t ;
      parms   df 7.5 arch0 .1 arch1 .2 garch1 .75 ;
      /* mean model */
      y = intercept ;
      /* variance model */
      h.y = arch0 + arch1 * xlag(resid.y **2, mse.y)  +
            garch1*xlag(h.y, mse.y);
      /* specify error distribution */
      errormodel y ~ t(h.y,df);
      /* fit the model */
      fit y / method=marquardt;
   run;
   quit;



   /* Estimate GARCH(1,1) with Cauchy distributed residuals */ 
   proc model data = cauchy ;  
      parms arch0 .1 arch1 .2 garch1 .75 intercept .5 ;
      control mse.y = &var2 ; /*defined in data generating step*/
      /* mean model */
      y = intercept ;
      /* variance model */
      h.y = arch0 + arch1 * xlag(resid.y ** 2, mse.y)  +
            garch1 * xlag(h.y, mse.y);
      /* specify error distribution */
      obj = log(h.y/((h.y**2+resid.y**2) * constant('pi')));
      obj = -obj ;
      errormodel y ~ general(obj);
      /* fit the model */
      fit y / method=marquardt;
   run;
   quit;


   /* Estimate GARCH(1,1) with generalized error distribution residuals */ 
   proc model data = normal ;
      parms nu 2 arch0 .1 arch1 .2 garch1 .75;
      control mse.y = &var2 ; /*defined in data generating step*/
      /* mean model */
      y = intercept ;
      /* variance model */
      h.y = arch0 + arch1 * xlag(resid.y ** 2, mse.y)  +
            garch1 * xlag(h.y, mse.y);
      /* specify error distribution */
      lambda = sqrt(2**(-2/nu)*gamma(1/nu)/gamma(3/nu)) ;
      obj = log(nu/lambda) -(1 + 1/nu)*log(2) - lgamma(1/nu)-
            .5*abs(resid.y/lambda/sqrt(h.y))**nu - .5*log(h.y) ;
      obj = -obj ;
      errormodel y ~ general(obj,nu);
      /* fit the model */
      fit y / method=marquardt;
   run;
   quit;

    

   /* Estimate GARCH-M Model with PROC AUTOREG */ 
   proc autoreg data= garchm ;
      model y =  /  garch=( p=1, q=1,  mean = sqrt);
   run; 
   quit; 



   /* Estimate GARCH-M Model with PROC MODEL */     
   proc model data = garchm ;
      parms  arch0 .1 arch1 .2 garch1 .75 gamma .5 ;
      h = arch0 + arch1*xlag(resid.y**2,mse.y) + garch1*xlag(h.y,mse.y);
      y = intercept + gamma*sqrt(h) ;
      h.y = h ;
      fit y / fiml method = marquardt;
   run;
   quit;


   /* Estimate EGARCH Model with PROC AUTOREG */ 
   proc autoreg data= egarch ;
      model y =  / garch=( q=1, p=1 , type = exp) ;
   run;
   quit;



   /* Estimate EGARCH Model with PROC MODEL */     
   proc model data = egarch ;
      parms  earch0 .1 earch1 .2 egarch1 .75 theta .65 ;
      /* mean model */
      y = intercept ;
      /* variance model */
      if (_obs_ = 1 ) then
         h.y = exp(earch0/(1. - egarch1)) ;
      else h.y = exp(earch0 + earch1*zlag(g) + egarch1*log(zlag(h.y))) ;
      g = theta*(-nresid.y) + abs(-nresid.y) - sqrt(2/constant('pi')) ;
      /* fit the model */
      fit y / fiml method = marquardt ;
   run;
   quit;



   /* Estimate Quadratic GARCH (QGARCH) Model */ 
   proc model data = qgarch ;
      parms arch0 .1 arch1 .2 garch1 .75 phi .2;
      /* mean model */
      y = intercept ;
      /* variance model */
      h.y = arch0 + arch1*xlag(resid.y**2,mse.y) + garch1*xlag(h.y,mse.y) +
            phi*xlag(-resid.y,mse.y);
      /* fit the model */
      fit y / method = marquardt fiml ;
   run ;
   quit ;


    

   /* Estimate GJR-GARCH Model */         
   proc model data = gjrgarch ;
      parms arch0 .1 arch1 .2 garch1 .75 phi .1;
      /* mean model */
      y = intercept ;
      /* variance model */
      if zlag(resid.y) > 0 then
         h.y = arch0 + arch1*xlag(resid.y**2,mse.y) + garch1*xlag(h.y,mse.y)  ;
      else
         h.y = arch0 + arch1*xlag(resid.y**2,mse.y) + garch1*xlag(h.y,mse.y) +
               phi*xlag(resid.y**2,mse.y) ;
      /* fit the model */
      fit y / method = marquardt fiml ;
   run ;
   quit ;




   /* Estimate Threshold Garch (TGARCH) Model */ 
   proc model data = tgarch ;
      parms arch0 .1 arch1_plus .1 arch1_minus .1 garch1 .75 ;
      /* mean model */
      y = intercept ;
      /* variance model */
      if zlag(resid.y) < 0 then
         h.y = (arch0 + arch1_plus*zlag(-resid.y) + garch1*zlag(sqrt(h.y)))**2 ;
      else
         h.y = (arch0 + arch1_minus*zlag(-resid.y) + garch1*zlag(sqrt(h.y)))**2 ;
      /* fit the model */
      fit y / method = marquardt fiml ;
   run ;
   quit ;