Skip to content

Commit 0bb1f0a

Browse files
author
Johan Dahlin
committed
added example 5
Former-commit-id: 85e621a [formerly 85e621a [formerly 5a8e5ae]] Former-commit-id: f7707cb608ff625bc2206530516cfab3b909b590 Former-commit-id: c6d1c7e
1 parent 3139d44 commit 0bb1f0a

2 files changed

Lines changed: 160 additions & 0 deletions

File tree

‎r/example5-sv.R‎

Lines changed: 159 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,159 @@
1+
##############################################################################
2+
#
3+
# Example of particle Metropolis-Hastings
4+
# in a stochastic volatility model
5+
#
6+
#
7+
# (c) 2015 Johan Dahlin
8+
# johan.dahlin (at) liu.se
9+
#
10+
##############################################################################
11+
12+
# Import helpers
13+
library("mvtnorm")
14+
15+
source("stateEstimationHelper.R")
16+
source("parameterEstimationHelper.R")
17+
18+
# Set the random seed to replicate results in tutorial
19+
set.seed(10)
20+
21+
##############################################################################
22+
# Define the model
23+
##############################################################################
24+
25+
# Here, we use the following model
26+
#
27+
# x[tt+1] = phi * x[tt] + sigma * v[tt]
28+
# y[tt] = beta * exp( xt[tt]/2 ) * e[tt]
29+
#
30+
# where v[tt] ~ N(0,1) and e[tt] ~ N(0,1)
31+
32+
# Set the number of time steps to simulate
33+
T = 500;
34+
35+
##############################################################################
36+
# Load data
37+
##############################################################################
38+
y = as.numeric( read.table("omxs30data.csv")[,1] )
39+
40+
##############################################################################
41+
# Parameter estimation using PMH
42+
##############################################################################
43+
44+
# The inital guess of the parameter
45+
initPar = c( 0, 0.9, 0.2);
46+
47+
# No. particles in the particle filter ( choose nPart ~ T )
48+
nPart = 500;
49+
50+
# The length of the burn-in and the no. iterations of PMH ( nBurnIn < nIter )
51+
nBurnIn = 2500;
52+
nIter = 7500;
53+
54+
# The standard deviation in the random walk proposal
55+
stepSize = matrix( c( 0.104402823, 0.008707826, -0.009090245,
56+
0.008707826, 0.071040325, -0.034589832,
57+
-0.009090245, -0.034589832, 0.086980437), ncol=3, nrow=3);
58+
stepSize = 0.8 * stepSize;
59+
60+
# Run the PMH algorithm
61+
res = pmh_sv_reparameterised(y,initPar,nPart,T,nIter,stepSize)
62+
63+
##############################################################################
64+
# Plot the results
65+
##############################################################################
66+
67+
# Extract the states after burn-in
68+
resTh = res$thhat[ nBurnIn:nIter, ];
69+
resXh = res$xhat[ nBurnIn:nIter, ];
70+
71+
# Estimate the KDE of the marginal posteriors
72+
kde1 = density( resTh[ , 1], kernel = "e", from = -1, to = 1 );
73+
kde2 = density( resTh[ , 2], kernel = "e", to = 1 );
74+
kde3 = density( resTh[ , 3], kernel = "e", from = 0 );
75+
76+
# Estimate the posterior mean and the corresponding standard deviation
77+
thhat = colMeans( resTh );
78+
thhatSD = apply( resTh, 2, sd );
79+
80+
# Estimate the log-volatility and the corresponding standad deviation
81+
xhat = colMeans( resXh );
82+
xhatSD = apply( resXh, 2, sd );
83+
84+
# Export plot to file
85+
cairo_pdf("example5-sv.pdf", height = 10, width = 8)
86+
87+
# Plot the parameter posterior estimate, solid black line indicate posterior mean
88+
# Plot the trace of the Markov chain after burn-in, solid black line indicate posterior mean
89+
#layout( matrix( c(1,1,1,2,2,2,3,4,5,6,7,8,9,10,11), 5, 3, byrow = TRUE ) );
90+
par( mar = c(4,5,0,0) );
91+
92+
# Grid for plotting the data and log-volatility
93+
gridy = seq( 1, length( y ), 1 );
94+
95+
plot( y, col = "#1B9E77", lwd = 1, type = "l", xlab = "time", ylab = "log-returns", bty = "n", ylim = c(-5,5) )
96+
polygon( c( gridy, rev( gridy ) ),c( y, rep( -5, length( y ) ) ), border = NA, col = rgb( t( col2rgb( "#1B9E77" ) ) / 256, alpha = 0.25) )
97+
98+
plot( xhat[-1], col = "#D95F02", lwd = 1.5, type = "l", xlab = "time", ylab = "log-volatility estimate", bty = "n", ylim = c(-2,2) )
99+
polygon( c( gridy, rev( gridy ) ), c( xhat[-1] - 1.96 * xhatSD[-1], rev( xhat[-1] + 1.96 * xhatSD[-1] ) ), border = NA,col = rgb( t( col2rgb( "#D95F02" ) ) /256, alpha = 0.25))
100+
101+
nPlot = 1500;
102+
grid = seq( nBurnIn, nBurnIn+nPlot-1, 1 );
103+
104+
# Mu
105+
hist( resTh[,1], breaks = floor( sqrt( nIter - nBurnIn ) ), col = rgb( t( col2rgb( "#7570B3" ) ) / 256, alpha = 0.25 ),border=NA,xlab=expression(mu),ylab="posterior estimate",main="",xlim=c(-1,1),freq=F)
106+
lines( seq(-1,1,0.01), dnorm(seq(-1,1,0.01),0,1),col="darkgrey")
107+
lines( kde1, lwd=2,col="#7570B3"); abline(v=thhat[1],lwd=1,lty="dotted");
108+
109+
plot(grid,resTh[1:nPlot,1],col='#7570B3', type="l",xlab="iteration",ylab=expression(mu),bty="n",ylim=c(-1,1),xlim=c(2500,4000))
110+
polygon(c(grid,rev(grid)),c(resTh[1:nPlot,1],rep(-1,length(resTh[1:nPlot,1]))),border=NA,col=rgb(t(col2rgb("#7570B3"))/256,alpha=0.25))
111+
abline(h=thhat[1],lwd=1,lty="dotted")
112+
113+
foo1=acf( resTh[,1], plot=F, lag.max=100)
114+
plot(foo1$lag,foo1$acf,col = '#7570B3', type="l",xlab="iteration",ylab=expression("ACF of "* mu),bty="n",lwd=2,ylim=c(-0.5,1))
115+
polygon(c(foo1$lag,rev(foo1$lag)),c(foo1$acf,rep(0,length(foo1$acf))),border=NA,col=rgb(t(col2rgb("#7570B3"))/256,alpha=0.25))
116+
abline(h=1.96/sqrt(nIter-nBurnIn),lty="dotted")
117+
abline(h=-1.96/sqrt(nIter-nBurnIn),lty="dotted")
118+
119+
# Phi
120+
hist(resTh[,2], breaks=floor(sqrt(nIter-nBurnIn)), col=rgb(t(col2rgb("#E7298A"))/256,alpha=0.25),border=NA,xlab=expression(phi),ylab="posterior estimate",main="",xlim=c(0.88,1.0),freq=F)
121+
lines(seq(0.88,1,0.01),dnorm(seq(0.88,1,0.01),0.95,0.05),col="darkgrey")
122+
lines(kde2,lwd=2,col="#E7298A"); abline(v=thhat[2],lwd=1,lty="dotted");
123+
plot(grid,resTh[1:nPlot,2],col='#E7298A', type="l",xlab="iteration",ylab=expression(phi),bty="n",ylim=c(0.9,1.1),xlim=c(2500,4000))
124+
125+
polygon(c(grid,rev(grid)),c(resTh[1:nPlot,2],rep(0.9,length(resTh[1:nPlot,2]))),border=NA,col=rgb(t(col2rgb("#E7298A"))/256,alpha=0.25))
126+
abline(h=thhat[2],lwd=1,lty="dotted")
127+
128+
foo2=acf( resTh[,2], plot=F, lag.max=100)
129+
plot(foo2$lag,foo2$acf,col = '#E7298A', type="l",xlab="iteration",ylab=expression("ACF of "* phi),bty="n",lwd=2,ylim=c(-0.5,1))
130+
polygon(c(foo2$lag,rev(foo2$lag)),c(foo2$acf,rep(0,length(foo2$acf))),border=NA,col=rgb(t(col2rgb("#E7298A"))/256,alpha=0.25))
131+
abline(h=1.96/sqrt(nIter-nBurnIn),lty="dotted")
132+
abline(h=-1.96/sqrt(nIter-nBurnIn),lty="dotted")
133+
134+
# Sigma[v]
135+
hist(resTh[,3], breaks=floor(sqrt(nIter-nBurnIn)), col=rgb(t(col2rgb("#66A61E"))/256,alpha=0.25),border=NA,xlab=expression(sigma[v]),ylab="posterior estimate",main="",xlim=c(0.0,0.4),freq=F)
136+
lines(seq(0,0.4,0.01),dgamma(seq(0,0.4,0.01),2,10),col="darkgrey")
137+
lines(kde3,lwd=2,col="#66A61E"); abline(v=thhat[3],lwd=1,lty="dotted");
138+
139+
plot(grid,resTh[1:nPlot,3],col='#66A61E', type="l",xlab="iteration",ylab=expression(sigma[v]),bty="n",ylim=c(0.0,0.4),xlim=c(2500,4000))
140+
polygon(c(grid,rev(grid)),c(resTh[1:nPlot,3],rep(0,length(resTh[1:nPlot,3]))),border=NA,col=rgb(t(col2rgb("#66A61E"))/256,alpha=0.25))
141+
abline(h=thhat[3],lwd=1,lty="dotted")
142+
143+
foo3=acf( resTh[,3], plot=F, lag.max=100)
144+
plot(foo3$lag,foo3$acf,col = '#66A61E', type="l",xlab="iteration",ylab=expression("ACF of "* sigma[v]),bty="n",lwd=2,ylim=c(-0.5,1))
145+
polygon(c(foo3$lag,rev(foo3$lag)),c(foo3$acf,rep(0,length(foo3$acf))),border=NA,col=rgb(t(col2rgb("#66A61E"))/256,alpha=0.25))
146+
abline(h=1.96/sqrt(nIter-nBurnIn),lty="dotted")
147+
abline(h=-1.96/sqrt(nIter-nBurnIn),lty="dotted")
148+
149+
#dev.off()
150+
151+
# Compute an estimate of the IACT using the first 100 ACF coefficients
152+
iact = 1 + 2 * c( sum(foo1$acf), sum(foo2$acf), sum(foo3$acf) )
153+
154+
# Estimate the covariance of the posterior to tune the proposal
155+
estCov = var( res$thhattansformed[ nBurnIn:nIter, ] )
156+
157+
##############################################################################
158+
# End of file
159+
##############################################################################

‎r/example5-sv.RData.REMOVED.git-id‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1 @@
1+
fa687828acff273072efc10811d6f4457957ee6d

0 commit comments

Comments
 (0)