/* This program generates some graphs of monte carlo data for use in learning
   the various Unit Root Testing cases.  */

/*  This segment of the program generates a Random Walk without Drift and Zero Starting Value.  */

data mc1;

x1 = 0;
iseed = 8888;

do t = 0 to 100;
  a = rannor(iseed);

x2 = 1.0*x1 + a;

if t > 0 then output;

x1 = x2;

end;

keep t x2;
run;

data mc1;
  set mc1;
  time = t;
  y = x2;
  run;

proc gplot data=mc1;
  symbol v=dot c=black i=join h=.8;
  title1 'Monte Carlo Random Walk Without Drift and Zero Starting Value';
  title2 'Case 1 Null Hypothesis of Unit Root';
  axis1 order=(0 to 100 by 10)
  label=(f=duplex 'Time');
  axis2 order=(-10 to 10 by 2.0)
  label=(f=duplex 'RW Series'); 
  plot y*time / haxis=axis1 vaxis=axis2 overlay vref = 0; 
  run;

/* This segment of the program generates the stationary AR(1) model with zero mean.  */

data mc2;

x1 = 0;
iseed = 12345;

do t = 0 to 100;
  a = rannor(iseed);

x2 = 0.6*x1 + a;

if t > 0 then output;

x1 = x2;

end;

keep t x2;
run;

data mc2;
  set mc2;
  time = t;
  y = x2;
  run;

proc gplot data=mc2;
  symbol v=dot c=black i=join h=.8;
  title1 'Monte Carlo AR(1) data with phi(0) = 0 and phi(1) = 0.6';
  title2 'Case 1 Alternative Hypothesis of Stationarity';
  axis1 order=(0 to 100 by 10)
  label=(f=duplex 'Time');
  axis2 order=(-5 to 5 by 1.0)
  label=(f=duplex 'AR(1) Series');
  plot y*time / haxis=axis1 vaxis=axis2 overlay vref = 0; 
  run;

  title;

/*  This segment of the program generates a Random Walk without Drift but a non-zero starting value.  */

data mc3;

x1 = 10;
iseed = 8888;

do t = 0 to 100;
  a = rannor(iseed);

x2 = 0.0 + 1.0*x1 + 2*a;

if t > 0 then output;

x1 = x2;

end;

keep t x2;
run;

data mc3;
  set mc3;
  time = t;
  y = x2;
  run;

proc gplot data=mc3;
  symbol v=dot c=black i=join h=.8;
  title1 'Monte Carlo Random Walk Without Drift But Non-zero Starting Value';
  title2 'Case 2 Null Hypothesis of Unit Root';
  axis1 order=(0 to 100 by 10)
  label=(f=duplex 'Time');
  axis2 order=(-10 to 40 by 10)
  label=(f=duplex 'RW Series');
  plot y*time / haxis=axis1 vaxis=axis2 overlay vref = 10;
  run;

/*  This segment of the program generates an AR(1) with mean of 10 and AR(1) coefficient of 0.8.  */

data mc4;

x1 = 10;
iseed = 12345;

do t = 0 to 100;
  a = rannor(iseed);

x2 = 4.0 + 0.6*x1 + a;

if t > 0 then output;

x1 = x2;

end;

keep t x2;
run;

data mc4;
  set mc4;
  time = t;
  y = x2;
  run;

proc gplot data=mc4;
  symbol v=dot c=black i=join h=.8;
  title1 'Monte Carlo AR(1) data with phi(0) = 4 and phi(1) = 0.6';
  title2 'Case 2 Alternative Hypothesis of Startionarity with Non-zero mean';
  axis1 order=(0 to 100 by 10)
  label=(f=duplex 'Time');
  axis2 order=(0 to 20 by 5)
  label=(f=duplex 'AR(1) Series');
  plot y*time / haxis=axis1 vaxis=axis2 overlay vref = 10; 
  run;

  title;

/*  This segment of the program generates a Random Walk with Drift.  */

data mc5;

x1 = 10;
iseed = 8888;

do t = 0 to 100;
  a = rannor(iseed);

x2 = 4.0 + 1.0*x1 + 10*a;

if t > 0 then output;

x1 = x2;

end;

keep t x2;
run;

data mc5;
  set mc5;
  time = t;
  y = x2;
  run;

title1 'Monte Carlo Random Walk with Drift Data';
title2 'Case 3 Null Hypothesis of Unit Root';
proc sgplot data=mc5 noautolegend;
   reg x=time y=y/lineattrs=(color=black) markerattrs=(size=5.0);
   run;

  title;

/*  This segment of the program generates Deterministic Trend data.  */

data mc6;

x1 = 0;
iseed = 8888;

do t = 0 to 100;
  a = rannor(iseed);

x2 = 10.0 + 4.0*t + 10*a;

if t > 0 then output;

end;

keep t x2;

run;

data mc6;
  set mc6;
  time = t;
  y = x2;
  run;

title1 'Monte Carlo Deterministic Trend Data';
title2 'Case 3 Alternative Hypothesis of Deterministic Trend Data';
proc sgplot data=mc6 noautolegend;
   reg x=time y=y/lineattrs=(color=black) markerattrs=(size=5.0);
   run;

 /*  Here we use the Augmented Dickey Fuller Tests and Phillips-Perron Tests to test for unit roots
     in the various cases.  */

Title 'Case 1 Test of Random Walk without Drift';
proc arima data = mc1;
    identify var = y stationarity = (adf=4);
    identify var = y stationarity = (phillips=4);
run; 

Title 'Case 1 Test of Stationary AR(1) Series';
proc arima data = mc2;
    identify var = y stationarity = (adf=4);
    identify var = y stationarity = (phillips=4);
run;

Title 'Case 2 Test of Random Walk without Drift';
proc arima data = mc3;
    identify var = y stationarity = (adf=4);
    identify var = y stationarity = (phillips=4);
run; 

Title 'Case 2 Test of Stationary AR(1) Series';
proc arima data = mc4;
    identify var = y stationarity = (adf=4);
    identify var = y stationarity = (phillips=4);
run;

Title 'Case 3 Test for Random Walk with Drift';
proc arima data = mc5;
    identify var = y stationarity = (adf=4);
    identify var = y stationarity = (phillips=4);
run; 

Title 'Case 3 Test for Trend Stationary Series';
proc arima data = mc6;
    identify var = y stationarity = (adf=4);
    identify var = y stationarity = (phillips=4);
run;