/*  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.   */

data blood;
  input obs units interv;
  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
 ;

/*  Let's examine the systematic dynamics of the data
    before the intervention.  */

proc arima data = blood(obs=134);
  identify var = units;
  
run;

/*  The pre-intervention data seems to be white noise except
    for marginally significant peaks at lag = 2 in the
    ACF and PACF.  DeLurgio chooses to ignore these peaks
    and attributes them to sampling variation.  We will do
    the same.  */

/*  Prepare the extrapolation data set for forecasting
    with the intervention.  */

data add;
input obs units interv;
datalines;
  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
  ;
  
data extrap;
   set blood;

proc append base = extrap data = add;

run;

/* Ignoring the intervention in the data leads us to a
   Box-Jenkins model with substantial memory in it.  */

proc arima data=blood;
  identify var = units;
  run;

/*  Here we consider several intervention models for the data.  */

proc arima data = blood;
  identify var = interv noprint;
  identify var = units crosscorr = (interv) noprint;
  estimate input = (interv) method=ml;
  estimate input = ((1)interv) method=ml;
  estimate input = (/(1)interv) method=ml;
  estimate p=1 input = (interv) method=ml;

  run;


/* Forecasting with the intervention model */

  proc arima data = extrap;
  identify var = interv noprint;
  identify var = units crosscorr = (interv) noprint;
  estimate input = (interv) method=ml noprint;
  forecast lead = 12;
  title 'Forecasts produced by intervention model';

/* Forecasting with an incorrect AR(1) model */

proc arima data = blood;
  identify var = units noprint;
  estimate p = 1 noprint;
  forecast lead = 12;
  title 'Forecasts produced by an incorrect AR(1) model';

run;