• Print  |
  • Feedback  |

FOCUS AREAS

SAS/ETS Examples

Bootstrapping Correct Critical Values in Tests for Structural Change

Contents | Back to Example

SAS Program 

   /*-----------------------------------------------------------------------------
     Example: Bootstrapping Correct Critical Values in Tests for Structural Change
     Requires: SAS/ETS, SAS/GRAPH, SAS/IML
     Version: 9.0                                                                       
     -----------------------------------------------------------------------------*/

   data gdp;
      set sashelp.citiqtr;
      keep date gdpq l_gdpq;
      l_gdpq = log(gdpq);
   run;

   title 'Log of Quarterly GDP';
   title2 '1980 - 1991';

   axis3 label=('Year')
         order=('01jan80'd to '01jan92'd by year);
   axis4 label=(angle=90 'Log GDP')
         order=(8.2 to 8.5 by .1);

   proc gplot data=gdp;
      symbol1 c=red i=join;
      format date year4.;
      plot l_gdpq * date/ vaxis=axis4 vminor=3
                          haxis=axis3 cframe=ligr;
   run;
   quit;


   proc iml;

   simnum = 10000;
   echo   = 500;

   use gdp;
   read all into data var{gdpq};
   read all into date var{date};
   ldata = log(data);

   nlag = 2;
   strtvls = ldata[1:nlag,];

   y    = ldata[3:nrow(ldata)];
   y_1  = ldata[2:nrow(ldata)-1];
   y_2  = ldata[1:nrow(ldata)-2];



   start battery(y,y_1,y_2) global(xr, k, m);
      batf = 0;
      n    = nrow(y);
      t    = cusum(j(n,1,1));
      xr   = j(n,1,1) || t || y_1 || y_2;
      m    = ncol(xr);
      er   = (i(n) - xr*inv(xr`*xr)*xr`)*y;
      rss  = er`*er;
      do i=(m-1) to n-2 by 1;
         di    = j(i-1,1,0) // j(n-i+1,1,1);
         xu    = j(n,1,1) || di || t || di#t || y_1 || y_2;
         k     = ncol(xu);
         eu    = (i(n) - xu*inv(xu`*xu)*xu`)*y;
         uss   = eu`*eu;
         fstat = ((rss-uss)/(k-m)) / (uss/(n-k));
         batf  = batf // fstat;
      end;
      batf = batf[2:nrow(batf),];
      return(batf);
   finish battery;



   start bootstrp(ndata, strtvls, ehat, betahat);
      ner     = nrow(ehat);
      rndehat = j(ner,1,1);
      do eloop=1 to ner by 1;
         rndehat[eloop] = ehat[ceil(ner*uniform(0))];
      end;
      ysim    = j(ndata,1,1);
      ysim[1:nrow(strtvls),] = strtvls;
      t       = cusum(j(ndata,1,1));
      do l=(nrow(strtvls)+1) to ndata by 1;
         ysim[l] = betahat[1] + 
                   betahat[2]*t[l-2] + 
                   betahat[3]*ysim[l-1] +
                   betahat[4]*ysim[l-2] + 
                   rndehat[l-2];
      end;
      return(ysim);
   finish bootstrp;



   start sortcm(x); 
      temp1=x;
      do i=1 to nrow(temp1) by 1;
         temp2 = temp1[i,];
         temp3 = temp2;
         temp2[,rank(temp2)] = temp3;
         temp1[i,] = temp2;
      end;
      return(temp1);
   finish sortcm;



   fdata   = battery(y, y_1, y_2);

   temp    = cusum(j(nrow(fdata)+m-2,1,1));
   brkpt   = temp[(m-1):nrow(temp),];
   ehat    = (i(nrow(y)) - xr*inv(xr`*xr)*xr`)*y;
   betahat = inv(xr`*xr)*xr`*y;



   f     = j(nrow(ehat)-2*(k-m),1,0);
   ndata = nrow(ldata);

   do simloop=1 to simnum by 1;
      ysim = bootstrp(ndata, strtvls, ehat, betahat);
      y   = ysim[3:ndata];
      y_1 = ysim[2:ndata-1];
      y_2 = ysim[1:ndata-2];
      f = f || battery(y, y_1, y_2);
      if mod(simloop,echo)=0 then print simloop;
   end;



   f       = f[,2:ncol(f)];
   fcolmax = sortcm(f[<>,]);
   fsort   = sortcm(f);
   cv95    = int(.95*simnum);

   brkpt   = date[brkpt,];
   fdata   = fdata;
   fstd_95 = j(nrow(fdata),1,finv(.95,k-m,nrow(ehat)-k));
   f_95    = fsort[,cv95];
   fmax_95 = j(nrow(fdata),1,fcolmax[,cv95]);



   create sasuser.bootout 
          var{brkpt fdata fstd_95 f_95 fmax_95};
   append;

   quit;

   /* -------  Graphics Output  ------- */

   title 'Bootstrapped Critical Values';

   axis1 label=('Break Point') minor=none
         order=('01jan80'd to '01jan91'd by year);
   axis2 label=(angle=90 'F-statistics')
         order=(0 to 15 by 3);

   symbol1 c=red   i=join;    /* for fdata   */
   symbol2 c=black i=spline;  /* for fstd_95 */
   symbol3 c=green i=spline;  /* for f_95    */
   symbol4 c=blue  i=spline;  /* for fmax_95 */
   
   proc gplot data=sasuser.bootout;
      format brkpt year4.;
      plot fdata   * brkpt = 1
           fstd_95 * brkpt = 2
           f_95    * brkpt = 3
           fmax_95 * brkpt = 4 / overlay
                                 vaxis=axis2
                                 vminor=1
                                 haxis=axis1
                                 legend
                                 cframe=ligr;
   run;
   quit;