*
* Hamp582.prg
* Cointegration analysis of exchange rate data.
* pp 582-599
*
cal 1973 1 12
open data exdata.rat
data(format=rats) 1973:1 1989:10 pc6it pzunew exritl
*
* Transform and normalize the data series
*
set italcpi = 100*log(pc6it/pc6it(1973:1))
set uscpi   = 100*log(pzunew/pzunew(1973:1))
set exrat   = -100*log(exritl/exritl(1973:1))
*
graph(header='Figure 19.2',key=attached,klabel=||'Italy CPI','US CPI','Ex Rate'||) 3
# italcpi
# uscpi
# exrat
*
* Dickey-Fuller tests on the variables
*
@dfunit(lags=12,trend) uscpi
@dfunit(lags=12,trend) italcpi
@dfunit(lags=12,trend) exrat
*
* Unit root tests on the hypothesized cointegrating vector
*
set ppp = uscpi-exrat-italcpi
graph(header='Figure 19.3 The real dollar-lira exchange rate')
# ppp
@dfunit(lags=12) ppp
@ppunit(lags=12) ppp
*
* IRF for ppp deviation.
*  1. Do a 12 lag autoregression and save the equation.
*  2. Do an IMPULSE instruction. We use the INPUT option because
*     we want a first period shock of 1.0, not one standard error.
*  3. Graph, use the number=0 option to label the entries as 0,...,71
*
linreg(define=pppeq) ppp
# constant ppp{1 to 12}
impulse(input,noprint) 1 72
# pppeq pppimp
# 1.0
graph(header='Figure 19.4. Impulse Response Function for Real Dollar-Lira Exchange Rate',$
  number=0)
# pppimp
*
* Residual-based tests using an estimated cointegrating vector
* pp 598-599
*
linreg uscpi / coresids
# constant exrat italcpi
linreg coresids / resids
# coresids{1}
mcov(lags=12,damp=1.0) / resids
# constant
*
* The mcov produces a raw sum, which needs to be divided by the number of observations
* to get the required "lambda**2" value. c0 can be obtained by dividing the sum of
* squared residuals by the same number of observations. The sigma for rho and the standard
* error of estimate can be pulled out of the regression variables.
*
* Note that the %nobs variable will be equal to the number of observations in the AR(1)
* regression, which will be T-1 in Hamilton's notation.
*
compute lambdasq=%cmom(1,1)/%nobs
compute c0      =%rss/%nobs
compute sigrho  =%stderrs(1)
compute sig     =sqrt(%seesq)
*
compute zp=%nobs*(%beta(1)-1)-.5*(%nobs*sigrho/sig)**2*(lambdasq-c0)
compute zt=sqrt(c0/lambdasq)*(%beta(1)-1)/sigrho-.5*(%nobs*sigrho/sig)*(lambdasq-c0)/sqrt(lambdasq)
disp 'Phillips-Ouliaris Tests'
disp 'Zp=' zp
disp 'Zt=' zt