1 Introduction

In this exercise, MCMC for the SV-AR(1) model will be performed via R package stochvol, while the SMC for the SV-AR(1) is performed via my own bootstrap filter code (conditional on static parameters)

1.1 Loadings packages and reading the data

library(stochvol)
#library(tidyquant)
#getSymbols("PBR",from='2000-08-10',to='2021-03-22',warnings=FALSE,auto.assign=TRUE)
#n0    = nrow(PBR)
#price = as.numeric(PBR[,6])

price = scan("pbr-081000-031921.txt")
n0    = length(price)
y0    = log(price[2:n0]/price[1:(n0-1)])
n0    = n0 - 1
y0    = y0-mean(y0)
n1    = 2000
y1    = y0[(n0-n1+1):n0]
n1    = length(y1)

par(mfrow=c(1,2))
ts.plot(price,main="PBR closing prices\n 08/10/00-03/22/21",ylab="")
ts.plot(y0,main="PBR log returns",ylab="")
abline(v=(n0-n1),col=2)

2 MCMC

Let us use stochvol to learn about log-volatilities and parameters for the whole time span. The package fits the following version of the SV-AR(1) model \[\begin{eqnarray*} y_t &=& \exp\{h_t/2\} \varepsilon_t \\ h_t &=& \mu + \phi (h_{t-1}-\mu) + \sigma \eta_t, \end{eqnarray*}\] where \(\varepsilon_t\) and \(\eta_t\) are iid standard normal with \(E(\varepsilon_t,\eta_{t'})=0\) for all \(t,t'\). The default prior hyperparameters are such that \[\begin{eqnarray*} \mu &\sim& N(0,100^2)\\ (\phi+1)/2 &\sim& Beta(5,1.5)\\ \sigma^2 &\sim& IG(1/2,1/2)\\ h_0 &\sim& N(\mu,\sigma^2/(1-\phi^2)). \end{eqnarray*}\] Finally, we also use their default for burn-in sample (\(1000\)) and number of draws for posterior inference (\(10,000\)).

ndraw = 10000
svfit = svsample(y0,draws=ndraw)
## Done!
## Summarizing posterior draws...
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(2,2))
hist(param[,1],main="mu",prob=TRUE,xlab="")
hist(param[,2],main="phi",prob=TRUE,xlab="")
hist(param[,3],main="sigma",prob=TRUE,xlab="")
ts.plot(qmcmc[,2],ylab="Standard deviation")

3 Bootstrap filter

We will now implement a sequential Monte Carlo algorithm assuming that \(\mu\), \(\phi\) and \(\sigma\) are all known. Notice that in my bootstrap filter below I consider the following version for the evolution of log-volatilities in the SV-AR(1) model \[ h_t = \alpha + \phi h_{t-1} + \sigma \eta_t, \] so is easy to see that \(\alpha=\mu(1-\phi)\). Since we are fixing the static parameters at some pre-specified values, we are actually obtainng sequential Monte Carlo approximations to \[ p(h_t|y_1,\ldots,y_t,\mu,\phi,\sigma), \qquad t=1,\ldots,n, \] while stochvol is performing smoothed posterior inference via MCMC, i.e. \[ p(h_1,\ldots,h_n,\mu,\phi,\sigma|y_1,\ldots,y_n). \]

mu    = mean(param[,1])
phi   = mean(param[,2])
sigma = mean(param[,3])
alpha = mu*(1-phi)

# Particle filter
set.seed(1234)
h0 = log(stdev[,n0]^2)
m0 = mean(h0)
C0 = var(h0)
M  = ndraw
h  = rnorm(M,m0,sqrt(C0))
hs = matrix(0,M,n1)
system.time({
for (t in 1:n1){
  # propagation
  h1 = rnorm(M,alpha+phi*h,sigma)
  # resampling
  w = dnorm(y1[t],0,exp(h1/2))
  h = sample(h1,replace=TRUE,size=M,prob=w)
  hs[,t] = h
}
})
##    user  system elapsed 
##   7.008   1.114   8.408
qstdev = t(apply(exp(hs/2),2,quantile,c(0.05,0.5,0.95)))

3.1 Output

par(mfrow=c(1,2))
ts.plot(qmcmc[(n0-n1+1):n0,],col=c(3,2,3),ylim=range(qstdev,qmcmc[(n0-n1+1):n0,]))
title("MCMC\nstochvol")
ts.plot(qstdev,col=c(3,2,3),ylim=range(qstdev,qmcmc[(n1+1):n0,]))
title("SMC\nBootstrap filter with fixed parameters)")

par(mfrow=c(1,2))
sd1 = stdev[,n0]
sd2 = exp(hs[,n1]/2)
hist(sd1,prob=TRUE,main="",xlab="Last observation",xlim=range(sd1,sd2))
title("MCMC\nstochvol")
hist(sd2,prob=TRUE,main="",xlab="Last observation",xlim=range(sd1,sd2))
title("SMC\nBootstrap filter with fixed parameters")

4 Comparing sequential MCMC and bootstrap filter

ts = seq(100,n1,by=100)
nt = length(ts)
stdev = matrix(0,nt,ndraw)
system.time({
for (t in 1:nt){
  svfit  = svsample(y1[1:ts[t]],draws=ndraw,quiet=TRUE)
  stdev[t,] = exp(svfit$latent[[1]]/2)[,ts[t]]
}
})
##    user  system elapsed 
## 191.211  34.541 231.998
qmcmc1 = t(apply(stdev,1,quantile,c(0.025,0.5,0.975)))
par(mfrow=c(1,1))
ts.plot(qstdev,col=c(3,2,3),ylim=range(qstdev,qmcmc[(n1+1):n0,]))
title("Bootstrap filter (with fixed parameters) vs Sequential MCMC")
points(ts,qmcmc1[,2],col=6,pch=16)
points(ts,qmcmc1[,1],col=4,pch=16)
points(ts,qmcmc1[,3],col=4,pch=16)

par(mfrow=c(4,5),mai=c(0.3,0.3,0.3,0.1))
for (t in 1:nt){
  den1 = density(stdev[t,])
  den2 = density(exp(hs[,ts[t]]/2))
  plot(den1$x,den1$y,xlab="",main="",ylab="",ylim=range(den1$y,den2$y),xlim=range(den1$x,den2$x),type="l")
  lines(den2$x,den2$y,col=2)
  mtext(side=3,line=0,paste("t=",ts[t],sep=""),cex=0.75)
  if (t==1){
    legend("topright",legend=c("BF","sMCMC"),col=2:1,lty=1,bty="n")
  }
}

5 Bootstrap filter for the whole data

mu    = mean(param[,1])
phi   = mean(param[,2])
sigma = mean(param[,3])
alpha = mu*(1-phi)

# Filtro de particulas
set.seed(1234)
m0 = -6.6
C0 = 0.5
M  = ndraw
h  = rnorm(M,m0,sqrt(C0))
hs = matrix(0,M,n0)
system.time({
for (t in 1:n0){
  # propagacao
  h1 = rnorm(M,alpha+phi*h,sigma)
  # resampling
  w = dnorm(y0[t],0,exp(h1/2))
  h = sample(h1,replace=TRUE,size=M,prob=w)
  hs[,t] = h
}
})
##    user  system elapsed 
##  17.959   3.377  21.928
qstdev0 = t(apply(exp(hs/2),2,quantile,c(0.05,0.5,0.95)))

par(mfrow=c(1,2),mai=c(0.8,0.8,0.8,0.5))
ts.plot(qmcmc,col=c(3,2,3),ylim=range(qstdev,qmcmc))
title("MCMC\nstochvol")
ts.plot(qstdev0,col=c(3,2,3),ylim=range(qstdev,qmcmc))
title("SMC\nBootstrap filter with fixed parameters")

par(mfrow=c(2,2),mai=c(0.4,0.4,0.3,0.3))
for (i in 1:4){
  ind = ((i-1)*1295+1):(i*1295)
  plot(ind,qmcmc[ind,2],ylim=range(qstdev0[ind,2],qmcmc[ind,2]),type="l",xlab="",ylab="")
  lines(ind,qstdev0[ind,2],col=2)
  if (i==1){
    legend("topright",legend=c("MCMC","SMC"),col=2:1,lty=1,bty="n")
  }
}