*
* Examples 17.6 from page 511, 17.7 from page 512
* Phillips-Perron tests
*
cal 1947 1 4
all 1989:1
open data gnptbill.txt
data(format=prn,org=obs)
*
* Direct calculations for 17.6
*
linreg tbill
# constant tbill{1}
*
* Doing things the hard way - computing the covariances and the weighted sums. Note that,
* because RATS uses series based at "1", the gamma-hat(0) is in entry 1 of covar, gamma-hat(1)
* is in entry 2, etc. METHOD=YULE is used to give the covariance calculations as shown in the
* text. The default method (METHOD=BURG) will give slightly different values.
*
corr(covariances,number=4,method=yule) %resids / covar
*
* Brute force calculation. Not recommended, as changing the number of lags requires an almost
* complete rewrite.
*
compute lambdasq=covar(1)+2*(4.0/5.0)*covar(2)+2*(3.0/5.0)*covar(3)+2*(2.0/5.0)*covar(4)+2*(1.0/5.0)*covar(5)
disp 'Method 1' lambdasq
*
* Doing this a bit more flexibly with a loop. Just be careful about the integer divide in
* doing the weights. The use of 2.0 forces a real divide, which is what you want.
*
compute lags=4
compute lambdasq=covar(1)
do i=2,lags+1
   compute lambdasq=lambdasq+2.0*(lags+2-i)/(lags+1)*covar(i)
end do i
disp 'Method 2' lambdasq
*
* And now, using the MCOV instead of the correlation and sums, since that's what
* MCOV is designed to do
*
mcov(lags=4,damp=1.0) / %resids
# constant
compute lambdasq=%cmom(1,1)/%nobs
disp 'Method 3' lambdasq
*
* However lambdasq is computed, the pieces needed for the final statistics are
*
*  %nobs = number of observations in the regression
*  %beta(2) = 2nd regression coefficient, that is, the one on tbill{1}
*  %xx(2,2) = 2,2 element of the regression inv(X'X) matrix, which is the same thing
*             as sigma-rho**2/s**2
*  %rss/%nobs = gamma-hat(0)
*
compute pprho=%nobs*(%beta(2)-1)-.5*(%nobs**2)*%xx(2,2)*(lambdasq-%rss/%nobs)
compute ppt  =sqrt((%rss/%nobs)/lambdasq)*(%beta(2)-1)/%stderrs(2)-$
   .5*(lambdasq-%rss/%nobs)*%nobs*sqrt(%xx(2,2))/sqrt(lambdasq)
disp 'PP rho statistic=' pprho
disp 'PP t statistic  =' ppt
*
* Using the PPUNIT procedure
*
@ppunit(nottest,lags=4) tbill
@ppunit(ttest,lags=4) tbill
*
* Example 17.7
*
set lgnp = 100.0*log(gnp)
set trend = t
*
* Direct calculation, using mcov to compute lambdasq
*
linreg lgnp
# constant lgnp{1} trend
mcov(lags=4,damp=1.0) / %resids
# constant
compute lambdasq=%cmom(1,1)/%nobs
compute pprho=%nobs*(%beta(2)-1)-.5*(%nobs**2)*%xx(2,2)*(lambdasq-%rss/%nobs)
compute ppt  =sqrt((%rss/%nobs)/lambdasq)*(%beta(2)-1)/%stderrs(2)-$
   .5*(lambdasq-%rss/%nobs)*%nobs*sqrt(%xx(2,2))/sqrt(lambdasq)
disp 'PP rho statistic=' pprho
disp 'PP t statistic  =' ppt
*
* Using PPUNIT. The results will be very slightly different as a result of using
* an asymptotically equivalent expression in the adjustment term.
*
@ppunit(trend,lags=4,nottest) lgnp
@ppunit(trend,lags=4,ttest) lgnp