questions over stochastic model

[chiade]

Dear Sir,

I am trying to utilize the "Replication file for Harvey, Ruiz and Shephard(1994), "Multivariate
Stochastic Variance Models", Review of Economic Studies, vol 61, no." I tested by inserting my numbers for a 3-variable model into the program u provided. I changed the dimensions to (3,3) and wonder wther it is correct?

I have some questions which I would like to clarify too. These questions were made in red. Many thanks for your kind replies.



* meanx2=mean of log chi-square,varx2=variance of log chi-square
*
compute meanx2=%digamma(0.5)-log(0.5)
compute varx2 =%trigamma(0.5)
*

dec vect[series] w(3)
set w(1) = log(svhat1)-meanx2
set w(2) = log(svhat2)-meanx2
set w(3) = log(svhat3)-meanx2
set w(4) = log(svhat4)-meanx2
*
The paper states that meanx2 = -1.27. The original program was written as w(1) = log(demean**2)-meanx2? since meanx2 is negative, is it subtracted to make make log e^2 positive?
in your program, is demean the standard deviation for dlogp?

I created a DLM model (not shown here), retrieving vhat and svhat . I subsequently used this svhat in place of demean**2 as seen above. However, there is no convergence. does it mean i am wrong and should use vhat instead?

By deriving the demean by:
diff(center) vhat1 / demean[/color][/color]
diff(center) vhat2 / demean[/color][/color]
diff(center) vhat3 / demean[/color][/color]
diff(center) vhat4 / demean[/color][/color][/color][/color]

dec packed rhosd(2,2) sigetaf(3,3)
dec symm sv(3,3)

What is the purpose of rhosd(2,2) and sigetaf(3,3)? does rhosd have to be one less than the size of sigetaf?

*
* %rhosw combines the packed lower triangular "rho" matrix with the
* known variances of the components to give the full covariance matrix.
* The sigma eta matrix is easier to estimate in square root form so it
* is also a packed lower triangle. This also makes the estimation
* process more similar to the factor method later.
*
function %rhosw
type symm %rhosw
local integer i j
dim %rhosw(3,3)
ewise %rhosw(i,j)=%if(i==j,varx2,varx2*rhosd(i-1,j))
end
*
function %ltouterxxF
type symm %ltouterxxF
local rect lt(3,3)
local integer i j
ewise lt(i,j)=%if(i>=j,sigetaf(i,j),0.0)
compute %ltouterxxF=lt*tr(lt)
end
*
* sigmaeps backs out the covariance matrix of epsilon from the estimate
* for eta
*
function %sigmaeps
type symm %sigmaeps
*
local integer i j n
local real rhowork rho factor
local parmset rootf
*
dim %sigmaeps(3,3)
do i=1,3
compute %sigmaeps(i,i)=%pi**2/2.0
do j=1,i-1
nonlin(parmset=rootf) rho
compute rho=rhosd(i-1,j)
find(noprint,parmset=rootf) root rhosd(i-1,j)-rhowork
compute rhowork=0.0
compute factor=2.0*rho**2
do n=1,100
compute rhowork=rhowork+factor/n
compute factor=factor*n/(n+.5)*rho**2
end do n
compute rhowork=rhowork*2.0/%pi**2
end find
compute %sigmaeps(i,j)=rho*%pi**2/2.0
end do j
end do i
end
*
*
nonlin rhosd sigetaf
compute [rect] c=%identity(3)
compute rhosd=%const(0.0)
compute sigetaf=%sqrt(%diag(||.001,.001,.001||))
why are these numbers (||4.2,16.1,3.4,19.4||)) chosen in your program?
dlm(y=%xt(w,t),c=c,startup=%(sv=%rhosw(),sw=.001*%ltouterxxF()),sw=sw,sv=sv,$
presample=diffuse,condition=1,method=bfgs,iters=100) 2 *
*
* Get the covariance matrix for epsilon
*
compute sigmaeps=%sigmaeps()
*
* Show the correlation matrices for epsilon and eta
*
disp %cvtocorr(sigmaeps)
disp %cvtocorr(sw)
*
* Table 4 analysis. Because the estimated sigmaeta is quite a bit
* different, these results are as well. The eigenvectors shown in this
* table are actually scaled by the square roots of eigenvalues as is
* done in table 5, so there is no need for a separate table.
*
eigen(scale) sw eigval eigvect
report(action=define)
report(atrow=1,atcol=1) "(a)" "Eigenvalues" eigval
report(atrow=2,atcol=2) "Eigenvectors" eigvect
report(atrow=6,atcol=2) "Percentage of Variance" eigval/%sum(eigval)
*
eigen(scale) %cvtocorr(sw) eigval eigvect
report(atrow=8,atcol=1) "(b)" "Eigenvalues" eigval
report(atrow=9,atcol=2) "Eigenvectors" eigvect
report(atrow=13,atcol=2) "Percentage of Variance" eigval/%sum(eigval)
report(action=show)
*
* Two factor model. This is parameterized somewhat differently than is
* shown in the text. It keeps four states, but with a deficient rank
* covariance matrix. This makes it a restriction on the original four
* state model.
*
dec rect theta(3,2)
I wonder whether changing to theta(3,2) is correct?
nonlin rhosd theta theta(1,2)=0.0
dlm(y=%xt(w,t),c=c,startup=%(sv=%rhosw(),sw=theta*tr(theta)),sv=sv,sw=sw,$
presample=diffuse,condition=1,type=smooth,$
pmethod=simplex,piter=10,method=bfgs) 2 * xstates
*
* Rotation to make 1st loading on yen=0
*
compute lambda1=%atan(-theta(2,1)/theta(2,2))
compute rotate1=||cos(lambda1),-sin(lambda1)|sin(lambda1),cos(lambda1)||
compute loading1=theta*rotate1
*
* Rotation to make 2nd loading on DM=0
*
compute lambda2=%atan(theta(1,1)/theta(2,1))
compute rotate2=||cos(lambda2),-sin(lambda2)|sin(lambda2),cos(lambda2)||
compute loading2=theta*rotate2
*
* Generate Table 6
*
report(action=define)
report(atrow=1,atcol=2,tocol=3,span) "Rotation 1"
report(atrow=1,atcol=4,tocol=5,span) "Rotation 2"
report(atrow=2,atcol=1,fillby=cols) "Pound/Dollar" "DM/Dollar" "Yen/Dollar" "SF/Dollar"
report(atrow=2,atcol=2) loading1
report(atrow=2,atcol=4) loading2
report(action=format,picture="*.###")
report(action=show)
*
* Invert the loadings to create the rotated factors
*
compute trans2=%ginv(loading2)
set factor21 2 * = %dot(%xrow(trans2,1),xstates(t))
set factor22 2 * = %dot(%xrow(trans2,2),xstates(t))
*
set dm1 2 * = factor21*loading2(2,1)
set yen2 2 * = factor22*loading2(3,2)
*
* Estimate the smoothed component from the univariate RW model for DM
*
nonlin swx
boxjenk(ar=1,ma=1,constant,noprint) w(2)
compute swx=.0161
dlm(y=w(2),c=1.0,sw=swx,sv=varx2,a=1.0,exact,$
method=bfgs,type=smooth,noprint) 2 * states
*
* The parameterization used for the factor model in the text generates
* factors which give roughly the correct level for the first two
* currencies. The other two require mean shifts, which is why there's a
* gap between the curves in 3(b). Since that decision was arbitrary (the
* two mean shifts could have been associated with any two currencies),
* we graph the standard deviations adjusted so the means match up (in
* log scale) with the means from the corresponding univariate process.
*
set dmuni 2 * = states(t)(1)
sstats(mean) 2 * dmuni-dm1>>adjust
set dm1 = exp(.5*(dm1+adjust))
set dmuni = exp(.5*dmuni)
*
graph(header="Figure 3(a) Smooth estimates of standard deviations for the first common trend") 2
# dm1
# dmuni
*
* Same for the yen
*
nonlin swx
boxjenk(ar=1,ma=1,constant,noprint) w(3)
compute swx=.0034
dlm(y=w(3),sw=swx,sv=varx2,c=1.0,a=1.0,exact,$
method=bfgs,type=smooth,noprint) 2 * states
*
set yenuni 2 * = states(t)(1)
sstats(mean) 2 * yenuni-yen2>>adjust
set yen2 = exp(.5*(yen2+adjust))
set yenuni = exp(.5*yenuni)
graph(header="Figure 3(b) Smooth estimates of standard deviations for the second common trend") 2
# yen2
# yenuni
*
* Estimation allowing for t-distribution. Add a diagonal adjustment
* matrix to the sv matrix. Since the elements of that are the variances
* of the log kappa's the adjustment has to be constrained to be
* non-negative.
*
dec vect zetad(3)
what is the purpose of zetad?
nonlin rhosd sigetaf zetad zetad>=0.0
compute [rect] c=%identity(3)
compute rhosd=%const(0.0)
compute sigetaf=%sqrt(%diag(||.001,.001,.001||))
dlm(y=%xt(w,t),c=c,startup=%(sv=%rhosw()+%diag(zetad),sw=.001*%ltouterxxF()),$
sw=sw,sv=sv,presample=diffuse,condition=1,$
pmethod=simplex,piter=0,method=bfgs,iters=100) 2 *

How do i get the covariance i.e. off -digaonal elements?

[TomDoan]

Dear Sir,

I am trying to utilize the "Replication file for Harvey, Ruiz and Shephard(1994), "Multivariate
Stochastic Variance Models", Review of Economic Studies, vol 61, no." I tested by inserting my numbers for a 3-variable model into the program u provided. I changed the dimensions to (3,3) and wonder wther it is correct?

I have some questions which I would like to clarify too. These questions were made in red. Many thanks for your kind replies.



* meanx2=mean of log chi-square,varx2=variance of log chi-square
*
compute meanx2=%digamma(0.5)-log(0.5)
compute varx2 =%trigamma(0.5)
*

dec vect[series] w(3)
set w(1) = log(svhat1)-meanx2
set w(2) = log(svhat2)-meanx2
set w(3) = log(svhat3)-meanx2
set w(4) = log(svhat4)-meanx2
*
The paper states that meanx2 = -1.27. The original program was written as w(1) = log(demean**2)-meanx2? since meanx2 is negative, is it subtracted to make make log e^2 positive?
in your program, is demean the standard deviation for dlogp?

If you're confused about this, you need to review the literature on (univariate) stochastic volatility models. Durbin and Koopman, for instance, has a section on it. meanx2=-1.27 is a two-digit rounding of the exact value, which is computed in the RATS program. You don't want log e^2 to be positive, you want it to be mean zero. The log of a chi-squared(1), which will be the error term in the actual (rather than approximate) measurement equation has a mean of roughly -1.27.

I created a DLM model (not shown here), retrieving vhat and svhat . I subsequently used this svhat in place of demean**2 as seen above. However, there is no convergence. does it mean i am wrong and should use vhat instead?

Don't do that. log(demean^2) is the dependent variable.

dec packed rhosd(2,2) sigetaf(3,3)
dec symm sv(3,3)

What is the purpose of rhosd(2,2) and sigetaf(3,3)? does rhosd have to be one less than the size of sigetaf?

See the top of page 255. The variances of the sigma-xi matrix are known, so all that needs to be estimated are the correlations (rhosd), which (as used in RATS) are the (N-1) x (N-1) lower submatrix of the correlation matrix.

compute sigetaf=%sqrt(%diag(||.001,.001,.001||))
why are these numbers (||4.2,16.1,3.4,19.4||)) chosen in your program?

I believe they're from univariate models.

dec rect theta(3,2)
I wonder whether changing to theta(3,2) is correct?

Yes, but a 2 factor model with 3 variables isn't especially interesting.

*
* Estimation allowing for t-distribution. Add a diagonal adjustment
* matrix to the sv matrix. Since the elements of that are the variances
* of the log kappa's the adjustment has to be constrained to be
* non-negative.
*
dec vect zetad(3)
what is the purpose of zetad?

That was just explained in the comments.

nonlin rhosd sigetaf zetad zetad>=0.0
compute [rect] c=%identity(3)
compute rhosd=%const(0.0)
compute sigetaf=%sqrt(%diag(||.001,.001,.001||))
dlm(y=%xt(w,t),c=c,startup=%(sv=%rhosw()+%diag(zetad),sw=.001*%ltouterxxF()),$
sw=sw,sv=sv,presample=diffuse,condition=1,$
pmethod=simplex,piter=0,method=bfgs,iters=100) 2 *

How do i get the covariance i.e. off -digaonal elements?

%RHOSW() is the first matrix in the sigma-xi on page 260, and SV is the sigma-eta matrix.