In this example, we consider the Markov switching stochastic volatility model.

Reference: C.M. Carvalho and H.F. Lopes. Simulation-based sequential analysis of Markov switching stochastic volatility models. Computational Statistics and Data analysis (2007) 4526-4542.

Statistical model

Let \(y_t\) be the response variable and \(x_t\) the unobserved log-volatility of \(y_t\). The stochastic volatility model is defined as follows for \(t\leq t_{max}\)

\[x_t|(x_{t-1},\alpha,\phi,\sigma,c_t) \sim \mathcal N (\alpha_{c_t} + \phi x_{t-1} , \sigma^2)\]

\[ y_t|x_t \sim \mathcal N (0, \exp(x_t)) \]

The regime variables \(c_t\) follow a two-state Markov process with transition probabilities

\[p_{ij}=Pr(c_t=j|c_{t-1}=i)\]

\(\mathcal N(m,\sigma^2)\) denotes the normal distribution of mean \(m\) and variance \(\sigma^2\).

Statistical model in BUGS language

Content of the file 'switch_stoch_volatility.bug':

model_file = 'switch_stoch_volatility.bug' # BUGS model filename
cat(readLines(model_file), sep = "\n")
var y[t_max],x[t_max],mu[t_max],mu_true[t_max],c[t_max],c_true[t_max]
data
{
  c_true[1] ~ dcat(pi[c0,])
  mu_true[1] <- alpha[1] * (c_true[1]==1) + alpha[2]*(c_true[1]==2) + phi*x0
  x_true[1] ~ dnorm(mu_true[1], 1/sigma^2)
  y[1] ~ dnorm(0, exp(-x_true[1]))
  for (t in 2:t_max)
  {
    c_true[t] ~ dcat(ifelse(c_true[t-1]==1,pi[1,],pi[2,]))
    mu_true[t] <- alpha[1]*(c_true[t]==1) + alpha[2]*(c_true[t]==2) + phi*x_true[t-1];
    x_true[t] ~ dnorm(mu_true[t], 1/sigma^2)
    y[t] ~ dnorm(0, exp(-x_true[t]))
  }
}

model
{
  c[1] ~ dcat(pi[c0,])
  mu[1] <- alpha[1] * (c[1]==1) + alpha[2]*(c[1]==2)+ phi*x0
  x[1] ~ dnorm(mu[1], 1/sigma^2)
  y[1] ~ dnorm(0, exp(-x[1]))
  for (t in 2:t_max)
  {
    c[t] ~ dcat(ifelse(c[t-1]==1, pi[1,], pi[2,]))
    mu[t] <- alpha[1] * (c[t]==1) + alpha[2]*(c[t]==2) + phi*x[t-1]
    x[t] ~ dnorm(mu[t], 1/sigma^2)
    y[t] ~ dnorm(0, exp(-x[t]))
  }
}
  1. Load rbiips package
require(rbiips)

General settings

par(bty='l')
light_blue = rgb(.7, .7, 1)
light_red = rgb(1, .7, .7)
hot_colors = colorRampPalette(c('black', 'red', 'yellow', 'white'))

Set the random numbers generator seed for reproducibility

set.seed(0)

Load model and data

Model parameters

sigma = .4; alpha = c(-2.5, -1); phi = .5; c0 = 1; x0 = 0; t_max = 100
pi = matrix(c(.9, .1, .1, .9), nrow=2, byrow=TRUE)
data = list(t_max=t_max, sigma=sigma,
            alpha=alpha, phi=phi, pi=pi, c0=c0, x0=x0)

Parse and compile BUGS model, and sample data

model = biips_model(model_file, data, sample_data=TRUE)
* Parsing model in: switch_stoch_volatility.bug
* Compiling data graph
  Declaring variables
  Resolving undeclared variables
  Allocating nodes
  Graph size: 1215
  Sampling data
  Reading data back into data table
* Compiling model graph
  Declaring variables
  Resolving undeclared variables
  Allocating nodes
  Graph size: 1218
data = model$data()

Plot the data

plot(1:t_max, data$y, type='l', col='blue', lwd=2,
     main='Observed data',
     xlab='Time', ylab='Log-return', xaxt = 'n')
Log-returns

Log-returns

Biips Sequential Monte Carlo

Run SMC

n_part = 5000 # Number of particles
variables = c('x') # Variables to be monitored
out_smc = biips_smc_samples(model, variables, n_part)
* Assigning node samplers
* Running SMC forward sampler with 5000 particles
  |--------------------------------------------------| 100%
  |**************************************************| 100 iterations in 7.13 s

Diagnosis of the algorithm

diag_smc = biips_diagnosis(out_smc)
* Diagnosis of variable: x[1:100] 
  Filtering: GOOD
  Smoothing: GOOD

Plot Smoothing ESS

plot(out_smc$x$s$ess, type='l', log='y', col='blue', lwd=2,
     xlab='Time', ylab='SESS',
     ylim=c(10, n_part))
lines(1:t_max, rep(30,t_max), lwd=2, lty=2)
SMC: SESS

SMC: SESS

Plot weighted particles

matplot(1:t_max, out_smc$x$s$values, type='n',
        xlab='Time', ylab='Log-volatility')
for (t in 1:t_max) {
  val = unique(out_smc$x$s$values[t,])
  weight = sapply(val, FUN=function(x) {
    ind = out_smc$x$s$values[t,] == x
    return(sum(out_smc$x$s$weights[t,ind]))
  })
  points(rep(t, length(val)), val, col='red', pch=20,
         cex=sapply(weight, FUN=function(x) {min(2, .02*n_part*x)}))
}
lines(1:t_max, data$x_true, col='green', lwd=2)
legend('topright', leg=c('Particles (smoothing)', 'True value'),
       col=c('red','green'), lwd=c(NA,2), pch=c(20,NA),
       bty='n')
SMC: Particles (smoothing)

SMC: Particles (smoothing)

Summary statistics

summ_smc = biips_summary(out_smc, probs=c(.025, .975))

Plot Filtering estimates

x_f_mean = summ_smc$x$f$mean
x_f_quant = summ_smc$x$f$quant

xx = c(1:t_max, t_max:1)
yy = c(x_f_quant[[1]], rev(x_f_quant[[2]]))
plot(xx, yy, type='n', xlab='Time', ylab='Log-volatility')

polygon(xx, yy, col=light_blue, border=NA)
lines(1:t_max, x_f_mean, col='blue', lwd=3)
lines(1:t_max, data$x_true, col='green', lwd=2)
legend('topright', leg=c('95% credible interval', 'Filtering mean estimate', 'True value'),
       col=c(light_blue,'blue','green'), lwd=c(NA,3,2), pch=c(15,NA,NA), pt.cex=c(2,1,1),
       bty='n')
SMC: Filtering estimates

SMC: Filtering estimates

Plot Smoothing estimates

x_s_mean = summ_smc$x$s$mean
x_s_quant = summ_smc$x$s$quant

xx = c(1:t_max, t_max:1)
yy = c(x_s_quant[[1]], rev(x_s_quant[[2]]))
plot(xx, yy, type='n', xlab='Time', ylab='Log-volatility')

polygon(xx, yy, col=light_red, border=NA)
lines(1:t_max, x_s_mean, col='red', lwd=3)
lines(1:t_max, data$x_true, col='green', lwd=2)
legend('topright', leg=c('95% credible interval', 'Smoothing mean estimate', 'True value'),
       col=c(light_red,'red','green'), lwd=c(NA,3,2), pch=c(15,NA,NA), pt.cex=c(2,1,1),
       bty='n')
SMC: Smoothing estimates

SMC: Smoothing estimates

Marginal filtering and smoothing densities

kde_smc = biips_density(out_smc)
time_index = c(5, 10, 15)
par(mfrow=c(2,2))
for (k in 1:length(time_index)) {
  tk = time_index[k]
  plot(kde_smc$x[[tk]], col=c('blue', 'red'), lwd=2,
       xlab=bquote(x[.(tk)]), ylab='Posterior density',
       main=paste('t=', tk, sep=''))
  points(data$x_true[tk], 0, col='green', pch=8, lwd=2)
}
plot(0, type='n', bty='n', xaxt='n', yaxt='n', xlab="", ylab="")
legend('center', leg=c('Filtering density', 'Smoothing density', 'True value'),
       col=c('blue', 'red', 'green'), pch=c(NA,NA,8), lty=c(1,1,NA), lwd=2,
       bty='n')
par(mfrow=c(1,1))
SMC: Marginal posteriors

SMC: Marginal posteriors

Biips Particle Independent Metropolis-Hastings

Parameters of the PIMH

n_burn = 2000
n_iter = 10000
thin = 1
n_part = 50

Run PIMH

obj_pimh = biips_pimh_init(model, variables)
* Initializing PIMH
biips_pimh_update(obj_pimh, n_burn, n_part) # Burn-in iterations
* Updating PIMH with 50 particles
  |--------------------------------------------------| 100%
  |**************************************************| 2000 iterations in 110.52 s
out_pimh = biips_pimh_samples(obj_pimh, n_iter, n_part, thin=thin) # Return samples
* Generating PIMH samples with 50 particles
  |--------------------------------------------------| 100%
  |**************************************************| 10000 iterations in 551.96 s

Some summary statistics

summ_pimh = biips_summary(out_pimh, probs=c(.025, .975))

Posterior mean and quantiles

x_pimh_mean = summ_pimh$x$mean
x_pimh_quant = summ_pimh$x$quant

xx = c(1:t_max, t_max:1)
yy = c(x_pimh_quant[[1]], rev(x_pimh_quant[[2]]))
plot(xx, yy, type='n', xlab='Time', ylab='Log-volatility')

polygon(xx, yy, col=light_blue, border=NA)
lines(1:t_max, x_pimh_mean, col='blue', lwd=3)
lines(1:t_max, data$x_true, col='green', lwd=2)
legend('topright', leg=c('95% credible interval', 'PIMH mean estimate', 'True value'),
       col=c(light_blue,'blue','green'), lwd=c(NA,3,2), pch=c(15,NA,NA), pt.cex=c(2,1,1),
       bty='n')
PIMH: Posterior mean and quantiles

PIMH: Posterior mean and quantiles

Trace of MCMC samples

time_index = c(5, 10, 15)
par(mfrow=c(2,2))
for (k in 1:length(time_index)) {
  tk = time_index[k]
  plot(out_pimh$x[tk,], col='blue', type='l',
       xlab='Iteration', ylab=bquote(x[.(tk)]),
       main=paste('t=', tk, sep=''))
  points(0, data$x_true[tk], col='green', pch=8, lwd=2)
}
plot(0, type='n', bty='n', xaxt='n', yaxt='n', xlab="", ylab="")
legend('center', leg=c('PIMH samples', 'True value'),
       col=c('blue', 'green'), pch=c(NA,8), lwd=c(1,2), lty=c(1,NA),
       bty='n')
PIMH: Trace samples

PIMH: Trace samples

Histograms of posteriors

par(mfrow=c(2,2))
for (k in 1:length(time_index)) {
  tk = time_index[k]
  hist(out_pimh$x[tk,], breaks=20, col='blue', border='white',
       xlab=bquote(x[.(tk)]), ylab='Number of samples',
       main=paste('t=', tk, sep=''))
  points(data$x_true[tk], 0, col='green', pch=8, lwd=2)
}
plot(0, type='n', bty='n', xaxt='n', yaxt='n', xlab="", ylab="")
legend('center', leg=c('Posterior density', 'True value'),
       col=c('blue', 'green'), pch=c(22,8), lwd=c(NA,2), lty=NA, pt.cex=c(2,1), pt.bg=c(4,NA),
       bty='n')
PIMH: Histograms marginal posteriors

PIMH: Histograms marginal posteriors

Kernel density estimates of posteriors

par(mfrow=c(2,2))
for (k in 1:length(time_index)) {
  tk = time_index[k]
  kde_pimh = density(out_pimh$x[tk,])
  plot(kde_pimh, col='blue', lwd=2,
       xlab=bquote(x[.(tk)]), ylab='Posterior density',
       main=paste('t=', tk, sep=''))
  points(data$x_true[tk], 0, col='green', pch=8, lwd=2)
}
plot(0, type='n', bty='n', xaxt='n', yaxt='n', xlab="", ylab="")
legend('center', leg=c('Posterior density', 'True value'),
       col=c('blue', 'green'), pch=c(NA,8), lwd=2, lty=c(1,NA),
       bty='n')
par(mfrow=c(1,1))
PIMH: KDE estimates marginal posteriors

PIMH: KDE estimates marginal posteriors

Biips Sensitivity analysis

We want to study sensitivity to the values of the parameter \(\alpha\)

Parameters of the algorithm

n_part = 50 # Number of particles
grid <- seq(-5,2,.2) # Grid of values for one component
A = rep(grid, times=length(grid)) # Values of the first component
B = rep(grid, each=length(grid)) # Values of the second component
param_values = list('alpha' = rbind(A, B))

Run sensitivity analysis with SMC

out_sens = biips_smc_sensitivity(model, param_values, n_part)
* Analyzing sensitivity with 50 particles
  |--------------------------------------------------| 100%
  |**************************************************| 1296 iterations in 71.27 s

Plot log-marginal likelihood and penalized log-marginal likelihood

# Avoid scaling problems by thresholding the values
thres = -40
zz = matrix(pmax(thres, out_sens$log_marg_like), nrow=length(grid))

require(lattice)
Le chargement a nécessité le package : lattice
levelplot(zz, row.values=grid, column.values=grid, bg='black',
          at=seq(thres, max(zz), len=100), col.regions=hot_colors(100),
          xlim=range(grid), ylim=range(grid),
          xlab=expression(alpha[1]), ylab=expression(alpha[2]),
          interpolate=TRUE, useRaster=TRUE, bty='l')
Sensitivity: Log-likelihood

Sensitivity: Log-likelihood