SV-AR(1) model

Brute-force MCMC vs Sequential MC

1 Simulating returns following a SV-AR(1) structure

library(stochvol)
set.seed(51413)
n     = 500
mu    = -7.3
phi   = 0.98
tau   = 0.18
alpha = mu*(1-phi)
true  = c(alpha,phi,tau)
h     = rep(mu,n)
for (t in 2:n)
  h[t] = alpha+phi*h[t-1]+rnorm(1,0,tau)
y = rnorm(n,0,exp(h/2))

par(mfrow=c(2,1))
ts.plot(exp(h/2),ylab="",main="Time-varying standard deviation")
ts.plot(y,ylab="",main="Simulated returns")

2 Particle filter (pure filter)

set.seed(1234)
m0 = mu
C0 = 1
M  = 10000
h  = rnorm(M,m0,sqrt(C0))
hs = matrix(0,M,n)
pftime = system.time({
for (t in 1:n){
  h1 = rnorm(M,alpha+phi*h,tau)
  w  = dnorm(y[t],0,exp(h1/2))
  h  = sample(h1,replace=TRUE,size=M,prob=w)
  hs[,t] = h
}
})
qsig.smc = t(apply(exp(hs/2),2,quantile,c(0.05,0.5,0.95)))

par(mfrow=c(1,1))
ts.plot(qsig.smc,col=c(3,2,3),ylim=range(qsig.smc),ylab="standard deviation")
lines(0.2*abs(y),type="h")

3 Markov chain Monte Carlo (MCMC)

ndraw = 10000
mcmctime = system.time({
  svfit = svsample(y,draws=ndraw,quiet=TRUE)
})
param = svfit$para[[1]][,1:3]
stdev = exp(svfit$latent[[1]]/2)
qmcmc = t(apply(stdev,2,quantile,c(0.025,0.5,0.975)))

par(mfrow=c(1,1))
ts.plot(qsig.smc,ylim=range(qsig.smc,qmcmc),ylab="Standard deviation")
lines(qmcmc[,1],col=2)
lines(qmcmc[,2],col=2)
lines(qmcmc[,3],col=2)

4 Brute-force MCMC

ts     = seq(5,n,by=5)
nt     = length(ts)
params = array(0,c(nt,ndraw,3))
sigs   = matrix(0,nt,ndraw)
bfmcmctime = system.time({
for (t in 1:nt){
  svfit       = svsample(y[1:ts[t]],draws=ndraw,quiet=TRUE)
  params[t,,] = svfit$para[[1]][,1:3]
  sigs[t,]   = exp(svfit$latent[[1]]/2)[,ts[t]]
}
})
alphas  = t(apply(params[,,1]*(1-params[,,2]),1,quantile,c(0.05,0.5,0.95)))
phis    = t(apply(params[,,2],1,quantile,c(0.05,0.5,0.95)))
taus    = t(apply(params[,,3],1,quantile,c(0.05,0.5,0.95)))
qsig.bf = t(apply(sigs,1,quantile,c(0.05,0.5,0.95)))

par(mfrow=c(2,2))
plot(ts,alphas[,2],ylim=range(alphas),main=expression(alpha),type="l",ylab="")
lines(ts,alphas[,1],lty=2)
lines(ts,alphas[,3],lty=2)
abline(h=true[1],col=2)

plot(ts,phis[,2],ylim=range(phis),main=expression(phi),type="l",ylab="")
lines(ts,phis[,1],lty=2)
lines(ts,phis[,3],lty=2)
abline(h=true[2],col=2)

plot(ts,taus[,2],ylim=range(taus),main=expression(tau),type="l",ylab="")
lines(ts,taus[,1],lty=2)
lines(ts,taus[,3],lty=2)
abline(h=true[3],col=2)

plot(ts,qsig.bf[,2],ylim=range(qsig.bf),main=expression(sigma[t]),type="l",ylab="")
lines(ts,qsig.bf[,1],lty=2)
lines(ts,qsig.bf[,3],lty=2)

5 Comparing SMC and brute-force-MCMC

par(mfrow=c(1,1))
ts.plot(qsig.smc,ylim=range(qsig.smc,qsig.bf),ylab="Standard deviation")
lines(ts,qsig.bf[,1],col=2,type="b",pch=16,cex=0.5)
lines(ts,qsig.bf[,1],col=2,type="b",pch=16,cex=0.5)
lines(ts,qsig.bf[,3],col=2,type="b",pch=16,cex=0.5)
legend("topright",legend=c(
  paste("Sequential Monte Carlo (time=",round(pftime[3],1),"secs)",sep=""),
  paste("Brute-force MCMC (time=",round((ts[2]-ts[1])*bfmcmctime[3],1),"secs)",sep="")),
  col=1:2,lty=1,bty="n",lwd=2)
title(paste("SMC is ",round((ts[2]-ts[1])*bfmcmctime[3]/pftime[3])," faster than Brute-force MCMC",sep=""))