/*  This is intervention data taken from the textbook Forecasting
    Pinciples and Applications by Stephen A. Delurgio, Irwin/Mcgraw-
    Hill, p.504, 1998 (Blood.dat).  This data exhibits a one-time level
    shift.   */

/* T=187, Demand at a community blood bank, intervention in period 135.  */


data blood;
input t demand step;
datalines;
   1  48.0512 0.0
   2  54.6040 0.0
   3  54.4152 0.0
   4  48.2036 0.0
   5  51.8904 0.0
   6  51.1796 0.0
   7  52.7856 0.0
   8  47.3224 0.0
   9  47.5040 0.0
  10  49.3248 0.0
  11  44.0460 0.0
  12  46.1816 0.0
  13  49.8044 0.0
  14  49.6500 0.0
  15  52.4868 0.0
  16  47.9396 0.0
  17  49.8144 0.0
  18  46.9932 0.0
  19  50.6812 0.0
  20  58.2292 0.0
  21  46.9724 0.0
  22  52.6488 0.0
  23  59.2996 0.0
  24  42.5796 0.0
  25  45.9340 0.0
  26  43.6144 0.0
  27  48.3856 0.0
  28  46.8232 0.0
  29  54.5104 0.0
  30  55.2688 0.0
  31  49.1388 0.0
  32  46.3228 0.0
  33  46.9800 0.0
  34  56.0104 0.0
  35  57.4948 0.0
  36  50.7712 0.0
  37  47.4812 0.0
  38  45.5748 0.0
  39  52.7724 0.0
  40  50.1060 0.0
  41  44.6716 0.0
  42  54.1820 0.0
  43  54.7084 0.0
  44  48.5996 0.0
  45  47.9940 0.0
  46  50.0444 0.0
  47  52.6608 0.0
  48  46.0724 0.0
  49  53.3328 0.0
  50  55.7360 0.0
  51  48.1464 0.0
  52  49.0988 0.0
  53  52.7856 0.0
  54  52.5848 0.0
  55  43.3000 0.0
  56  54.0040 0.0
  57  51.1072 0.0
  58  50.7196 0.0
  59  46.3604 0.0
  60  49.8252 0.0
  61  50.0808 0.0
  62  45.4580 0.0
  63  51.9388 0.0
  64  49.4192 0.0
  65  51.7264 0.0
  66  51.4972 0.0
  67  51.6768 0.0
  68  56.8900 0.0
  69  47.2976 0.0
  70  49.4820 0.0
  71  51.0196 0.0
  72  55.9668 0.0
  73  51.5204 0.0
  74  52.6864 0.0
  75  44.8484 0.0
  76  45.2304 0.0
  77  55.0664 0.0
  78  48.4280 0.0
  79  47.4304 0.0
  80  46.9304 0.0
  81  46.9828 0.0
  82  53.2140 0.0
  83  49.3760 0.0
  84  47.7768 0.0
  85  46.8632 0.0
  86  47.0004 0.0
  87  52.9332 0.0
  88  51.3660 0.0
  89  51.2872 0.0
  90  47.6016 0.0
  91  47.9992 0.0
  92  44.0712 0.0
  93  51.4236 0.0
  94  55.0836 0.0
  95  52.5000 0.0
  96  52.9540 0.0
  97  51.7560 0.0
  98  40.1956 0.0
  99  48.7160 0.0
 100  47.8500 0.0
 101  43.6180 0.0
 102  45.8724 0.0
 103  48.7644 0.0
 104  52.9196 0.0
 105  51.3944 0.0
 106  52.2768 0.0
 107  60.8528 0.0
 108  47.9732 0.0
 109  50.1516 0.0
 110  42.7968 0.0
 111  51.9736 0.0
 112  52.8548 0.0
 113  50.4712 0.0
 114  45.4760 0.0
 115  52.0192 0.0
 116  54.4736 0.0
 117  52.6600 0.0
 118  37.5564 0.0
 119  48.9820 0.0
 120  49.6440 0.0
 121  49.2392 0.0
 122  44.8048 0.0
 123  56.3584 0.0
 124  52.2124 0.0
 125  48.6380 0.0
 126  52.2732 0.0
 127  49.1876 0.0
 128  46.6468 0.0
 129  44.1940 0.0
 130  54.4896 0.0
 131  50.9604 0.0
 132  47.8004 0.0
 133  51.4256 0.0
 134  50.9356 0.0
 135  63.1568 1.0
 136  61.4968 1.0
 137  59.7200 1.0
 138  57.9160 1.0
 139  59.1512 1.0
 140  56.5008 1.0
 141  61.6096 1.0
 142  63.7836 1.0
 143  56.7316 1.0
 144  61.4992 1.0
 145  59.1636 1.0
 146  63.5348 1.0
 147  55.3004 1.0
 148  51.8636 1.0
 149  58.2244 1.0
 150  64.9836 1.0
 151  64.9592 1.0
 152  54.9152 1.0
 153  60.8076 1.0
 154  60.1860 1.0
 155  54.8640 1.0
 156  57.5832 1.0
 157  62.8244 1.0
 158  67.1144 1.0
 159  62.4016 1.0
 160  67.0232 1.0
 161  58.0780 1.0
 162  61.1660 1.0
 163  54.1880 1.0
 164  59.3728 1.0
 165  55.7932 1.0
 166  59.8828 1.0
 167  57.6148 1.0
 168  62.1868 1.0
 169  60.5244 1.0
 170  59.2156 1.0
 171  59.8224 1.0
 172  51.1768 1.0
 173  59.5136 1.0
 174  62.7992 1.0
 175  61.2560 1.0
 176  62.2316 1.0
 177  64.9992 1.0
 178  64.6568 1.0
 179  66.2020 1.0
 180  61.3168 1.0
 181  61.7592 1.0
 182  59.7084 1.0
 183  66.9688 1.0
 184  61.6948 1.0
 185  53.5616 1.0
 186  54.9176 1.0
 187  64.1828 1.0
 188 . 1.0
 189 . 1.0
 190 . 1.0
 191 . 1.0
 192 . 1.0
 193 . 1.0
 194 . 1.0
 195 . 1.0
 196 . 1.0
 197 . 1.0
 198 . 1.0
 199 . 1.0
 ;

proc arima data = blood;
  identify var = demand crosscorr = step noprint;
  estimate input = step plot;
  forecast lead=12 out=fore;

run;

data combine;
  merge blood fore;

run;

goptions cback=white colors=(black) border reset=(axis symbol);

axis1 offset = (1 cm) label = ('time')
      order = (0 to 210 by 10);
axis2 label = (angle=90 'Blood Demand')
      order = (20 to 80 by 10);

symbol1 i=join l=1;
symbol2 h=1.5 pct v=F;
symbol3 l=20 i=join;

proc gplot data=combine;
  plot demand*t=1
       forecast*t=2
       l95*t = 3
       u95*t = 3 /overlay
                     haxis=axis1
                     vaxis=axis2
                     vminor=4
                     href = 135
                     lh=2;

run;

/* Let's figure out the systematic dynamics of the blood demand before the intervention
   t = 135.  */

data blooda;
  set blood;
  if _n_ < 135;

run;
 
proc arima data = blooda;
  identify var=demand;
  estimate p = 0 q = 0;
  estimate p = 1;
  estimate q = 1;
  estimate p = (2);
  estimate q = (2);

run;

/*  It appears that the MA model y(t) = phi0 + a(t) - thetha2*a(t-2) might be a good model
    for the systematic dynamics of the blood demand data before t = 135.  */

proc arima data = blood;
  identify var = demand crosscorr = step noprint;
  estimate input = step q=(2) plot;

  /*  These are two overfitting models.  */

  estimate input = ((1)step) q=(2) plot;
  estimate input = (/(1)step) q=(2) plot;
  
run;

/*  It looks like the best model for the blood demand data is a b = 0, r = 0, s = 0 
    intervention model with errors following a "gap" MA(2) model, namely,
    a(t) - thetha2*a(t-2). */

/*  So let's do our forecasting one last time with the slightly more sophisticated
    intervention model having the gap MA(2) error process.  */

proc arima data = blood;
  identify var = demand crosscorr = step noprint;
  estimate input = step q=(2) noprint;
  forecast lead = 12 out=fore;

data combine;
   merge blood fore;

run;

goptions cback=white colors=(black) border reset=(axis symbol);

axis1 offset = (1 cm) label = ('time')
      order = (0 to 210 by 10);
axis2 label = (angle=90 'Blood Demand')
      order = (20 to 80 by 10);

symbol1 i=join l=1;
symbol2 h=1.5 pct v=F;
symbol3 l=20 i=join;

proc gplot data=combine;
  plot demand*t=1
       forecast*t=2
       l95*t = 3
       u95*t = 3 /overlay
                     haxis=axis1
                     vaxis=axis2
                     vminor=4
                     href = 135
                     lh=2;

run;