*
* Replication file for Harvey, Ruiz and Shephard(1994), "Multivariate
* Stochastic Variance Models", Review of Economic Studies, vol 61, no.
* 2, pp 247-264.
*
* Univariate Models (results from table 3)
*
open data xrates.xls
data(format=xls,org=columns) 1 946 usxuk usxger usxjpn usxsui
*
* meanx2=mean of log chi-square,varx2=variance of log chi-square.
* These are exact (rather than rounded) values.
*
compute meanx2=%digamma(0.5)-log(0.5)
compute varx2 =%trigamma(0.5)
*
* Table 1
*
report(action=define)
report(atrow=1,atcol=2) "Dlogp" "(Dlogp)^2" "Log Dlogp^2"
dofor p = usxuk usxger usxjpn usxsui
   report(row=new,atcol=1) %l(p)
   set dlogp = log(p{0}/p{1})
   diff(center) dlogp
   set dlogp2 = dlogp^2
   set logdlogp2 = log(dlogp2)
   corr(number=10,method=yule,qstats,noprint) dlogp
   report(row=current,atcol=2) %qstat
   corr(number=10,method=yule,qstats,noprint) dlogp2
   report(row=current,atcol=3) %qstat
   corr(number=10,method=yule,qstats,noprint) logdlogp2
   report(row=current,atcol=4) %qstat
end dofor
report(action=format,picture="*.##")
report(action=show)
*
*  Table 2
*
* The <<gamma>> and <<phi>> variables are very highly correlated when
* <<phi>> is near 1, which makes estimation by general hill-climbing
* procedures somewhat tricky. Instead, we reparameterize this to use
* gammax=gamma(1-phi) in place of <<gamma>>. Note that the log
* likelihoods in the article differ from the ones shown here by the
* constant -T/2 log(2pi)
*
report(action=define)
report(atcol=1,atrow=2) "(a)"
report(atcol=1,atrow=8) "(b)"
report(atcol=2,atrow=2,fillby=cols) "phi" "" "sigma" "" "gammax" "log L" "sigma" "" "log L" "Q(10)"
compute col=3
dofor p = usxuk usxger usxjpn usxsui
   set dlogp = log(p{0}/p{1})
   diff(center) dlogp / demean
   nonlin phi sw gammax
   set ysq = log(demean^2)-meanx2
   *
   * Get initial guess values from ARMA(1,1) model
   *
   boxjenk(ar=1,ma=1,constant,noprint) ysq
   *
   * phi is taken directly as the AR(1) coefficient. The variance sw is
   * backed out by matching first order autocorrelations in the MA term.
   * gamma is chosen to reproduce the mean of the ysq series
   *
   compute phi=%beta(2),sw=-phi*varx2*(1+%beta(3)^2)/%beta(3)-(1+phi^2)*varx2
   compute sw=%if(sw<0,.1,sw)
   compute gammax=%mean
   *
   * Estimate the unconstrained model
   *
   dlm(y=ysq,sw=sw,sv=varx2,c=1.0,a=phi,z=gammax*(1-phi),presample=ergodic,$
      method=bfgs,type=filter) 2 * states
   report(atcol=col,atrow=2,fillby=cols) phi %stderrs(1) sw %stderrs(2) gammax %logL
   *
   * Estimate the RW model
   *
   nonlin sw
   dlm(y=ysq,sw=sw,sv=varx2,c=1.0,a=1.0,exact, $
       method=bfgs,type=filter,vhat=e,svhat=s) 2 * states
   report(atcol=col,atrow=8,fillby=cols) sw %stderrs(1) %logL
   *
   * Compute the LB Q statistic on the standardized residuals
   *
   set estd 2 * = e(t)(1)/sqrt(s(t)(1,1))
   corr(number=10,method=yule,qstats,noprint) estd
   report(atcol=col,atrow=11) %qstat
   compute col=col+1
end dofor p
report(action=format,align=decimal,picture="*.####")
report(action=format,atrow=3,torow=3,special=parens)
report(action=format,atrow=5,torow=5,special=parens)
report(action=format,atrow=9,torow=9,special=parens)
report(action=format,atrow=7,torow=7,picture="*.##")
report(action=format,atrow=10,torow=11,picture="*.##")
report(action=show)
*
* Table 3
* These are actually done with lags=8 the way they are counted by
* DFUNIT. (LAGS=n on DFUNIT means n lags on the difference, which would
* mean an n+1 lag autoregression).
*
report(action=define)
dofor p = usxuk usxger usxjpn usxsui
   set dlogp = log(p{0}/p{1})
   diff(center) dlogp / demean
   set ysq = log(demean^2)
   @dfunit(lags=8,noprint) ysq
   report(col=new,atrow=1,fillby=cols) %l(p) %cdstat
end dofor
report(action=format,picture="*.##")
report(action=show)
*
* Figure 1. Re-estimate the model for the US/UK data.
* Add the TYPE=SMOOTH option to get the Kalman smoothed states
*
set dlogp = log(usxuk{0}/usxuk{1})
diff(center) dlogp / demean
nonlin phi sw gammax
set ysq = log(demean^2)-meanx2
boxjenk(ar=1,ma=1,constant,noprint) ysq
compute phi=%beta(2),sw=-phi*varx2*(1+%beta(3)^2)/%beta(3)-(1+phi^2)*varx2
compute sw=%if(sw<0,.1,sw)
compute gammax=%mean
*
dlm(y=ysq,sw=sw,sv=varx2,c=1.0,a=phi,z=gammax*(1-phi),presample=ergodic,$
   method=bfgs,type=smooth) 2 * states
*
set absval   = abs(dlogp)
set smoothed = exp(.5*states(t)(1))
graph(footer="Figure 1") 2
# absval
# smoothed