Nothing
##############################################################################
# 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.
#' @param syntheticData If TRUE, data is not downloaded from the Internet. This is only used when running tests of the package.
#' @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."
#' Journal of Statistical Software, Code Snippets,
#' 88(2): 1--41, 2019.
#' @author
#' Johan Dahlin \email{uni@@johandahlin.com}
#' @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), syntheticData=FALSE) {
# Set the random seed to replicate results in tutorial
set.seed(10)
##############################################################################
# Load data
##############################################################################
if (syntheticData) {
y <- rnorm(10)
} else {
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(1500, 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)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.