/*
okun.def

 g  u

  Obs: 98, quarterly (1985Q2 - 2009Q3)  

	g = percentage change in U.S. Gross Domestic Product, seasonally adjusted.

        u = U.S. Civilian Unemployment Rate  (Seasonally adjusted)

	The variable DU used in Chapter 9 is defined as U(t)-U(t-1).

    Variable |       Obs        Mean    Std. Dev.       Min        Max
-------------+--------------------------------------------------------
           g |        98    1.276531    .6469279       -1.4        2.5
           u |        98    5.704082    1.132638        3.9        9.6


Data Source:  Federal Reserve Bank of St Louis
*/

data Okun;
input g	u;
datalines;
1.4	7.3
2	7.2
1.4	7
1.5	7
0.9	7.2
1.5	7
1.2	6.8
1.5	6.6
1.6	6.3
1.7	6
2.5	5.8
1.3	5.7
2.2	5.5
1.7	5.5
2.1	5.3
2.1	5.2
1.7	5.2
1.5	5.2
0.9	5.4
2.3	5.3
1.6	5.3
0.9	5.7
-0.1	6.1
0.6	6.6
1.4	6.8
1.2	6.9
1	7.1
1.6	7.4
1.7	7.6
1.5	7.6
1.6	7.4
0.8	7.1
1.2	7.1
1	6.8
1.9	6.6
1.5	6.6
1.9	6.2
1.2	6
1.6	5.6
0.8	5.5
0.7	5.7
1.3	5.7
1.2	5.6
1.3	5.5
2.1	5.5
1.2	5.3
1.7	5.3
1.4	5.2
1.7	5
1.6	4.9
1.1	4.7
1.1	4.6
1.1	4.4
1.7	4.5
2	4.4
1.3	4.3
1.1	4.3
1.6	4.2
2.2	4.1
1.1	4
2.5	3.9
0.7	4
1.1	3.9
0.3	4.2
1.3	4.4
0	4.8
0.7	5.5
1.2	5.7
1	5.8
0.9	5.7
0.6	5.9
1.1	5.9
1.1	6.1
2.2	6.1
1.4	5.8
1.6	5.7
1.6	5.6
1.5	5.4
1.6	5.4
1.9	5.3
1.1	5.1
1.8	5
1.4	4.9
2.1	4.7
1.2	4.7
0.8	4.6
1.2	4.4
1.4	4.5
1.5	4.5
1.3	4.7
1.1	4.8
0.3	4.9
0.9	5.4
0.3	6.1
-1.4	6.9
-1.2	8.1
-0.2	9.3
0.8	9.6
;

data okun;
  set okun;
  difu = u - lag(u);
  difu_1 = lag(difu);
  difu_2 = lag2(difu);
  g_1 = lag(g);
  g_2 = lag2(g);
  time = _n_;
run;

symbol1 interpol=join
        value=dot;

title 'Difference in Unemployment versus Time';
proc gplot data = okun;
   plot difu*time;
run;

title 'Growth versus Time';
proc gplot data = okun;
   plot g*time;
run;

/*  Are the difu and g series ergotic?  Let's look at the ACF's of the series
    to see if they are quickly dampling.  */

title 'Autocorrelation Function of difu series';
proc arima data = okun;
   identify var = difu;
run;

title 'Autocorrelation Function of g series';
proc arima data = okun;
   identify var = g;
run;

/*  Both series appear to be ergotic (stationary) because their ACF's are quickly damping.
    Let's proceed with our analysis. */

/* The atheoretic approach:  Pure autoregrssion  */ 

title 'Estimation of Three AR models';
proc arima data=okun;
  identify var=u(1);
    estimate p = 1;
	estimate p = 2;
	estimate p = 3;
run;

title 'Forecasting 4 steps ahead using AR(1) model';
proc arima data=okun;
  identify var=u(1) noprint;
    estimate p = 1 noprint;
	forecast lead = 4;
run;

data okun_modified;
  set okun;
  if _n_ > 3;
run;

/* Here we are going to pick a best ARDL(p,q) by choosing the one that
   has the smallest AIC and SBC measres while checking to make sure
   that the coefficients are statistically significant and that the 
   residuals of the model are white noise.  Notice that in the above data
   step we made sure that all of the competing models used the same number
   of observations, otherwise, the AIC and SBC measures would not be 
   comparable.  */

title 'Trying ARDL(p,q) models with p = 0,1,2 and q =  none,0,1,2';
title2 'The ARDL(1,1) model appears to be the best';
title3 'because it has the smallest AIC and SBC of the competing models';
proc autoreg data = okun_modified;
  model difu = / nlag = 4 method=ml backstep slstay=0.05;
  model difu = difu_1 / nlag = 4 method=ml backstep slstay=0.05;
  model difu = difu_1 difu_2 / nlag = 4 method=ml backstep slstay=0.05;
  model difu = g / nlag = 4 method=ml backstep slstay=0.05;
  model difu = g g_1 / nlag = 4 method=ml backstep slstay=0.05;
  model difu = g g_1 g_2/ nlag = 4 method=ml backstep slstay=0.05;
  model difu = difu_1 g / nlag = 4 method=ml backstep slstay=0.05;
  model difu = difu_1 difu_2 g / nlag = 4 method=ml backstep slstay=0.05;
  model difu = difu_1 g g_1 / nlag = 4 method=ml backstep slstay=0.05;
  model difu = difu_1 difu_2 g g_1 / nlag = 4 method=ml backstep slstay=0.05;
  model difu = difu_1 g g_1 g_2 / nlag = 4 method=ml backstep slstay=0.05;
  model difu = difu_1 difu_2 g g_1 g_2 / nlag = 4 method=ml backstep slstay=0.05;
run;

title 'ARDL(1,1) model: Using HAC standard errors';
title2 'See p. 357 in the Hill et al textbook';
title3 'for discussion of HAC standard errors';
proc model data = okun;
  endo difu;
  exog difu_1 g g_1;
  parms a0 d1 b0 d1;
  difu = a0 + d1*difu_1 + b0*g + b1*g_1;
  label a0 = "intercept";
  fit difu / gmm kernel = (bart,3,0);
run;

title 'Apply the Chow Test to see how stable the ARDL(1,1) equation is';
title2 'Model appears to be stable';
proc autoreg data = okun;
  model difu = difu_1 g g_1 / chow=(49);
run;

/*  The ARDL(1,1) model appears to be stable given the Chow test results. The Chow F-statistic
    is equal to 1.38 which has a p-value of 0.2469 which is greater than 0.05.  */

/*  Now we are going to pursue one of the original approaches to modeling time series
    regressions = the Cochrane-Orcutt transformation.  */

title 'AR(1) estimation of Okun equation of form';
title2 'difu = alpha + beta*growth + e with e = rho*e(-1) + v';
proc autoreg data=okun;
   model difu = g / nlag=1 dwprob;
run; 
 
title ' Estimation of the nonlinear version of the AR(1) Okun equation';
proc reg data = okun;
   model difu = difu_1 g g_1 / covb;
run;

/*  Here we conduct the Delta Method test of the nonlinear restriction
    implied by the AR(1) error specification in the Okun equation.
    Implied restriction: coef(difu_1)*coef(g) + coef(g_1) = 0 */

title 'Test Statistics for Common Factor Restriction Implied by AR(1) error specification';
title2 'Common Factor Restriction coef(difu_1)*coef(g) + coef(g_1) = 0 is NOT supported by the data.';
title3 'The Chi-square with one df is 35.19 with p-value = 3.0E-09.  Go with the ARDL(1,1) model.';
title4 'See page 215 - 220 in Hill et al book for discussion of the Delta method'; 
proc iml;

/* define necessary matrices for the delta method */

 /* v = cov(difu_1 g g_1) obtained from "covb" in above proc reg statement */

v = 
{0.0071525856 0.0007035581 0.0018709475, 
0.0007035581 0.0009423936 -0.000248285, 
0.0018709475 -0.000248285 0.0013560385};

g = {-0.18408,
     0.35012,
     1}; 

varest = g`*v*g;

seest = sqrt(varest);

est = 0.35012*(-0.18408) + -0.09916; *coef(difu_1)*coef(g) + coef(g_1) = 0;

z = est/seest;

chi_square = z**2;

print est;
print seest;
print z;
print chi_square;

end;

/*  Here we report the dynamic mulitplers of the ARDL(1,1) model.  See pp. 378 - 382 in the Hill et al book.  */

data multiplier;
  input lag int_mult cum_mult;
datalines;
  0 -0.1841 -0.1841
  1 -0.1636 -0.3477
  2 -0.0573 -0.4050
  3 -0.0201 -0.4251
  4 -0.0070 -0.4321 
  5 -0.0024 -0.4345
  6 -0.0009 -0.4354
  7 -0.0003 -0.4357
  ;

title 'The Interim (sometimes called impact and delay) multipliers for ARDL(1,1) Model';
proc gplot data = multiplier;
   plot int_mult*lag;
run;

title 'Cumulative Interim Multipliers for ARDL(1,1) Model';
proc gplot data = multiplier;
   plot cum_mult*lag;
run;

/*  Therefore we see that the cumulative multiplier is approaching the equilibrium multiplier value
    of (coef(g) + coef(g_1))/(1 - coef(difu_1)) = (-0.18408 -0.09916)/(1 - 0.35012) = -0.4358.  */

/*  Now we move to the issue of forecasting with the ARDL(1,1) model.  */

/*  First we need a formula to provide future values of the growth rate (g).  Let us assume that the g
    series follows an AR(2) process: g(t) = alpha + beta1*g(t-1) + beta2*g(t-2) + e(t)  */ 

title 'Getting forecasts of g to use in forecasting unemployment with the ARDL(1,1) model';
proc arima data = okun;
  identify var = g;
  estimate p = 2;
  forecast lead = 4;
run;

/* Using the future values of g we can produce the future changes in unemployment by recursively using the 
   the equation difu(t) =   */
 
title 'Using Proc ARIMA Transfer Function Capability to Forecast Unemployment using ARDL(1,1) model';
proc arima data= okun;
  identify var=g noprint;
  estimate p=2 method=ml;
  identify var=u(1)crosscor=(g) noprint;
  estimate input=((1)/(1)g) method=ml;
  forecast lead=4 nooutall;
run;