pmhtutorial (version 1.5)

particleMetropolisHastingsSVmodelReparameterised: Particle Metropolis-Hastings algorithm for a stochastic volatility model model

Description

Estimates the parameter posterior for \(\theta=\{\mu,\phi,\sigma_v\}\) in a stochastic volatility model of the form \(x_t = \mu + \phi ( x_{t-1} - \mu ) + \sigma_v v_t\) and \(y_t = \exp(x_t/2) e_t\), where \(v_t\) and \(e_t\) denote independent standard Gaussian random variables, i.e. \(N(0,1)\). In this version of the PMH, we reparameterise the model and run the Markov chain on the parameters \(\vartheta=\{\mu,\psi, \varsigma\}\), where \(\phi=\tanh(\psi)\) and \(sigma_v=\exp(\varsigma)\).

Usage

particleMetropolisHastingsSVmodelReparameterised(y, initialTheta,
  noParticles, noIterations, stepSize)

Arguments

y

Observations from the model for \(t=1,...,T\).

initialTheta

An inital value for the parameters \(\theta=\{\mu,\phi,\sigma_v\}\). The mean of the log-volatility process is denoted \(\mu\). The persistence of the log-volatility process is denoted \(\phi\). The standard deviation of the log-volatility process is denoted \(\sigma_v\).

noParticles

The number of particles to use in the filter.

noIterations

The number of iterations in the PMH algorithm.

stepSize

The standard deviation of the Gaussian random walk proposal for \(\theta\).

Value

The trace of the Markov chain exploring the posterior of \(\theta\).

References

Dahlin, J. & Schon, T. B. "Getting Started with Particle Metropolis-Hastings for Inference in Nonlinear Dynamical Models." Journal of Statistical Software, Code Snippets, 88(2): 1--41, 2019.

Examples

# NOT RUN {
  # Get the data from Quandl
  library("Quandl")
  d <- Quandl("NASDAQOMX/OMXS30", start_date="2012-01-02",
              end_date="2014-01-02", type="zoo")
  y <- as.numeric(100 * diff(log(d$"Index Value")))

  # Estimate the marginal posterior for phi
  pmhOutput <- particleMetropolisHastingsSVmodelReparameterised(
    y, initialTheta = c(0, 0.9, 0.2), noParticles=500,
    noIterations=1000, stepSize=diag(c(0.05, 0.0002, 0.002)))

  # Plot the estimate
  nbins <- floor(sqrt(1000))
  par(mfrow=c(3, 1))
  hist(pmhOutput$theta[,1], breaks=nbins, main="", xlab=expression(mu),
    ylab="marginal posterior", freq=FALSE, col="#7570B3")
  hist(pmhOutput$theta[,2], breaks=nbins, main="", xlab=expression(phi),
    ylab="marginal posterior", freq=FALSE, col="#E7298A")
  hist(pmhOutput$theta[,3], breaks=nbins, main="",
    xlab=expression(sigma[v]), ylab="marginal posterior",
    freq=FALSE, col="#66A61E")
# }