Stochastic Volatility 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.

[bekkdcc]

Dear Tom, my prof insist of tt run the SV.rpf code. I can not able to decide where should I put the X term in the attached form. I looked for codes and also other pdf docs. I will be so pleased if you help me. I attached the doc. form an haighlighted my problems and questions.
Thanks in advance

[TomDoan]

I'm confused. What is X1? Is it a series with a unit coefficient? (Basically an input mean series for Y)? Is it a series with a non-zero coefficient that you forgot to include? Is it a coefficient (fixed but unknown mean of Y1)?

[bekkdcc]

for example Y is export rate and X is exchange rate (Basically an input mean series for Y)

[TomDoan]

Can't you just subtract x from y and model the difference?

[bekkdcc]

can I do that?
Because the model is like that:

Yt=b0+b1x1+ε_t √(h_t )

for example Y is export rate and X is exchange rate and b1 is the coefficient of the X1 series.

[TomDoan]

That's not what you wrote in the first place, and that isn't what you're describing. If it's an "input mean series" then yes, you can subtract it off. If the mean is an unknown linear function of X (which is what you just wrote), then you have to estimate the coefficients for that as part of the model, though you can do a two-step estimator where you use the residuals from a linear regression of y on 1 and x.

[bekkdcc]

Dear Tom,

How can I do that in the SV.rpf code? Can you describe me? Because ı couldnt find the right place inside the code and how ı should re-wrote the code “ to estimate the coefficient for that as part of..”
I will be so pleased if you could help me or show an example code.
Thanks in advance

[TomDoan]

I don't have an example of that. However, the standard procedure in the literature for estimating this using the state-space approximation is to remove the mean from the y series first and use log(e^2) or log(e^2) minus the theoretical mean of the error term as the Y component in the state-space model. So following that would mean regressing the series on 1 and X1 and using the residuals from that.

[bekkdcc]

Dear Tom,

I think I managed to change the code as I wanted.

I ran two different codes together, but could you check the logic accuracy? The first one is Commandeur & Koopman, An Introduction to State Space Time Series Analysis. * Chapter 5. UK data. and the second one is SV.RPF Estimation of a stochastic volatility model by (approximate) state-space methods

I am sending the code, the output and also my a doc. that explain my model equation and an explanation of the procedures I am trying to do in the code.

I will be so gratefull, if you can look and reply me whether its true or not, and what are the mistakes.

Thanks in advance

[TomDoan]

I'm not sure what the point is of using the C&K program for anything. (And particularly not taking your data and using C&K variable names---log ksi is log of Killed and Seriously Injured, which I doubt is what you're modeling.) Note also that you have some signs wrong in your description of U, and you aren't including a constant among the regressors, which is probably wrong.

Do NOT try to take your residuals and convert them into a form similar to the original data in SV.RPF. In SV.RPF, this removes the mean from the growth in the exchange rate, and the remainder of the program is done with DEMEAN. DEMEAN is the analogue to your residuals U. Just replace the top of what's in SV with your data calculations and LINREG to extract the residuals, and go from there. Don't try to turn this into two separate programs.


set dlogp = log(usxuk{0}/usxuk{1})
diff(center) dlogp / demean

[bekkdcc]

Dear Tom,

As you say, I completely canceled the first part.

1. did I use the linreg code in the right place or where should it be? :

2. and Can "diff(center) dlogp / demean" correctly express the error for the model that I am using or if not, how should I re-write it?

* SV.RPF
*
open data dk.xls
calendar(m) 2010
data(org=obs,format=xls) 2010:01 2018:12 usd ith
print
*
* 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)
*
* 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.
*

set dlogk = log(usd{0}/usd{1})
set dlogp = log(ith)
print
linreg dlogp
#constant dlogk{1} dlogp{1} dlogp{2} dlogp{3}
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,print) 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(method=bfgs,sw=sw,sv=varx2,y=ysq,type=filter,c=1.0, $
sx0=sw/(1-phi^2),x0=gammax,a=phi,z=gammax*(1-phi)) 2 * states
*
* Estimate the RW model
*
nonlin sw
dlm(method=bfgs,sw=sw,sv=varx2,y=ysq,type=filter,c=1.0, $
presample=diffuse,a=1.0,vhat=e,svhat=s) 2 * states

3. There is no problem in the results of "Estimate the unconstrained model" but there is an error code in the results of "Estimate the RW model" .It says : DLM2. No Observations Produce Valid Output. Check Data and Initial Values. What can I do for that?

[TomDoan]

No. "Demeaning" a series is equivalent to taking the residuals from a linear regression on a constant alone. Instead of demeaning the series, you just want the residuals from your LINREG. After your LINREG, delete the DIFF(CENTER) ... and replace it with

set demean = %resids

You can't start your estimation with entry 2, as you're losing entries due to the lags. Use the same range as the LINREG tells you it uses.

BTW, is

set dlogp = log(ith)

correct? That doesn't look like DLOG of anything.

[bekkdcc]

Dear Tom,

Thanks for your help. It seems perfect now, but I have some question for the last time

1. I corrected the set dlogp = log(ith) but what is the abbreviation BTW ?

2. There is a code to narrow down the sample range and analyze it, but I couldn't remember what was the code beginning with .....:

calendar(m) 2010
data(org=obs,format=xls) 2010:01 2018:12 usd ith
*
*
............ 2011:04 2015:09

4. phi seems statistically insignificant, is that mean SW model not appropriate or something else?

3. You can't start your estimation with entry 2, as you're losing entries due to the lags. Where is the entry 2 and how can change it?

Thanks in advance