/*  Analysis of the Series M data set of Box-Jenkins using an 
    equal lag-length VAR.  */  
     

OPTIONS PAGESIZE =60 LINESIZE=72 NODATE;
data lead;
  input obs x y;
/* obs lead sales */
datalines;
  1 10.01 200.1
  2 10.07 199.5
  3 10.32 199.4
  4  9.75 198.9
  5 10.33 199.0
  6 10.13 200.2
  7 10.36 198.6
  8 10.32 200.0
  9 10.13 200.3
 10 10.16 201.2
 11 10.58 201.6
 12 10.62 201.5
 13 10.86 201.5
 14 11.20 203.5
 15 10.74 204.9
 16 10.56 207.1
 17 10.48 210.5
 18 10.77 210.5
 19 11.33 209.8
 20 10.96 208.8
 21 11.16 209.5
 22 11.70 213.2
 23 11.39 213.7
 24 11.42 215.1
 25 11.94 218.7
 26 11.24 219.8
 27 11.59 220.5
 28 10.96 223.8
 29 11.40 222.8
 30 11.02 223.8
 31 11.01 221.7
 32 11.23 222.3
 33 11.33 220.8
 34 10.83 219.4
 35 10.84 220.1
 36 11.14 220.6
 37 10.38 218.9
 38 10.90 217.8
 39 11.05 217.7
 40 11.11 215.0
 41 11.01 215.3
 42 11.22 215.9
 43 11.21 216.7
 44 11.91 216.7
 45 11.69 217.7
 46 10.93 218.7
 47 10.99 222.9
 48 11.01 224.9
 49 10.84 222.2
 50 10.76 220.7
 51 10.77 220.0
 52 10.88 218.7
 53 10.49 217.0
 54 10.50 215.9
 55 11.00 215.8
 56 10.98 214.1
 57 10.61 212.3
 58 10.48 213.9
 59 10.53 214.6
 60 11.07 213.6
 61 10.61 212.1
 62 10.86 211.4
 63 10.34 213.1
 64 10.78 212.9
 65 10.80 213.3
 66 10.33 211.5
 67 10.44 212.3
 68 10.50 213.0
 69 10.75 211.0
 70 10.40 210.7
 71 10.40 210.1
 72 10.34 211.4
 73 10.55 210.0
 74 10.46 209.7
 75 10.82 208.8
 76 10.91 208.8
 77 10.87 208.8
 78 10.67 210.6
 79 11.11 211.9
 80 10.88 212.8
 81 11.28 212.5
 82 11.27 214.8
 83 11.44 215.3
 84 11.52 217.5
 85 12.10 218.8
 86 11.83 220.7
 87 12.62 222.2
 88 12.41 226.7
 89 12.43 228.4
 90 12.73 233.2
 91 13.01 235.7
 92 12.74 237.1
 93 12.73 240.6
 94 12.76 243.8
 95 12.92 245.3
 96 12.64 246.0
 97 12.79 246.3
 98 13.05 247.7
 99 12.69 247.6
100 13.01 247.8
101 12.90 249.4
102 13.12 249.0
103 12.47 249.9
104 12.47 250.5
105 12.94 251.5
106 13.10 249.0
107 12.91 247.6
108 13.39 248.8
109 13.13 250.4
110 13.34 250.7
111 13.34 253.0
112 13.14 253.7
113 13.49 255.0
114 13.87 256.2
115 13.39 256.0
116 13.59 257.4
117 13.27 260.4
118 13.70 260.0
119 13.20 261.3
120 13.32 260.4
121 13.15 261.6
122 13.30 260.8
123 12.94 259.8
124 13.29 259.0
125 13.26 258.9
126 13.08 257.4
127 13.24 257.7
128 13.31 257.9
129 13.52 257.4
130 13.02 257.3
131 13.25 257.6
132 13.12 258.9
133 13.26 257.8
134 13.11 257.7
135 13.30 257.2
136 13.06 257.5
137 13.32 256.8
138 13.10 257.5
139 13.27 257.0
140 13.64 257.6
141 13.58 257.3
142 13.87 257.5
143 13.53 259.6
144 13.41 261.1
145 13.25 262.9
146 13.50 263.3
147 13.58 262.8
148 13.51 261.8
149 13.77 262.2
150 13.40 262.7
;

data yseries;
  set lead;
  y1=y;
  y2=x;

 /* Select Equal Lag-Length VAR using SBC criterion. Maximum Lag
    Length set at 6. */

proc varmax data=yseries(obs=120);
  model y1-y2 / dify(1) minic=(type=sbc p=(0:6) q=(0:0));
  run;
 /* Select Equal Lag-Length VAR using AICC criterion. Maximum Lag
    Length set at 10. */

proc varmax data=yseries(obs=120);
  model y1-y2 / dify(1) minic = (type=aicc p=(0:10) q=(0:0));

  run;

  /* As we can see the SBC and AICC criteria choose different lag
     lengths.  The SBC chooses lag length = 5 while the AICC chooses
     lag length = 8.  Let's go with lag = 8 since the residuals
     of such a VAR appear to be white noise. */ 

  /* Here we use Proc Varmax to conduct Granger Causal tests.
     Conclusion: There is one-way causality.  The leading indicator (y2)
     Granger-causes sales (y1)while sales do not Granger-cause the
     leading indicator. */ 

proc varmax data=yseries(obs=120);
  model y1-y2 / dify(1) p=8;
  causal group1=(y1) group2=(y2);
  causal group1=(y2) group2=(y1);

  run;
   
/*  Here we do the Granger Causal Tests "manually."  */
 
data yseries;
  set yseries;
  dy1=dif(y1);
  dy11=lag1(dy1);
  dy12=lag2(dy1);
  dy13=lag3(dy1);
  dy14=lag4(dy1);
  dy15=lag5(dy1);
  dy16=lag6(dy1);
  dy17=lag7(dy1);
  dy18=lag8(dy1);
  dy2=dif(y2);
  dy21=lag1(dy2);
  dy22=lag2(dy2);
  dy23=lag3(dy2);
  dy24=lag4(dy2);
  dy25=lag5(dy2);
  dy26=lag6(dy2);
  dy27=lag7(dy2);
  dy28=lag8(dy2);
 
  /* Test that y2 (i.e. x = leading indicator) Granger causes y1 (i.e. y = sales). */

proc autoreg data=yseries;
  model dy1 = dy11 dy12 dy13 dy14 dy15 dy16 dy17 dy18 
              dy21 dy22 dy23 dy24 dy25 dy26 dy27 dy28  / lagdep;
  test dy21, dy22, dy23, dy24, dy25, dy26, dy27, dy28;
title 'Test that y2 (leading indicator) Granger causes y1 (sales).';

  run;

  /* Test that y1 (i.e. y = sales) Granger causes y2 (i.e. x = leading indicator). */


proc autoreg data=yseries;
  model dy2 = dy11 dy12 dy13 dy14 dy15 dy16 dy17 dy18 
              dy21 dy22 dy23 dy24 dy25 dy26 dy27 dy28 / lagdep;
  test dy11, dy12, dy13, dy14, dy15, dy16, dy17, dy18; 
title 'Test that y1 (sales) Granger causes y2 (leading indicator).';

  run;

title ' ';
proc varmax data=yseries(obs=120);
  model y1-y2 / dify(1) p = 8;
  output out=result lead=1;

  run;

data result;
  set result;
  if _n_ = 121;
  keep for1 std1 for2 std2;

proc print data = result;

  run;

  /* Here we are using the diagnose option to print out the error
     cross-correlation matrices so that we can see if our lag length
     is sufficient to produce white noise in the residuals.  Lag 
     length 8 does better than lag length 5.  */ 
  
 proc varmax data=yseries(obs=120);
  model y1-y2 / dify(1) p=5 print = (diagnose);

  run;

proc varmax data=yseries(obs=120);
  model y1-y2 / dify(1) p=8 print = (diagnose);

  run;