R/example5-sv.R

Defines functions example5_sv

Documented in example5_sv

##############################################################################
# Parameter estimation using particle Metropolis-Hastings in a reparameterised version of a
# stochastic volatility model with a proposal adapted from a pilot run.
#
# Johan Dahlin <uni (at) johandahlin.com.nospam>
# Documentation at https://github.com/compops/pmh-tutorial
# Published under GNU General Public License
##############################################################################

#' Parameter estimation in a simple stochastic volatility model
#' 
#' @description
#' Minimal working example of parameter estimation in a stochastic volatility 
#' model using the particle Metropolis-Hastings algorithm with a bootstrap 
#' particle filter providing an unbiased estimator of the likelihood. The 
#' code estimates the parameter posterior for three parameters using 
#' real-world data. 
#' @details 
#' The Particle Metropolis-Hastings (PMH) algorithm makes use of a Gaussian 
#' random walk as the proposal for the parameters. The data are scaled 
#' log-returns from the OMXS30 index during the period from January 2, 2012 
#' to January 2, 2014.
#' 
#' This version of the code makes use of a proposal that is tuned using a 
#' pilot run. Furthermore the model is reparameterised to enjoy better mixing 
#' properties by making the parameters unrestricted to a certain part of the 
#' real-line.
#' @param noBurnInIterations The number of burn-in iterations in the PMH 
#' algorithm. Must be smaller than \code{noIterations}.
#' @param noIterations The number of iterations in the PMH algorithm. 100 
#' iterations takes about a minute on a laptop to execute.
#' @param noParticles The number of particles to use when estimating the likelihood.
#' @param initialTheta The initial guess of the parameters theta.
#' @return 
#' The function returns the estimated marginal parameter posteriors for each 
#' parameter, the trace of the Markov chain and the resulting autocorrelation 
#' function. The data is also presented with an estimate of the 
#' log-volatility.
#' 
#' The function returns a list with the elements:
#' \itemize{
#' \item{thhat: The estimate of the mean of the parameter posterior.}
#' \item{xhat: The estimate of the mean of the log-volatility posterior.}
#' \item{thhatSD: The estimate of the standard deviation of the parameter 
#' posterior.}
#' \item{xhatSD: The estimate of the standard deviation of the log-volatility 
#' posterior.}
#' \item{iact: The estimate of the integrated autocorrelation time for each 
#' parameter.}
#' \item{estCov: The estimate of the covariance of the parameter posterior.}
#' }
#' @references 
#' Dahlin, J. & Schon, T. B. "Getting started with particle 
#' Metropolis-Hastings for inference in nonlinear dynamical models." 
#' pre-print, arXiv:1511.01707, 2017.
#' @author 
#' Johan Dahlin <uni (at) johandahlin.com.nospam>
#' @note 
#' See Section 6.3.2 in the reference for more details.
#' @example ./examples/example5
#' @keywords 
#' misc
#' @export
#' @importFrom grDevices col2rgb
#' @importFrom grDevices rgb
#' @importFrom graphics abline
#' @importFrom graphics hist
#' @importFrom graphics layout
#' @importFrom graphics lines
#' @importFrom graphics par
#' @importFrom graphics plot
#' @importFrom graphics points
#' @importFrom graphics polygon
#' @importFrom stats acf
#' @importFrom stats density
#' @importFrom stats sd
#' @importFrom stats var

example5_sv <- function(noBurnInIterations=2500, noIterations=7500, noParticles=500,
                         initialTheta=c(0, 0.9, 0.2)) {

  # Set the random seed to replicate results in tutorial
  set.seed(10)

  ##############################################################################
  # Load data
  ##############################################################################
  d <-
    Quandl::Quandl(
      "NASDAQOMX/OMXS30",
      start_date = "2012-01-02",
      end_date = "2014-01-02",
      type = "zoo"
    )
  y <- as.numeric(100 * diff(log(d$"Index Value")))


  ##############################################################################
  # PMH
  ##############################################################################
  stepSize <- matrix(
    c(
      0.041871682,-0.001200581,-0.002706803,-0.001200581,
      0.054894707,-0.056321320,-0.002706803,-0.056321320,
      0.087342276
    ),
    ncol = 3,
    nrow = 3
  )
  stepSize <- 2.562^2 / 3 * stepSize

  res <- particleMetropolisHastingsSVmodelReparameterised(y, initialTheta, noParticles, noIterations, stepSize)

  ##############################################################################
  # Plot the results
  ##############################################################################
  noIterationsToPlot <- min(c(100, noIterations-noBurnInIterations))
  iact <- makePlotsParticleMetropolisHastingsSVModel(y, res, noBurnInIterations, 
                                                     noIterations, noIterationsToPlot)

  ##############################################################################
  # Compute and save the results
  ##############################################################################  

  # Extract the states after burn-in
  resTh <- res$theta[noBurnInIterations:noIterations, ]
  resXh <- res$xHatFiltered[noBurnInIterations:noIterations, ]
  
  # Estimate the posterior mean and the corresponding standard deviation
  thhat   <- colMeans(resTh)
  thhatSD <- apply(resTh, 2, sd)
  
  # Estimate the log-volatility and the corresponding standad deviation
  xhat    <- colMeans(resXh)
  xhatSD  <- apply(resXh, 2, sd)

  # Estimate the covariance of the posterior to tune the proposal
  estCov <- var(resTh)

  list(thhat=thhat, xhat=xhat, thhatSD=thhatSD, xhatSD=xhatSD, iact=iact, estCov=estCov)
}

Try the pmhtutorial package in your browser

Any scripts or data that you put into this service are public.

pmhtutorial documentation built on Oct. 7, 2017, 5:02 p.m.