/* This program uses the data from
Enders, Walter and C.W.J. Granger.,
Unit-Root Tests and Asymmetric Adjustment with an Example Using the Term Structure of Interest Rates.
Journal of Business and Economic Statistics 16, July, 1998. pp. 304-11.

The technical details of the program are described on pages 128 - 132 of the RATS Progamming Manual.

Note that the consistent estimate of the threshold is obtained by Chan's method and differs from that
reported in the original paper. Results for Chan's method of estimating the threshold are also reported in:
Enders, Walter. "Applied Econometric Time Series", 2nd ed., (Wiley: Hoboken, NJ), 2004. pp. 429 - 34.
*/

cal 1958 1 4
all 1994:1
open data c:\granger.xls
data(format=xls,org=obs) /
table
source c:\winrats\bic.src      ; * As described on page 206 of my RATS Programming Manual, bic.src is a procedure to calculate the aic and sbc
source c:\winrats\bic1.src     ; * bic1 adds 1 regressor to account for the estimated threshold
                                * For example, the aic is calculated compute aic = (%nobs)*log(%rss) + 2*(%nreg+1)

set spread = r_10 - r_short
dif spread / ds
dif r_10 / drl
dif r_short / drs

gra(hea='Figure 1: The Interest Rates',patterns,key=below) 2 ; # r_10 ; # r_short

* This is the linear Dickey Fuller test reported in Table 5 of the paper
lin ds / resids; # constant spread{1} ds{1}
@bic

*** Now estimate the TAR model
linreg spread / y  ; # constant ; * demean the spread
set flag = %if(y{1}<0,0,1) ; * set the indicator function
set z_plus = flag*y{1} ; * create z_plus and z_minus
set z_minus = (1-flag)*y{1}
dif y / dy
lin dy 1958:3 *  ; * Estimate the TAR model reported in Table 5
# z_plus z_minus ; * To obtain the results for 1 lagged change include dy{1} in the regression @bic

exclude        ; * Perform the F-test for rho1 = rho2 = 0
# z_plus z_minus

* Now test rho1 = rho2
restrict 1
# 1 2
# 1 -1 0

*** Now estimate the M-TAR model reported in Table 5
set flag = %if(dy{1}<0,0,1)
set zplus = flag*y{1}
set zminus = (1-flag)*y{1}

linreg dy 58:3 *
# zplus zminus
* Obtain the aic, bic, and f-tests using the program steps above

***
* Now use the TAR-consistent specification
* Create the bounds for the lower and upper 15% bounds for the threshold
compute low = fix(%nobs*.15) , high = fix(%nobs*.85)
set thresh_test = spread
order thresh_test
compute thresh = thresh_test(low)

compute rss_test = 1000000.0
do i = low,high
   set yy = spread - thresh_test(i)
   set flag = %if(yy{1}<0,0,1)
   set minus = 1 - flag
   set z_plus = flag*yy{1}
   set z_minus = minus*yy{1}

   lin(noprint) ds
   # z_plus z_minus ds{1}

   if %rss < rss_test {
   	compute rss_test = %rss
   	compute thresh = thresh_test(i)
   }
end do i
dis 'We have found the attractor'
dis ' Threshold = ' thresh

   set flag = %if(spread{1}-thresh<0.,0,1)
   set minus = 1 - flag
   set z_plus = flag*(spread{1}-thresh)
   set z_minus = minus*(spread{1}-thresh)

   lin ds
   # z_plus z_minus ds{1}
   exclude        ; * Perform the F-test for rho1 = rho2 = 0
   # z_plus z_minus

   * Now test rho1 = rho2
   restrict 1
   # 1 2
   # 1 -1 0

* Now use the MTAR-consistent specification
compute rss_test = 1000000.0
set thresh_test = spread
order thresh_test
compute thresh = thresh_test(low)

do i = low,high
   set yy = spread - thresh_test(i)
   set flag = %if(ds{1}<0,0,1)
   set minus = 1 - flag
   set z_plus = flag*yy{1}
   set z_minus = minus*yy{1}

   lin(noprint) ds
   # z_plus z_minus ds{1}

   if %rss < rss_test {
   	compute rss_test = %rss
   	compute thresh = thresh_test(i)
   }
end do i
dis 'We have found the attractor'
dis ' Threshold = ' thresh

   set flag = %if(ds{1}<0.,0,1)
   set minus = 1 - flag
   set z_plus = flag*(spread{1}-thresh)
   set z_minus = minus*(spread{1}-thresh)

   lin ds
   # z_plus z_minus ds{1}
   exclude        ; * Perform the F-test for rho1 = rho2 = 0
   # z_plus z_minus

   * Now test rho1 = rho2
   restrict 1
   # 1 2
   # 1 -1 0


**** Now use OLS to obtain the residuals from the equilibrium relationship
lin spread / resids ; # constant
* Estimate the error-correction model reported as the second set of equations in Table 6
system 1 to 2
var drl drs
lags 1 to 2
det constant resids{1}
end(system)
estimate(outsigma=v)

*Calculate the multivariate aic and sbc reported in Table 6
compute aic = %nobs * %logdet + 2*12
compute sbc = %nobs * %logdet + 12*log(%nobs)
display 'aic = ' aic 'sbc = ' sbc

** Instead use the MTAR values of z_plus and z_minus
* Estimate the error-correction model reported as the second set of equations in Table 6
system 1 to 2
var drl drs
lags 1 to 2
det z_plus z_minus
end(system)
estimate(outsigma=v)

*Calculate the multivariate aic and sbc reported in Table 6
compute aic = %nobs * %logdet + 2*14
compute sbc = %nobs * %logdet + 14*log(%nobs)
display 'aic = ' aic 'sbc = ' sbc