* ============================================================
*                    COMPUTER LAB # 3
*               ECO 648, UNCG, Prof. Bearse
*
*  DATA: Download c.txt from our ECO 648 homepage
*        under miscellaneous.
* 
*  In this lab, we use RATS to estimate the autocorrelation
*  and partial autocorrelation functions of several
*  macroeconomic time series. 
*
*  As we have discussed in lecture, the sample autocorrelation
*  function (SACF) and the sample partial autocorrelation
*  function (SPACF) can be useful tools for deciding on an
*  appropriate ARMA(p,q) model of a macroeconomic time series.
*  In particular, we'd like to find an ARMA(p,q) model that whose
*  theoretical autocorrelation function (ACF) and theoretical
*  partial autocorrelation function (PACF) is similar to the
*  SACF and SPACF.
*
*               A QUICK SUMMARY OF THE ACF AND PACF
*                      OF AN ARMA(p,q)
*         --------------------------------------------------------                           
*         1) For a pure AR(p),
*              a) the autocorrelation function exhibits
*                 gradual decay (goes to zero only in the limit)
*              b) the partial autocorrelation function equals zero
*                 after lag p
*         2) For a pure MA(q),
*              a) the autocorrelation function equals zero after
*                 lag q
*              b) the partial autocorrelation function exhibits gradual
*                 decay (goes to zero only in the limit)
*         3) For ARMA(p,q) with p>0 and q>0, (THIS IS TOUGH!)
*              a) the autocorrelation function exhibits gradual decay
*                 after lag q (goes to zero only in the limit)
*              b) the partial autocorrelation function exhibits gradual
*                 decay after lag p (goes to zero only in the limit)
*               
* ============================================================
*
* STEP 1: Set Calendar and Load data 
*
CALENDAR 1959 1 12
ALLOCATE 40 1999:1
*
OPEN DATA c.txt
DATA(FORMAT=free,ORG=obs) / dates c
*            Monthly data on real personal consumption exp.
*            from 1959:1 through 1999:1. (Measured in
*            billions of Chained 1992 dollars. Seasonally
*            adjusted.)

set c = LOG(c)
*   We will work with the natural logarithm of
*   the series.

OPEN PLOT g1.rgf
GRAPH(header='Real Personal Consumption') 1
# c

STATISTICS c

*       It is good practice to examine a
*       plot of your data and compute summary
*       statistics of it.
*
*       The data exhibits a clear trend. Trending data
*       is not consistent with the assumption of a constant
*       mean. Consequently, the data on consumption does not appear
*       to be stationary.
*
*       The STATISTICS command created a variable
*                       "%NOBS"
*       which is the number of observations in
*       our data set.

* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
*             COMPUTING THE SACF AND SPACF
* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

COMPUTE numlags = FIX(%NOBS/4)
*       When computing the SACF and SPACF, we must decide how many
*       auto and partial correlations to estimate. The rule of
*       thumb is the sample size divided by 4.
*       The FIX() instruction removes the decimal
*       portion of %/NOBS/4 to ensure that the number
*       of lags estimated (numlags) is an integer.

CORRELATE(PRINT,NUMBER=numlags,STDERRS=sesacf,$
          PARTIAL=spacf,QSTATS) c / sacf
*
*    The CORRELATE statement is used to compute the SACF
*    and SPACF of our time series. We tell RATS to save
*       1) the SACF under the name sacf
*       2) the SPACF under the name spacf
*       3) the standard errors of the sacf
*          under the name sesacf
* 
*    Including the option QSTATS tells RATS to compute the
*    Ljung-Box Q statistics.
*    

* print(number=0) 1 numlags+1 sacf sesacf
* print(number=0) 1 numlags+1 spacf

* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
*         PLOTTING THE SACF AND SPACF ALONG WITH
*         TESTS OF SIGNIFICANCE
* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

SET sigsacf 1 numlags+1 = abs(sacf/sesacf)>1.96
*     
*     In large samples, the SACF is Normally distributed.
*     Thus, we reject
*             H0: sacf at lag j = 0
*     in favor of
*             H1: sacf at lag j not equal to 0
*     at the 0.05 level of significance if
*            abs(sacf/sesacf) > 1.96
*

SET sigspacf 1 numlags+1 = abs(sqrt(%NOBS)*spacf)>1.96
*
*     In large samples, the SPACF is Normally distributed
*     with standard error
*                     1/sqrt(T)
*     where T is the sample size.
*     Thus, we reject
*             H0: spacf at lag j = 0
*     in favor of
*             H1: spacf at lag j not equal to 0
*     at the 0.05 level of significance if
*            abs(spacf/sespacf) > 1.96
*       <==>
*            abs( spacf / (1/sqrt(T)) ) > 1.96
*       <==>
*            abs( sqrt(T)*spacf ) > 1.96
*
*     NOTE: In the CORRELATE statement, RATS defined a variable
*           %NOBS that was equal to the number of observations
*           of our data series.*

OPEN PLOT g2.rgf
GRAPH(HEADER='SACF',STYLE=BARGRAPH,NUMBER=1,MAX=1.0,$
      MIN=-1.0,SHADING=sigsacf) 1    
# sacf 2 numlags+1

open plot g3.rgf
GRAPH(HEADER='SPACF',STYLE=BARGRAPH,NUMBER=1,MAX=1.0,$
      MIN=-1.0,SHADING=sigspacf) 1
# spacf 2 numlags+1

*
*   The SACF we saw here is characteristic of
*   a nonstationary series. It decays toward
*   zero very gradually. This is consistent
*   with a theoretical ACF that does not decay
*   at all.

* ++++++++++++++++++++++++++++++++++++++++++++++
*           DETRENDED CONSUMPTION
* ++++++++++++++++++++++++++++++++++++++++++++++
*
*   Let's try to make the series stationary.
*   Based on our plot of the series, the 
*   nonstationarity appears to be due to
*   the trend in the data. Thus, let's try
*   to eliminate the trend.
*
*   One possibility is that consumption is
*   stationary around a deterministic linear
*   time trend;
*      i.e.,
*          c_t = constant + delta*t + ctilda_t
*   where
*          ctilda_t is a stationary series.
*   We can obtain the "detrended" series
*              ctilda
*   as the residuals from a linear regression
*   of c on a constant and time.

SET trend = t
*   This SET statement creates a linear time trend

LINREG c / ctilda
# CONSTANT trend
*                LINREG dep. variable / residual series
*                # list of predictors

OPEN PLOT g4.rgf
GRAPH(header='Detrended Real Personal Consumption') 1
# ctilda

COMPUTE numlags = FIX(%NOBS/4)
CORRELATE(PRINT,NUMBER=numlags,STDERRS=sesacf,$
          PARTIAL=spacf,QSTATS) ctilda / sacf
SET sigsacf 1 numlags+1 = abs(sacf/sesacf)>1.96
SET sigspacf 1 numlags+1 = abs(sqrt(%NOBS)*spacf)>1.96

OPEN PLOT g5.rgf
GRAPH(HEADER='SACF',STYLE=BARGRAPH,NUMBER=1,MAX=1.0,$
      MIN=-1.0,SHADING=sigsacf) 1    
# sacf 2 numlags+1

open plot g6.rgf
GRAPH(HEADER='SPACF',STYLE=BARGRAPH,NUMBER=1,MAX=1.0,$
      MIN=-1.0,SHADING=sigspacf) 1
# spacf 2 numlags+1

*
*   Notice that the SACF for ctilda also decays
*   very slowly. We could make the
*   argument that removing a linear time trend
*   did not make consumption stationary.
*

* ++++++++++++++++++++++++++++++++++++++++++++++
*           DIFFERENCED CONSUMPTION
* ++++++++++++++++++++++++++++++++++++++++++++++
*
*
*   One possibility is that consumption has
*   a stochastic trend caused by a unit root
*   in the AR polynomial. If this is the case,
*   then first differencing consumption should
*   make it stationary.
*      i.e.,
*          c_t = constant + c_(t-1) + eps_t
*          eps_t ~ WN(0,sigmasq)
*   implying
*          dc_t = c_t - c_(t-1)
*               = constant + eps_t
*
SET dc = c(t) - c(t-1)
*   This set statement creates a new series
*   dc which the first difference of c

OPEN PLOT g7.rgf
GRAPH(header='Differenced Real Personal Consumption') 1
# dc

COMPUTE numlags = FIX(%NOBS/4)
CORRELATE(PRINT,NUMBER=numlags,STDERRS=sesacf,$
          PARTIAL=spacf,QSTATS) dc / sacf
SET sigsacf 1 numlags+1 = abs(sacf/sesacf)>1.96
SET sigspacf 1 numlags+1 = abs(sqrt(%NOBS)*spacf)>1.96

OPEN PLOT g8.rgf
GRAPH(HEADER='SACF',STYLE=BARGRAPH,NUMBER=1,MAX=1.0,$
      MIN=-1.0,SHADING=sigsacf) 1    
# sacf 2 numlags+1

open plot g9.rgf
GRAPH(HEADER='SPACF',STYLE=BARGRAPH,NUMBER=1,MAX=1.0,$
      MIN=-1.0,SHADING=sigspacf) 1
# spacf 2 numlags+1

*
*  The SACF of differenced consumption decays quickly.
*  If we were to build an ARMA model, I'd probably
*  choose to build it for dc.
*

end