| Title: | Minimal Working Examples for Particle Metropolis-Hastings |
|---|---|
| Description: | Routines for state estimate in a linear Gaussian state space model and a simple stochastic volatility model using particle filtering. Parameter inference is also carried out in these models using the particle Metropolis-Hastings algorithm that includes the particle filter to provided an unbiased estimator of the likelihood. This package is a collection of minimal working examples of these algorithms and is only meant for educational use and as a start for learning to them on your own. |
| Authors: | Johan Dahlin |
| Maintainer: | Johan Dahlin <[email protected]> |
| License: | GPL-2 |
| Version: | 1.5 |
| Built: | 2024-10-14 06:24:19 UTC |
| Source: | CRAN |
Minimal working example of state estimation in a linear Gaussian state space model using Kalman filtering and a fully-adapted particle filter. The code estimates the bias and mean squared error (compared with the Kalman estimate) while varying the number of particles in the particle filter.
example1_lgss()example1_lgss()
The Kalman filter is a standard implementation without an input. The particle filter is fully adapted (i.e. takes the current observation into account when proposing new particles and computing the weights).
Returns a plot with the generated observations y and the difference in the state estimates obtained by the Kalman filter (the optimal solution) and the particle filter (with 20 particles). Furthermore, the function returns plots of the estimated bias and mean squared error of the state estimate obtained using the particle filter (while varying the number of particles) and the Kalman estimates.
The function returns a list with the elements:
y: The observations generated from the model.
x: The states generated from the model.
kfEstimate: The estimate of the state from the Kalman filter.
pfEstimate: The estimate of the state from the particle filter with 20 particles.
See Section 3.2 in the reference for more details.
Johan Dahlin [email protected]
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.
example1_lgss()example1_lgss()
Minimal working example of parameter estimation in a linear Gaussian state space model using the particle Metropolis-Hastings algorithm with a fully-adapted particle filter providing an unbiased estimator of the likelihood. The code estimates the parameter posterior for one parameter using simulated data.
example2_lgss(noBurnInIterations = 1000, noIterations = 5000, noParticles = 100, initialPhi = 0.5)example2_lgss(noBurnInIterations = 1000, noIterations = 5000, noParticles = 100, initialPhi = 0.5)
noBurnInIterations |
The number of burn-in iterations in the PMH algorithm.
This parameter must be smaller than |
noIterations |
The number of iterations in the PMH algorithm. 100 iterations takes about ten seconds on a laptop to execute. 5000 iterations are used in the reference below. |
noParticles |
The number of particles to use when estimating the likelihood. |
initialPhi |
The initial guess of the parameter phi. |
The Particle Metropolis-Hastings (PMH) algorithm makes use of a Gaussian random walk as the proposal for the parameter. The PMH algorithm is run using different step lengths in the proposal. This is done to illustrate the difficulty when tuning the proposal and the impact of a too small/large step length.
Returns the estimate of the posterior mean.
See Section 4.2 in the reference for more details.
Johan Dahlin [email protected]
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.
example2_lgss(noBurnInIterations=200, noIterations=1000)example2_lgss(noBurnInIterations=200, noIterations=1000)
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.
example3_sv(noBurnInIterations = 2500, noIterations = 7500, noParticles = 500, initialTheta = c(0, 0.9, 0.2), stepSize = diag(c(0.1, 0.01, 0.05)^2), syntheticData = FALSE)example3_sv(noBurnInIterations = 2500, noIterations = 7500, noParticles = 500, initialTheta = c(0, 0.9, 0.2), stepSize = diag(c(0.1, 0.01, 0.05)^2), syntheticData = FALSE)
noBurnInIterations |
The number of burn-in iterations in the PMH
algorithm. Must be smaller than |
noIterations |
The number of iterations in the PMH algorithm. 100 iterations takes about a minute on a laptop to execute. |
noParticles |
The number of particles to use when estimating the likelihood. |
initialTheta |
The initial guess of the parameters theta. |
stepSize |
The step sizes of the random walk proposal. Given as a covariance matrix. |
syntheticData |
If TRUE, data is not downloaded from the Internet. This is only used when running tests of the package. |
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 somewhat well-tuned proposal as a pilot run to estimate the posterior covariance and therefore increase the mixing of the Markov chain.
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:
thhat: The estimate of the mean of the parameter posterior.
xhat: The estimate of the mean of the log-volatility posterior.
thhatSD: The estimate of the standard deviation of the parameter posterior.
xhatSD: The estimate of the standard deviation of the log-volatility posterior.
iact: The estimate of the integrated autocorrelation time for each parameter.
estCov: The estimate of the covariance of the parameter posterior.
theta: The trace of the chain exploring the parameter posterior.
See Section 5 in the reference for more details.
Johan Dahlin [email protected]
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.
## Not run: example3_sv(noBurnInIterations=200, noIterations=1000) ## End(Not run)## Not run: example3_sv(noBurnInIterations=200, noIterations=1000) ## End(Not run)
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.
example4_sv(noBurnInIterations = 2500, noIterations = 7500, noParticles = 500, initialTheta = c(0, 0.9, 0.2), syntheticData = FALSE)example4_sv(noBurnInIterations = 2500, noIterations = 7500, noParticles = 500, initialTheta = c(0, 0.9, 0.2), syntheticData = FALSE)
noBurnInIterations |
The number of burn-in iterations in the PMH
algorithm. Must be smaller than |
noIterations |
The number of iterations in the PMH algorithm. 100 iterations takes about a minute on a laptop to execute. |
noParticles |
The number of particles to use when estimating the likelihood. |
initialTheta |
The initial guess of the parameters theta. |
syntheticData |
If TRUE, data is not downloaded from the Internet. This is only used when running tests of the package. |
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 run of
example3_sv and therefore have better mixing
properties.
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:
thhat: The estimate of the mean of the parameter posterior.
xhat: The estimate of the mean of the log-volatility posterior.
thhatSD: The estimate of the standard deviation of the parameter posterior.
xhatSD: The estimate of the standard deviation of the log-volatility posterior.
iact: The estimate of the integrated autocorrelation time for each parameter.
estCov: The estimate of the covariance of the parameter posterior.
See Section 6.3.1 in the reference for more details.
Johan Dahlin [email protected]
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.
## Not run: example4_sv(noBurnInIterations=200, noIterations=1000) ## End(Not run)## Not run: example4_sv(noBurnInIterations=200, noIterations=1000) ## End(Not run)
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.
example5_sv(noBurnInIterations = 2500, noIterations = 7500, noParticles = 500, initialTheta = c(0, 0.9, 0.2), syntheticData = FALSE)example5_sv(noBurnInIterations = 2500, noIterations = 7500, noParticles = 500, initialTheta = c(0, 0.9, 0.2), syntheticData = FALSE)
noBurnInIterations |
The number of burn-in iterations in the PMH
algorithm. Must be smaller than |
noIterations |
The number of iterations in the PMH algorithm. 100 iterations takes about a minute on a laptop to execute. |
noParticles |
The number of particles to use when estimating the likelihood. |
initialTheta |
The initial guess of the parameters theta. |
syntheticData |
If TRUE, data is not downloaded from the Internet. This is only used when running tests of the package. |
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.
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:
thhat: The estimate of the mean of the parameter posterior.
xhat: The estimate of the mean of the log-volatility posterior.
thhatSD: The estimate of the standard deviation of the parameter posterior.
xhatSD: The estimate of the standard deviation of the log-volatility posterior.
iact: The estimate of the integrated autocorrelation time for each parameter.
estCov: The estimate of the covariance of the parameter posterior.
See Section 6.3.2 in the reference for more details.
Johan Dahlin [email protected]
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.
## Not run: example5_sv(noBurnInIterations=200, noIterations=1000) ## End(Not run)## Not run: example5_sv(noBurnInIterations=200, noIterations=1000) ## End(Not run)
Generates data from a specific linear Gaussian state space model of the form
and , where and denote independent standard
Gaussian random variables, i.e. .
generateData(theta, noObservations, initialState)generateData(theta, noObservations, initialState)
theta |
The parameters |
noObservations |
The number of time points to simulate. |
initialState |
The initial state. |
The function returns a list with the elements:
x: The latent state for .
y: The observation for .
Johan Dahlin [email protected]
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.
Estimates the filtered state and the log-likelihood for a linear Gaussian
state space model of the form
and , where and denote
independent standard Gaussian random variables, i.e..
kalmanFilter(y, theta, initialState, initialStateCovariance)kalmanFilter(y, theta, initialState, initialStateCovariance)
y |
Observations from the model for |
theta |
The parameters |
initialState |
The initial state. |
initialStateCovariance |
The initial covariance of the state. |
The function returns a list with the elements:
xHatFiltered: The estimate of the filtered state at time .
logLikelihood: The estimate of the log-likelihood.
See Section 3 in the reference for more details.
Johan Dahlin [email protected]
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.
# Generates 500 observations from a linear state space model with # (phi, sigma_e, sigma_v) = (0.5, 1.0, 0.1) and zero initial state. theta <- c(0.5, 1.0, 0.1) d <- generateData(theta, noObservations=500, initialState=0.0) # Estimate the filtered state using Kalman filter kfOutput <- kalmanFilter(d$y, theta, initialState=0.0, initialStateCovariance=0.01) # Plot the estimate and the true state par(mfrow=c(3, 1)) plot(d$x, type="l", xlab="time", ylab="true state", bty="n", col="#1B9E77") plot(kfOutput$xHatFiltered, type="l", xlab="time", ylab="Kalman filter estimate", bty="n", col="#D95F02") plot(d$x-kfOutput$xHatFiltered, type="l", xlab="time", ylab="difference", bty="n", col="#7570B3")# Generates 500 observations from a linear state space model with # (phi, sigma_e, sigma_v) = (0.5, 1.0, 0.1) and zero initial state. theta <- c(0.5, 1.0, 0.1) d <- generateData(theta, noObservations=500, initialState=0.0) # Estimate the filtered state using Kalman filter kfOutput <- kalmanFilter(d$y, theta, initialState=0.0, initialStateCovariance=0.01) # Plot the estimate and the true state par(mfrow=c(3, 1)) plot(d$x, type="l", xlab="time", ylab="true state", bty="n", col="#1B9E77") plot(kfOutput$xHatFiltered, type="l", xlab="time", ylab="Kalman filter estimate", bty="n", col="#D95F02") plot(d$x-kfOutput$xHatFiltered, type="l", xlab="time", ylab="difference", bty="n", col="#7570B3")
Creates diagnoistic plots from runs of the particle Metropolis-Hastings algorithm.
makePlotsParticleMetropolisHastingsSVModel(y, res, noBurnInIterations, noIterations, nPlot)makePlotsParticleMetropolisHastingsSVModel(y, res, noBurnInIterations, noIterations, nPlot)
y |
Observations from the model for |
res |
The output from a run of particleMetropolisHastings, particleMetropolisHastingsSVmodel or particleMetropolisHastingsSVmodelReparameterised. |
noBurnInIterations |
The number of burn-in iterations in the PMH algorithm. Must be smaller than noIterations. |
noIterations |
The number of iterations in the PMH algorithm. |
nPlot |
Number of steps in the Markov chain to plot. |
The function returns plots similar to the ones in the reference as well as the estimate of the integrated autocorrelation time for each parameter.
Johan Dahlin [email protected]
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.
Estimates the filtered state and the log-likelihood for a linear Gaussian
state space model of the form
and , where and denote
independent standard Gaussian random variables, i.e..
particleFilter(y, theta, noParticles, initialState)particleFilter(y, theta, noParticles, initialState)
y |
Observations from the model for |
theta |
The parameters |
noParticles |
The number of particles to use in the filter. |
initialState |
The initial state. |
The function returns a list with the elements:
xHatFiltered: The estimate of the filtered state at time .
logLikelihood: The estimate of the log-likelihood.
particles: The particle system at each time point.
weights: The particle weights at each time point.
See Section 3 in the reference for more details.
Johan Dahlin [email protected]
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.
# Generates 500 observations from a linear state space model with # (phi, sigma_e, sigma_v) = (0.5, 1.0, 0.1) and zero initial state. theta <- c(0.5, 1.0, 0.1) d <- generateData(theta, noObservations=500, initialState=0.0) # Estimate the filtered state using a Particle filter pfOutput <- particleFilter(d$y, theta, noParticles = 50, initialState=0.0) # Plot the estimate and the true state par(mfrow=c(3, 1)) plot(d$x[1:500], type="l", xlab="time", ylab="true state", bty="n", col="#1B9E77") plot(pfOutput$xHatFiltered, type="l", xlab="time", ylab="paticle filter estimate", bty="n", col="#D95F02") plot(d$x[1:500]-pfOutput$xHatFiltered, type="l", xlab="time", ylab="difference", bty="n", col="#7570B3")# Generates 500 observations from a linear state space model with # (phi, sigma_e, sigma_v) = (0.5, 1.0, 0.1) and zero initial state. theta <- c(0.5, 1.0, 0.1) d <- generateData(theta, noObservations=500, initialState=0.0) # Estimate the filtered state using a Particle filter pfOutput <- particleFilter(d$y, theta, noParticles = 50, initialState=0.0) # Plot the estimate and the true state par(mfrow=c(3, 1)) plot(d$x[1:500], type="l", xlab="time", ylab="true state", bty="n", col="#1B9E77") plot(pfOutput$xHatFiltered, type="l", xlab="time", ylab="paticle filter estimate", bty="n", col="#D95F02") plot(d$x[1:500]-pfOutput$xHatFiltered, type="l", xlab="time", ylab="difference", bty="n", col="#7570B3")
Estimates the filtered state and the log-likelihood for a stochastic
volatility model of the form and , where and
denote independent standard Gaussian random variables, i.e. .
particleFilterSVmodel(y, theta, noParticles)particleFilterSVmodel(y, theta, noParticles)
y |
Observations from the model for |
theta |
The parameters |
noParticles |
The number of particles to use in the filter. |
The function returns a list with the elements:
xHatFiltered: The estimate of the filtered state at time .
logLikelihood: The estimate of the log-likelihood.
See Section 5 in the reference for more details.
Johan Dahlin [email protected]
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.
## 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 filtered state using a particle filter theta <- c(-0.10, 0.97, 0.15) pfOutput <- particleFilterSVmodel(y, theta, noParticles=100) # Plot the estimate and the true state par(mfrow=c(2, 1)) plot(y, type="l", xlab="time", ylab="log-returns", bty="n", col="#1B9E77") plot(pfOutput$xHatFiltered, type="l", xlab="time", ylab="estimate of log-volatility", bty="n", col="#D95F02") ## End(Not run)## 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 filtered state using a particle filter theta <- c(-0.10, 0.97, 0.15) pfOutput <- particleFilterSVmodel(y, theta, noParticles=100) # Plot the estimate and the true state par(mfrow=c(2, 1)) plot(y, type="l", xlab="time", ylab="log-returns", bty="n", col="#1B9E77") plot(pfOutput$xHatFiltered, type="l", xlab="time", ylab="estimate of log-volatility", bty="n", col="#D95F02") ## End(Not run)
Estimates the parameter posterior for a linear Gaussian state
space model of the form and
, where and denote
independent standard Gaussian random variables, i.e..
particleMetropolisHastings(y, initialPhi, sigmav, sigmae, noParticles, initialState, noIterations, stepSize)particleMetropolisHastings(y, initialPhi, sigmav, sigmae, noParticles, initialState, noIterations, stepSize)
y |
Observations from the model for |
initialPhi |
The mean of the log-volatility process |
sigmav |
The standard deviation of the state process |
sigmae |
The standard deviation of the observation process
|
noParticles |
The number of particles to use in the filter. |
initialState |
The inital state. |
noIterations |
The number of iterations in the PMH algorithm. |
stepSize |
The standard deviation of the Gaussian random walk proposal
for |
The trace of the Markov chain exploring the marginal posterior for
.
See Section 4 in the reference for more details.
Johan Dahlin [email protected]
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.
# Generates 100 observations from a linear state space model with # (phi, sigma_e, sigma_v) = (0.5, 1.0, 0.1) and zero initial state. theta <- c(0.5, 1.0, 0.1) d <- generateData(theta, noObservations=100, initialState=0.0) # Estimate the marginal posterior for phi pmhOutput <- particleMetropolisHastings(d$y, initialPhi=0.1, sigmav=1.0, sigmae=0.1, noParticles=50, initialState=0.0, noIterations=1000, stepSize=0.10) # Plot the estimate nbins <- floor(sqrt(1000)) par(mfrow=c(1, 1)) hist(pmhOutput, breaks=nbins, main="", xlab=expression(phi), ylab="marginal posterior", freq=FALSE, col="#7570B3")# Generates 100 observations from a linear state space model with # (phi, sigma_e, sigma_v) = (0.5, 1.0, 0.1) and zero initial state. theta <- c(0.5, 1.0, 0.1) d <- generateData(theta, noObservations=100, initialState=0.0) # Estimate the marginal posterior for phi pmhOutput <- particleMetropolisHastings(d$y, initialPhi=0.1, sigmav=1.0, sigmae=0.1, noParticles=50, initialState=0.0, noIterations=1000, stepSize=0.10) # Plot the estimate nbins <- floor(sqrt(1000)) par(mfrow=c(1, 1)) hist(pmhOutput, breaks=nbins, main="", xlab=expression(phi), ylab="marginal posterior", freq=FALSE, col="#7570B3")
Estimates the parameter posterior for in
a stochastic volatility model of the form and , where and
denote independent standard Gaussian random variables, i.e.
.
particleMetropolisHastingsSVmodel(y, initialTheta, noParticles, noIterations, stepSize)particleMetropolisHastingsSVmodel(y, initialTheta, noParticles, noIterations, stepSize)
y |
Observations from the model for |
initialTheta |
An inital value for the parameters
|
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 |
The trace of the Markov chain exploring the posterior of .
See Section 5 in the reference for more details.
Johan Dahlin [email protected]
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.
## 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 <- particleMetropolisHastingsSVmodel(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)## 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 <- particleMetropolisHastingsSVmodel(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)
Estimates the parameter posterior for in
a stochastic volatility model of the form and , where and
denote independent standard Gaussian random variables, i.e.
. In this version of the PMH, we reparameterise the model and
run the Markov chain on the parameters , where and .
particleMetropolisHastingsSVmodelReparameterised(y, initialTheta, noParticles, noIterations, stepSize)particleMetropolisHastingsSVmodelReparameterised(y, initialTheta, noParticles, noIterations, stepSize)
y |
Observations from the model for |
initialTheta |
An inital value for the parameters
|
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 |
The trace of the Markov chain exploring the posterior of .
See Section 5 in the reference for more details.
Johan Dahlin [email protected]
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.
## 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)## 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)