particleMetropolisHastingsSVmodelReparameterised {pmhtutorial}

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.

Note

See Section 5 in the reference for more details.

Author(s)

Johan Dahlin [e-mail]

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")

## End(Not run)