/*
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 'Plot of nonstationary time series U';
title2 'Times Series is slow-turning';
title3 'Data needs to be differenced'; 
proc gplot data = okun;
   plot u*time;
run;

title 'Getting the autocorrelation function for U';
title2 'It appears to be slowly damping and thus U is nonstationary';
proc arima data = okun;
  identify var = u;
run;

title 'Difference in Unemployment versus Time';
title2 'Data is now in stationary form';
proc gplot data = okun;
   plot difu*time;
run;

title 'Autocorrelation Function of difu series';
title2 'The ACF of difu is quickly damping';
proc arima data = okun;
   identify var = difu;
run;

title 'Growth versus Time';
title2 'Data is quickly turning and therefore stationary';
proc gplot data = okun;
   plot g*time;
run;

title 'Autocorrelation Function of g series';
title2 'The ACF of g is quickly damping';
proc arima data = okun;
   identify var = g;
run;

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

/* The atheoretic approach:  Pure autoregression  */ 

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

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

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

/*  As it turns out, the AR(1) model is the best model.  */

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 'ARDL(p,q) model with p = ___ and q = ___';
proc autoreg data = okun_modified;
  model difu = / nlag = 4 method=ml backstep slstay=0.05;
run;

title 'ARDL(p,q) model with p = 1 and q = ___';
proc autoreg data = okun_modified;
  model difu = difu_1 / nlag = 4 method=ml backstep slstay=0.05;
run;

title 'ARDL(p,q) model with p = 2 and q = ___';
proc autoreg data = okun_modified;
  model difu = difu_1 difu_2 / nlag = 4 method=ml backstep slstay=0.05;
run;

title 'ARDL(p,q) model with p = ___ and q = 0';
proc autoreg data = okun_modified;
  model difu = g / nlag = 4 method=ml backstep slstay=0.05;
run;

title 'ARDL(p,q) model with p = ___ and q = 1';
proc autoreg data = okun_modified;
  model difu = g g_1 / nlag = 4 method=ml backstep slstay=0.05;
run;

title 'ARDL(p,q) model with p = ___ and q = 2';
proc autoreg data = okun_modified;
  model difu = g g_1 g_2/ nlag = 4 method=ml backstep slstay=0.05;
run;

title 'ARDL(p,q) model with p = 1 and q = 0';
proc autoreg data = okun_modified;
  model difu = difu_1 g / nlag = 4 method=ml backstep slstay=0.05;
run;

title 'ARDL(p,q) model with p = 2 and q = 0';
proc autoreg data = okun_modified;
  model difu = difu_1 difu_2 g / nlag = 4 method=ml backstep slstay=0.05;
run;

title 'ARDL(p,q) model with p = 1 and q = 1';
proc autoreg data = okun_modified;
  model difu = difu_1 g g_1 / nlag = 4 method=ml backstep slstay=0.05;
run;

title 'ARDL(p,q) model with p = 2 and q = 1';
proc autoreg data = okun_modified;
  model difu = difu_1 difu_2 g g_1 / nlag = 4 method=ml backstep slstay=0.05;
run;

title 'ARDL(p,q) model with p = 1 and q = 2';
proc autoreg data = okun_modified;
  model difu = difu_1 g g_1 g_2 / nlag = 4 method=ml backstep slstay=0.05;
run;

title 'ARDL(p,q) model with p = 2 and q = 2';
proc autoreg data = okun_modified;
  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.  We accept the null hypothesis
    of no structural break.  */

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

/* Here we need to determine the correct model B-J model for g.  */
/* As it turns out, the AR(2) (p=2) model is the best based on the
   AIC and SBC goodness-of-fit measures. */ 

title 'Estimation of Several ARMA models for g to help via ARDL to forecast u';
title2 'AR(1) Model for g';
proc arima data = okun;
  identify var = g;
  estimate p=1;
run;

title 'Estimation of Several ARMA models for g to help via ARDL to forecast u';
title2 'AR(2) Model for g';
proc arima data = okun;
  identify var = g;
  estimate p=2;
run;

title 'Estimation of Several ARMA models for g to help via ARDL to forecast u';
title2 'MA(1) Model for g';
proc arima data = okun;
  identify var = g;
  estimate q=1;
run;
 
title 'Estimation of Several ARMA models for g to help via ARDL to forecast u';
title2 'MA(2) Model for g';
proc arima data = okun;
  identify var = g;
  estimate q=2;
run;

title 'Estimation of Several ARMA models for g to help via ARDL to forecast u';
title2 'ARMA(1,1) Model for g';
proc arima data = okun;
  identify var = g;
  estimate p=1 q=1;
run;
  
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) =   */

/*  Now let us turn to the more direct approach with forecasting with a Transfer Function Model  */
 
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;