Stochastic Volatility Models:WINBUGS and R

This mourning I am editing my working paper, "Testing the Commodity Price Leverage Effect through Non-Linear State Space Stochastic Volatility Models" for publication and want to mention my use of WINBUGS 1.4 (Windows Bayesian Inference Using Gibbs Sampling). This package enables me to estimate stochastic volatility models as Yu(2002, 2005) did and I can embed the application in R. The software can be downloaded at
http://www.mrc/-bsu.cam.ac.uk/bugs/ .

An example of a stochastic volatility in WinBUGS and R is posted on my older notes site at
http://www.thecromwellworkshop.com/NotesFromWorkshop/BayesianStochasticVolatilityModels.aspx. However, I am in the process of migrating these older posts and reproduce the information here.

A: The Stochastic Model for WinBugs

model
{
mu~dnorm(0,0.04)
phistar~dbeta(20,1.5)
itau2~dgamma(2.5,0.025)
rho~dunif(-1,1)
#beta<-exp(mu/2)
phi<-2*phistar-1
pi <- 3.141592654
tau<-sqrt(1/itau2)
alpha<-mu*(1-phi)
theta0 ~ dnorm(mu,itau2)
thmean[1] <- mu + phi*(theta0-mu)
theta[1] ~ dnorm(thmean[1],itau2)I(-100,100)
for (i in 2:N) {
thmean[i] <- mu + phi*(theta[i-1]-mu)
theta[i] ~ dnorm(thmean[i],itau2)I(-80,80)
}
Ymean[1]<-0
Yisigma2[1]<-1/(exp(theta[1]))
Y[1]~dnorm(Ymean[1],Yisigma2[1])
loglike[1]<- (-0.5*log(2*pi) + 0.5*log(Yisigma2[1]) - 0.5*Yisigma2[1]*Y[1]*Y[1]) for (i in 2:N) {
Ymean[i]<-rho/tau*exp(0.5*theta[i])*(theta[i]-mu-phi*(theta[i-1]-mu))
Yisigma2[i] <- 1/(exp(theta[i])*(1-rho*rho))
Y[i]~ dnorm(Ymean[i],Yisigma2[i])
loglike[i]<- (-0.5*log(2*pi) + 0.5*log(Yisigma2[i]) - 0.5*Yisigma2[i]*Y[i]*Y[i])
}
node1<- -sum(loglike[])
}

B: Running the above model (A) in R with BRUGS.

Here is the SVModel in SVModel.txt. Just copy and paste into the R GUI at the prompt. Make sure that you download the libraries and the 10,000 iterations takes some time so you might want to initially reduce the amount when testing the algorithm.

library(BRugs)
#
# Check the syntax of the model
#
modelCheck("SV1.txt")
#
# Load the Data
#
modelData("dSugar.txt")
#
# Specify the number of chains
#
modelCompile(numChains=1)
#
# Specify the initial conditions for the chains
#
#modelInits("SVinit2.txt")
modelGenInits()
#
# Initial burin-in
#
modelUpdate(1000)
#
# Monitor the variables
#
samplesSet(c("phi","mu","itau2","rho"))
#
# Do more iterations
#
modelUpdate(10000)
#
# variable statistics
#
samplesStats("*")
#
#Deviance Information Criterior
#
dicSet()
dicStats()
#
# Build MCMC objects for CODA
#
MC.phi<-buildMCMC("phi")
MC.itau2<-buildMCMC("itau2")
MC.rho<-buildMCMC("rho")
MC.mu<-buildMCMC("mu")

#
# Raferty and Lewis diagnostics
#
raftery.diag(MC.phistar, q=0.025, r=0.005, s=0.95, converge.eps=0.001)
raftery.diag(MC.itau2, q=0.025, r=0.005, s=0.95, converge.eps=0.001)
raftery.diag(MC.rho, q=0.025, r=0.005, s=0.95, converge.eps=0.001)
raftery.diag(MC.mu, q=0.025, r=0.005, s=0.95, converge.eps=0.001)
#
# Gewek Convergence diagnostics
#
geweke.diag(MC.phistar, frac1=0.1, frac2=0.5)
geweke.diag(MC.itau2, frac1=0.1, frac2=0.5)
geweke.diag(MC.rho, frac1=0.1, frac2=0.5)
geweke.diag(MC.mu, frac1=0.1, frac2=0.5)
#


C: Here is the output in R.

R version 2.4.0 (2006-10-03)
Copyright (C) 2006 The R Foundation for Statistical Computing
ISBN 3-900051-07-0

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
Natural language support but running in an English locale
R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(BRugs)
Welcome to BRugs running on OpenBUGS version 2.2.0 beta
> #
> # Check the syntax of the model
> #
> modelCheck("SV1.txt")
model is syntactically correct
> #
> # Load the Data
> #
> modelData("dSugar.txt")
data loaded
> #
> # Specify the number of chains
> #
> modelCompile(numChains=1)
model compiled
> #
> # Specify the initial conditions for the chains
> #
> #modelInits("SVinit2.txt")
> modelGenInits()
initial values generated, model initialized
> #
> # Initial burin-in
> #
> modelUpdate(1000)
1000 updates took 212 s
> #
> # Monitor the variables
> #
> samplesSet(c("phi","mu","itau2","rho"))
monitor set for variable 'phi'
monitor set for variable 'mu'
monitor set for variable 'itau2'
monitor set for variable 'rho'
> #
> # Do more iterations
> #
> modelUpdate(10000)
10000 updates took 1949 s
> #
> # variable statistics
> #
> samplesStats("*")
mean sd MC_error val2.5pc median val97.5pc start sample
itau2 3.5300 0.70520 0.062810 2.4280 3.4170 5.03200 1001 10000
mu -8.1460 0.06697 0.002150 -8.2740 -8.1470 -8.01000 1001 10000
phi 0.8258 0.03279 0.002703 0.7589 0.8264 0.88310 1001 10000
rho -0.1139 0.03989 0.001499 -0.1907 -0.1138 -0.03554 1001 10000
> #
> #Deviance Information Criterior
> #
> dicSet()
deviance set
> dicStats()
Dbar Dhat DIC pD
Y 0 0 0 0
theta 0 0 0 0
total 0 0 0 0
> #
> # Build MCMC objects for CODA
> #
> MC.phi<-buildMCMC("phi")
Loading required package: coda
Loading required package: lattice
Warning message:
package 'lattice' was built under R version 2.4.1
> MC.itau2<-buildMCMC("itau2")
> MC.rho<-buildMCMC("rho")
> MC.mu<-buildMCMC("mu")
>
> #
> # Raferty and Lewis diagnostics
> #
> raftery.diag(MC.phistar, q=0.025, r=0.005, s=0.95, converge.eps=0.001)
Error in inherits(x, "mcmc.list") : object "MC.phistar" not found
> raftery.diag(MC.itau2, q=0.025, r=0.005, s=0.95, converge.eps=0.001)

Quantile (q) = 0.025
Accuracy (r) = +/- 0.005
Probability (s) = 0.95

Burn-in Total Lower bound Dependence
(M) (N) (Nmin) factor (I)
itau2 144 160832 3746 42.9
> raftery.diag(MC.rho, q=0.025, r=0.005, s=0.95, converge.eps=0.001)

Quantile (q) = 0.025
Accuracy (r) = +/- 0.005
Probability (s) = 0.95

Burn-in Total Lower bound Dependence
(M) (N) (Nmin) factor (I)
rho 14 16032 3746 4.28


> raftery.diag(MC.mu, q=0.025, r=0.005, s=0.95, converge.eps=0.001)

Quantile (q) = 0.025
Accuracy (r) = +/- 0.005
Probability (s) = 0.95

Burn-in Total Lower bound Dependence
(M) (N) (Nmin) factor (I)
mu 3 4095 3746 1.09
> #
> # Gewek Convergence diagnostics
> #
> geweke.diag(MC.phistar, frac1=0.1, frac2=0.5)
Error in inherits(x, "mcmc.list") : object "MC.phistar" not found
> geweke.diag(MC.itau2, frac1=0.1, frac2=0.5)
[[1]]

Fraction in 1st window = 0.1
Fraction in 2nd window = 0.5

itau2
-3.931
> geweke.diag(MC.rho, frac1=0.1, frac2=0.5)
[[1]]

Fraction in 1st window = 0.1
Fraction in 2nd window = 0.5

rho
2.779

> geweke.diag(MC.mu, frac1=0.1, frac2=0.5)
Fraction in 1st window = 0.1
Fraction in 2nd window = 0.5
mu
-3.410

As I build the web services API to re-run the empirical section for my article mentioned, I plan to combine the two into a tutorial post for http://www.codeproject.com/ and discuss some of these features at length. I like the approach by http://www.lokad.com/ for on-line demand forecasting and using web services API for developers, http://www.lokad.com/Developers.ashx . This company located in Paris, France provides a good working model on how to develop and use an SOA architecture for time series forecasting purposes.

References

Yu, J. (2002) Forecasting volatility in the New Zealand stock market. Applied Financial Economics 12, 193—202.

Yu, J. (2005) On leverage in a stochastic volatility model. Journal ofEconometrics 127, 165—178.