/* This program looks at the RW with drift data first producting correct (minimum mean square error)
   predictions and correct prediction confidence interval.  Then the incorrect predictions are generated
   a deterministic trend model with the size of the prediction appraching a limit as the forecast horizon
   goes to infinity.  As a followup, ADF and Phillips-Perron unit root tests are conducted to see
   if the trend is stochastic or deterministic.  The tests support the RW with drift model.  */  

data mc;

y1 = 0;
iseed = 2468;

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

y2 = 1.0 + 1.0*y1 + 2*a;

if t > 0 then output;

y1 = y2;

end;

keep t y2;

run;

title 'Monte Carlo Random Walk With Drift Data';

proc gplot data=mc;
  symbol v=dot c=black i=join h=.8;
  title1 'Monte Carlo Random Walk With Drift Data';
  title2 'X=Time Y=RW Series';
  axis1 order=(0 to 100 by 10)
  label=(f=duplex 'Obs');
  axis2 order=(-10 to 120 by 10)
  label=(f=duplex 'RW Series');
  plot y2*t / haxis=axis1 vaxis=axis2;
  run;

data mc;
  set mc;
  y=y2;
  time = _N_;
  drop y2;
  run;

  title;

proc arima data=mc;
  identify var = y(1);
  estimate p = 0;
  forecast lead=20 id=time out=pred1;
  run;

proc gplot data=pred1;
      title1 'Forecasts Based on ARIMA(0,1,0) Box-Jenkins Model';
      title2 'X=Time Y=RW Series';
      symbol1 i=none   v=star h=0.5;
      symbol2 i=spline v=circle h=0.5;
      symbol3 i=spline l=5;
      plot y * time = 1 
           forecast * time = 2 
           ( l95 u95 ) * time = 3 /
           overlay href = 100.5;
     run;

	 title;

  data b;
   y = .;
   do time = 101 to 120; output; end;
run;

data post;
   merge mc b;
   by time;
run;

proc autoreg data = post;
   model y = time / nlag = 12 method=ml backstep slstay = 0.05;
   output out = result p = yhat pm = ytrend lcl = lcl ucl = ucl;
run;

proc gplot data=result;
title1 'Forecasts Based on Proc Autoreg Linear Trend with AR(1,8) Errors';
title2 'X=Time Y=RW Series';
plot y*time=1 yhat*time=2 ytrend*time=3
lcl*time=3 ucl*time=3 /
overlay href=100.5;
symbol1 v=star i=none;
symbol2 v=circle i=join;
symbol3 v=none i=join;
run;

title;

/* Conducting Various Unit Root Tests - Augmented Dickey-Fuller test and Phillips-Perron test */

proc arima data = mc;
    identify var = y stationarity = (adf=4);
    identify var = y stationarity = (phillips=4);
run;

data mc;
  set mc;
   y_1 = lag(y);
   ydif = y - lag1(y);
   ydif1 = lag1(ydif);
   ydif2 = lag2(ydif);
   ydif3 = lag3(ydif);
   ydif4 = lag4(ydif);
   run;

   /*  Backward Examination of the lagged differences (augmenting terms) for significance
       The actual number of lags should be zero given the Monte Carlo Experiment.  In the
       backward elimination process one could either choose zero lags or two lags.  In the
       case of two augmenting terms the second augmenting term is significant at
       the 10% level but not at the 5% level.  Therefore, one should look at the lag = 2 and
       lag = 0 choices in the Trend case.  In either case the null of a unit root is supported
       by both the ADF and Phillips-Perron tests.  */

proc reg data=mc;
   model ydif =  y_1 time ydif1 ydif2 ydif3 ydif4;
   model ydif =  y_1 time ydif1 ydif2 ydif3;
   model ydif =  y_1 time ydif1 ydif2;
   model ydif =  y_1 time ydif1;
   model ydif =  y_1 time;
   run;