/* This program looks at the deterministic trend producting correct (minimum mean square error)
   predictions and correct prediction confidence interval.  Then the incorrect predictions are generated
   by a stochastic trend model (ARIMA(0,1,0).  As a followup the ADF and Phillips-Perron unit root tests
   are conducted to see if the trend is stochastic or deterministic.  The tests support the deterministic
   trend model.  */  

data mc;

x1 = 0;
iseed = 7654321;

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

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

if t > 0 then output;

end;

keep x2;
run;

data mc;
  set mc;
  y = x2;
  time = _N_;
  run;

title 'Deterministic Trend Data';

proc gplot data=mc;
  symbol v=dot c=black i=join h=.8;
  title1 'Deterministic Trend Data';
  title2 'X=Time Y=Trend Series';
  axis1 order=(0 to 100 by 10)
  label=(f=duplex 'Obs');
  axis2 order=(0 to 450 by 25)
  label=(f=duplex 'Trend Series');
  plot y*time / haxis=axis1 vaxis=axis2;

  run;
  title;

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

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

proc reg data=post;
  model y = time;
  output out=result p = yhat ucl = ucl lcl = lcl r=det_resid; 
  run;

proc gplot data=result;
title1 'Forecasts Based on Deterministic Trend';
title2 'X=Time Y=RW Series';
plot y*time=1 yhat*time=2 
lcl*time=3 ucl*time=3 /
overlay href=100.5;
symbol1 v=star i=none h=0.5;
symbol2 v=circle i=join h=0.5;
symbol3 v=none i=join h=0.5;
run;

title;

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

proc gplot data=pred;
      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;

/* Conducting Various Unit Root Tests - Augmented Dickey-Fuller test and Phillips 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 three lags or two lags.  In the
       case of three augmenting terms the third augmenting term is significant at
       the 10% level but not at the 5% level.  Therefore, one should look at the lag = 3 and
       lag = 2 choices in the Trend case.  In either case the null of a unit root is rejected
       by both the ADF and Phillips tests and therefore we accept the alternative hypothesis of
       a deterministic trend.  */

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;