/*  rgdp = real gross domestic product, quarterly, seasonally adjusted in Billions
           of chained 2009 Dollars.  The series label is GDPC96.
    qual_spread = Quarterly average of the daily time series on Moody's Seasoned Baa
                  Corporate Cond Yield Relative to yield on 10-year Treasury Constant
                  Maturity.  This is what is referred to as the quality interest rate
                  spread.
    All of this data was downloaded from the FRED Economic Data base of the Federal
    Reserve Bank of St. Louis.  The shorter of the two time series is qual_spread with
    the first quarter observation being 1986 Q1.  
*/
 
data spread;
input date	yymmdd10. rgdp	qual_spread;
format date yymmdd10.;
datalines;
1986-01-01	7784.109	2.45
1986-04-01	7819.843	2.67
1986-07-01	7898.555	2.88
1986-10-01	7939.450	2.83
1987-01-01	7994.965	2.47
1987-04-01	8084.710	2.02
1987-07-01	8158.035	2.03
1987-10-01	8292.687	2.25
1988-01-01	8339.324	2.33
1988-04-01	8449.510	2.07
1988-07-01	8498.282	1.98
1988-10-01	8610.850	1.56
1989-01-01	8697.711	1.43
1989-04-01	8766.105	1.60
1989-07-01	8831.544	1.78
1989-10-01	8850.230	1.91
1990-01-01	8947.129	1.68
1990-04-01	8981.727	1.64
1990-07-01	8983.945	1.71
1990-10-01	8907.363	2.20
1991-01-01	8865.564	2.19
1991-04-01	8934.366	1.80
1991-07-01	8977.252	1.74
1991-10-01	9016.444	2.06
1992-01-01	9122.950	1.89
1992-04-01	9223.545	1.75
1992-07-01	9313.208	2.09
1992-10-01	9406.528	2.13
1993-01-01	9424.065	2.13
1993-04-01	9480.106	2.15
1993-07-01	9526.337	2.01
1993-10-01	9653.509	1.94
1994-01-01	9748.156	1.77
1994-04-01	9881.384	1.51
1994-07-01	9939.655	1.51
1994-10-01	10052.518	1.37
1995-01-01	10086.878	1.40
1995-04-01	10122.122	1.62
1995-07-01	10208.772	1.73
1995-10-01	10281.246	1.75
1996-01-01	10348.691	1.80
1996-04-01	10529.379	1.58
1996-07-01	10626.778	1.51
1996-10-01	10739.057	1.57
1997-01-01	10820.908	1.51
1997-04-01	10984.150	1.49
1997-07-01	11124.013	1.52
1997-10-01	11210.329	1.53
1998-01-01	11321.246	1.67
1998-04-01	11431.046	1.66
1998-07-01	11580.587	1.92
1998-10-01	11770.686	2.58
1999-01-01	11864.675	2.41
1999-04-01	11962.524	2.20
1999-07-01	12113.075	2.22
1999-10-01	12323.336	2.09
2000-01-01	12359.095	1.86
2000-04-01	12592.530	2.42
2000-07-01	12607.676	2.42
2000-10-01	12679.338	2.65
2001-01-01	12643.283	2.84
2001-04-01	12710.303	2.76
2001-07-01	12670.106	2.94
2001-10-01	12705.269	3.16
2002-01-01	12822.258	2.88
2002-04-01	12893.002	2.92
2002-07-01	12955.769	3.36
2002-10-01	12964.016	3.60
2003-01-01	13031.169	3.20
2003-04-01	13152.089	2.85
2003-07-01	13372.357	2.58
2003-10-01	13528.710	2.38
2004-01-01	13606.509	2.26
2004-04-01	13706.247	2.06
2004-07-01	13830.828	2.15
2004-10-01	13950.376	2.01
2005-01-01	14099.081	1.67
2005-04-01	14172.695	1.81
2005-07-01	14291.757	1.77
2005-10-01	14373.438	1.85
2006-01-01	14546.119	1.73
2006-04-01	14589.585	1.67
2006-07-01	14602.633	1.70
2006-10-01	14716.930	1.65
2007-01-01	14726.022	1.62
2007-04-01	14838.664	1.64
2007-07-01	14938.467	1.90
2007-10-01	14991.784	2.24
2008-01-01	14889.450	3.08
2008-04-01	14963.357	3.11
2008-07-01	14891.643	3.34
2008-10-01	14576.985	5.58
2009-01-01	14375.018	5.49
2009-04-01	14355.558	4.65
2009-07-01	14402.477	3.15
2009-10-01	14541.901	2.86
2010-01-01	14604.845	2.57
2010-04-01	14745.933	2.69
2010-07-01	14845.458	2.99
2010-10-01	14939.001	3.04
2011-01-01	14881.301	2.63
2011-04-01	14989.555	2.65
2011-07-01	15021.149	3.04
2011-10-01	15190.255	3.21
2012-01-01	15291.035	3.16
2012-04-01	15362.415	3.27
2012-07-01	15380.802	3.23
2012-10-01	15384.254	2.87
2013-01-01	15457.169	2.86
2013-04-01	15500.224	2.84
2013-07-01	15614.441	2.70
2013-10-01	15761.517	2.61
2014-01-01	15724.948	2.35
2014-04-01	15901.450	2.20
2014-07-01	16068.849	2.24
2014-10-01	16151.431	2.46
2015-01-01	16177.317	2.53
2015-04-01	16333.558	2.67
2015-07-01	16394.200	3.02
;

data spread;
  set spread;
  lrgdp = log(rgdp);
  g = (lrgdp - lag(lrgdp))*400;
  dqs = qual_spread - lag(qual_spread);
  time = _n_; 
  run;

symbol1 interpol=join
        value=dot;

title 'Plot of nonstationary time series rgdp';
title2 'Times Series is growing exponentially';
proc gplot data = spread;
   plot rgdp*time;
run;

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

title 'Plot of nonstationary time series lrgdp';
title2 'Times Series is growing linearly';
proc gplot data = spread;
   plot lrgdp*time;
run;

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

title 'Plot of g';
title2 'Times Series is flat and quickly turning';
proc gplot data = spread;
   plot g*time;
run;

title 'Getting the autocorrelation function for g';
title2 'It appears to be quicky damping and thus g is stationary';
proc arima data = spread;
  identify var = g;
run;

title 'Plot of nonstationary quality spread';
title2 'Times Series is slowly turning';
proc gplot data = spread;
   plot qual_spread*time;
run;

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

title 'Plot of difference of quality spread';
title2 'Times Series is flat and quickly turning';
proc gplot data = spread;
   plot dqs*time;
run;

title 'Getting the autocorrelation function for difference in quality spread';
title2 'It appears to be slowly damping and thus difference in quality spread is stationary';
proc arima data = spread;
  identify var = dqs;
run;

/*  Both the g and dqs series appear to be stationary.  */
    
/* Now we are going to use formal tests of whether a time series needs to be differenced
   to make it stationary.  The battery of tests are called the Augmented Dickey-Fuller
   Unit Root tests.  The major drawback of the SAS version of the ADF tests is that
   there is no choice of optimal number of augminting terms to use. */ 

title 'Unit Root Test for lrgdp';
proc arima data=spread;
    identify var=lrgdp stationarity=(adf=6);
run;

title 'Unit Root Test for g';
proc arima data=spread;
	identify var=g stationarity=(adf=6);
run;

title 'Unit Root Test for qual_spread';
proc arima data=spread;
    identify var=qual_spread stationary=(adf=6);
run;

title 'Unit Root Test for difference in quality spread';
proc arima data=spread;
    identify var=dqs stationary=(adf=6);
run;