*
* Example of Gibbs sampling with DCC model. This uses the same data set
* as the GARCHMV.RPF example.
*
open data g10xrate.xls
data(format=xls,org=columns) 1 6237 usxjpn usxfra usxsui
*
set xjpn = 100.0*log(usxjpn/usxjpn{1})
set xfra = 100.0*log(usxfra/usxfra{1})
set xsui = 100.0*log(usxsui/usxsui{1})
compute n=3
dec vect[strings] vlabels(n)
compute vlabels=||"Japan","France","Switzerland"||
***********************************************************************
garch(p=1,q=1,mv=dcc,rvectors=rv,hmatrices=h) / xjpn xfra xsui
compute gstart=%regstart(),gend=%regend()
*
* Proposal density for the increment in independence chain MH is based
* upon the covariance matrix from the GARCH estimates. This is a
* multivariate t with <> degrees of freedom. The scale on FXX and
* the value for NUXX may need to be adjusted in different applications
* to get a reasonable acceptance rate.
*
compute fxx=.5*%decomp(%xx)
compute nuxx=10.0
*
* Start at the ML estimates
*
compute logplast=%logl
compute blast=%beta
compute accept=0
*
compute nburn=1000,nkeep=1000
*
* This keeps track of the DCC correlations estimates. This is a
* complicated structure since we have, for each draw, one full series
* for each off-diagonal. (For simplicity, this is also saving the
* diagonals, which are, of course, 1's). With a bigger model (longer
* time series or more variables) or more draws, you could run into
* memory issues.
*
dec symm[series[vect]] dcccorr(n,n)
do i=1,n
do j=1,i
gset dcccorr(i,j) gstart gend = %zeros(nkeep,1)
end do j
end do i
*
dec series[symm] hlast
gset hlast gstart gend = h
*
infobox(action=define,progress,lower=-nburn,upper=nkeep) "Independence M-H"
do draw=-nburn,nkeep
compute btest=blast+%ranmvt(fxx,nuxx)
*
* Evaluate the GARCH model at the proposal
*
garch(p=1,q=1,mv=dcc,initial=btest,method=evaluate,rvectors=rv,hmatrices=htest) / xjpn xfra xsui
compute logptest=%logl
*
* Carry out the Metropolis test
*
compute alpha=exp(logptest-logplast)
if alpha>1.0.or.%uniform(0.0,1.0)0 {
do i=1,n
do j=1,i
do t=gstart,gend
compute dcccorr(i,j)(t)(draw)=%cvtocorr(hlast(t))(i,j)
end do t
end do j
end do i
}
end do draw
infobox(action=remove)
*
dec vect xcorr(nkeep)
dec symm[series] lower(n,n) upper(n,n) median(n,n)
clear(zeros) lower upper median
*
* Generate the median and 90% confidence band for the DCC correlations.
* For this data set, the period to period movement is large enough
* relative to the spread that the three lines basically end up on top of
* each other. If you graph over a shorter period, you can see at least
* some spread between the upper and lower bounds.
*
do j=2,n
do k=1,j-1
do t=gstart,gend
ewise xcorr(i)=dcccorr(j,k)(t)(i)
compute frac=%fractiles(xcorr,||.05,.50,.95||)
compute lower(j,k)(t)=frac(1)
compute upper(j,k)(t)=frac(3)
compute median(j,k)(t)=frac(2)
end do tme
graph(header="DCC Estimates for "+vlabels(k)+" with "+vlabels(j)) 3
# median(j,k) gstart gend
# lower(j,k) gstart gend 2
# upper(j,k) gstart gend 2
end do k
end do j