/* This program generates the Level Shift Data  */

data mc;

iseed = 13367;

do t = 1 to 50;
  a = rannor(iseed);

x2 =  a;

output;

end;

iseed = 745;

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

x2 = 5 + a;

output;

end;

keep t x2;

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

proc gplot data=mc;
  symbol v=dot c=black i=join h=.8;
  title1 'Monte Carlo Level Shift Data';
  title2 'X=Time Y=Level Shift Data';
  axis1 order=(0 to 100 by 10)
  label=(f=duplex 'Obs');
  axis2 order=(-10 to 10 by 2)
  label=(f=duplex 'y');
  plot y*time / haxis=axis1 vaxis=axis2;
  run;

title;

proc arima data=mc;
  identify var = y;
  estimate p = 2;
  forecast lead=20 id=time out=pred1;
  run;

data pred1;
    set pred1;
    fore1 = forecast;
	drop forecast;
	run;

proc gplot data=pred1;
      title1 'Forecasts Based on AR(2) Model based on ALL of the data';
      title2 'X=Time Y=Forecasting While Ignoring the Level Shifta';
      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 
           fore1 * time = 2
           ( l95 u95 ) * time = 3 /
           overlay href = 100.5;
      run;

title;
  
proc arima data=mc(where = (time > 50));
  identify var = y;
  estimate p = 0;
  forecast lead=20 id=time out=pred2;
  run;

  data pred2;
     set pred2;
     fore2 = forecast;
     drop forecast;
     run;

data combine;
    merge pred1 pred2;
       by time;
       run;

proc print data=combine;
   run; 
  
proc gplot data=combine;
      title1 'Forecasts Based on AR(0) on second half of data versus AR(2) on all of the data';
      title2 'X=Time Y=Second Half of Level Shift Data';
      symbol1 i=none   v=star h=0.5;
      symbol2 i=spline v=circle h=0.5;
	  symbol3 i = spline v=diamond h=0.5;
      symbol4 i=spline l=5;
      plot y * time = 1 
           fore1 * time = 2
           fore2 * time = 3  
           ( l95 u95 ) * time = 4 /
           overlay href = 100.5;
           where time > 50;
     run;