/* This is intervention data taken from the textbook Forecasting Principles
   and Applications by Stephen A. DeLurgio, p. 508, Irwin/McGraw-Hill, 1998.
   This is a temporary intervention.  */
/* week gallons intervention
   n=110, Weekly demand for bottled water prior to, during, and
   after flood from periods 75 to 84.  Water in thousands of gallons. */

data water;
 input week kgals interv;
 datalines;
   1  89.5460  0.0
   2  92.1816  0.0
   3  96.3044  0.0
   4  96.6500  0.0
   5  99.9868  0.0
   6  95.9396  0.0
   7  98.3144  0.0
   8  95.9932  0.0
   9 100.1812  0.0
  10 108.2292  0.0
  11  97.4724  0.0
  12 103.6488  0.0
  13 110.7996  0.0
  14  94.5796  0.0
  15  98.4340  0.0
  16  96.6144  0.0
  17 101.8856  0.0
  18 100.8232  0.0
  19 109.0104  0.0
  20 110.2688  0.0
  21 104.6388  0.0
  22 102.3228  0.0
  23 103.4800  0.0
  24 113.0104  0.0
  25 114.9948  0.0
  26 108.7712  0.0
  27 105.9812  0.0
  28 104.5748  0.0
  29 112.2724  0.0
  30 110.1060  0.0
  31 105.1716  0.0
  32 115.1820  0.0
  33 116.2084  0.0
  34 110.5996  0.0
  35 110.4940  0.0
  36 113.0444  0.0
  37 116.1608  0.0
  38 110.0724  0.0
  39 117.8328  0.0
  40 120.7360  0.0
  41 113.6464  0.0
  42 115.0988  0.0
  43 119.2856  0.0
  44 119.5848  0.0
  45 110.8000  0.0
  46 122.0040  0.0
  47 119.6072  0.0
  48 119.7196  0.0
  49 115.8604  0.0
  50 119.8252  0.0
  51 120.5808  0.0
  52 116.4580  0.0
  53 123.4388  0.0
  54 121.4192  0.0
  55 124.2264  0.0
  56 124.4972  0.0
  57 125.1768  0.0
  58 130.8900  0.0
  59 121.7976  0.0
  60 124.4820  0.0
  61 126.5196  0.0
  62 131.9668  0.0
  63 128.0204  0.0
  64 129.6864  0.0
  65 122.3484  0.0
  66 123.2304  0.0
  67 133.5664  0.0
  68 127.4280  0.0
  69 126.9304  0.0
  70 126.9304  0.0
  71 127.4828  0.0
  72 134.2140  0.0
  73 130.8760  0.0
  74 129.7768  0.0
  75 154.3632  1.0
  76 155.0004  1.0
  77 161.4332  1.0
  78 160.3660  1.0
  79 160.7872  1.0
  80 157.6016  1.0
  81 158.4992  1.0
  82 155.0712  1.0
  83 162.9236  1.0
  84 167.0836  1.0
  85 140.0000  1.0
  86 140.9540  1.0
  87 140.2560  1.0
  88 129.1956  1.0
  89 138.2160  1.0
  90 137.8500  1.0
  91 134.1180  1.0
  92 136.8724  1.0
  93 140.2644  1.0
  94 144.9196  1.0
  95 143.8944  1.0
  96 145.2768  1.0
  97 143.5000  1.0
  98 141.9732  1.0
  99 144.6516  1.0
 100 137.7968  1.0
 101 147.4736  1.0
 102 148.8548  1.0
 103 146.9712  1.0
 104 142.4760  1.0
 105 149.5192  1.0
 106 152.4736  1.0
 107 151.1600  1.0
 108 149.0000  1.0
 109 148.4820  1.0
 110 149.6440  1.0
 ;

 data water;
   set water;
 pulse = (week = 75);

 run;

 /*  Let's look at the data's systematic dynamics before
     the intervention (Obs. 1 - 74). */

 proc arima data = water (obs=74);
   identify var = kgals(1);
   estimate p = 1;
   estimate p = 2;
   estimate q = 1;
   estimate q = 2;
   estimate p = 1 q = 1;

 /*  An MA(1) model for kgals(1) seems to fit the data the best.  */  

   run;


 proc arima data = water;
  identify var = pulse noprint;
  identify var = kgals(1) crosscorr = (pulse) noprint;
  estimate q = 1 input = ((1)/(1,2)pulse) method=ml;
  estimate q = 1 input = (/(1,2)pulse) method=ml;
  estimate q = 1 input = (/(1)pulse) method = ml;
  estimate q = 1 input = (/(2)pulse) method = ml;
  estimate q = 1 input = (pulse) method = ml;

/* This last model seems pretty good.  Its residuals are white
   noise.  For comparison let's try an AR(2)error instead of
   the MA(1) error. */ 

  estimate p = 2 input = (pulse) method=ml;
  
/* This model comes in a close second in my opinion.  The AIC
   criterion is lower but the former model has the lower SBC
   measure.  From the perspective of parsimony, we can go for 
   the MA(1) error specification.  Besides, since the pre-
   intervention dynamics seem to be MA(1) why not stick with
   it.  */
 
  run;

/* Now let's forecast 12 steps ahead assuming that there will
   be no intervention during the forecast horizon.  */

 data add1;
  input week kgals interv;
  if week >= 111 then pulse = 0;
 datalines;
  111 . 1.0
  112 . 1.0
  113 . 1.0
  114 . 1.0
  115 . 1.0
  116 . 1.0
  117 . 1.0
  118 . 1.0
  119 . 1.0
  120 . 1.0
  121 . 1.0
  121 . 1.0
;

 data extrap1;
  set water;

  proc append base = extrap1 data = add1;

  run;

/*  Here we generate forecasts assuming that there will be no
    no intervention during the forecast horizon.  */

proc arima data = extrap1;
  identify var = pulse noprint;
  identify var = kgals(1) crosscorr = (pulse) noprint;
  estimate q = 1 input = (pulse) method = ml noprint;
  forecast lead = 12;

  run;

/*  Now let's assume that we want to forecast with an intervention.  */

data add2;
  input week kgals interv;
  pulse = (week = 111);
 datalines;
  111 . 1.0
  112 . 1.0
  113 . 1.0
  114 . 1.0
  115 . 1.0
  116 . 1.0
  117 . 1.0
  118 . 1.0
  119 . 1.0
  120 . 1.0
  121 . 1.0
  121 . 1.0
;

data extrap2;
  set water;


proc append base = extrap2 data = add2;

  run;

/*  Here we generate forecasts assuming that there will be an
    at the beginning (obs = 111) of the forecast horizon.  */

proc arima data = extrap2;
  identify var = pulse noprint;
  identify var = kgals(1) crosscorr = (pulse) noprint;
  estimate q = 1 input = (pulse) method = ml;
  forecast lead = 12;

  run;