Learn R Programming
SMC (version 1.1)

sequentialMonteCarlo: The sequential Monte Carlo (SMC) algorithm

Description

Function for the doing sequential Monte Carlo algorithm given the propagation rule over time (via propagateFunc). This is the most general interface for implementing a new SMC strategy, by providing a new propagation rule.

See the sections Details, Required Functions and Optional Functions for explanation on the arguments and the return values of the arguments that are themselves functions.

Usage

sequentialMonteCarlo(nStreams,                      nPeriods,                      dimPerPeriod,                  propagateFunc,              resampCriterionFunc = NULL,    resampFunc          = NULL,    summaryFunc         = NULL,    nMHSteps            = 0,       MHUpdateFunc        = NULL,    nStreamsPreResamp   = NULL,    returnStreams       = FALSE,   returnLogWeights    = FALSE,   verboseLevel        = 0,       ...)

Arguments

nStreams
integer $>$ 0.
nPeriods
integer $>$ 0.
dimPerPeriod
integer $>$ 0.
propagateFunc
function of six arguments (currentPeriod, nStreamsToGenerate, lag1Streams, lag1LogWeights, startingStreams, ...).
resampCriterionFunc
function of four arguments (currentPeriod, currentStreams, currentLogWeights, ...).
resampFunc
function of four arguments (currentPeriod, currentStreams, currentLogWeights, ...).
summaryFunc
function of four arguments (currentPeriod, currentStreams, currentLogWeights, ...).
nMHSteps
integer $>=$ 0.
MHUpdateFunc
function of six arguments (currentPeriod, nMHSteps, currentStreams, lag1Streams, lag1LogWeights, ...).
nStreamsPreResamp
integer $>$ 0.
returnStreams
logical.
returnLogWeights
logical.
verboseLevel
integer, a value $>=$ 2 produces a lot of output.
...
optional arguments to be passed to propagateFunc, resampCriterionFunc, resampFunc, summaryFunc and MHUpdateFunc.

Value

This function returns a list with the following components:
draws
a list with the following components: summary, propUniqueStreamIds, streams, logWeights, acceptanceRates. See the section Note for more details.
nStreams
the nStreams argument.
nPeriods
the nPeriods argument.
dimPerPeriod
the dimPerPeriod argument.
nStreamsPreResamp
the nStreamsPreResamp argument.
nMHSteps
the nMHSteps argument.
filterType
type of the filter: “sequentialMonteCarlo”.
time
the time taken by the run.

Required function: propagateFunc

Arguments:
The following argument(s) require some explanation:
nStreamsToGenerate
the number of streams to generate for propagating from currentPeriod - 1 to currentPeriod. This function is usally called by setting nStreamsToGenerate to nStreamsPreResamp.
lag1Streams
a matrix of dimension nStreams $x$ dimPerPeriod of streams for currentPeriod - 1.
lag1LogWeights
a vector of length nStreams of log weights corresponding to the streams in the argument matrix lag1Streams.
startingStreams
a matrix of dimension nStreams $x$ dimPerPeriod to be used for currentPeriod = 1. If this is NULL, then the function should provide a way to generate streams for currentPeriod = 1.
Return value:
a named list with the following components:
currentStreams
a matrix of dimension nStreamsToGenerate $x$ dimPerPeriod. The rows of this matrix contain the propagated (updated) streams for period currentPeriod, given the argument lag1Streams matrix and the argument lag1LogWeights vector for currentPeriod - 1.
currentLogWeights
the propagated (updated) log weights vector of length nStreamsToGenerate, associated with the streams in the rows of the returned currentStreams matrix.

Optional function: resampCriterionFunc

Arguments:
The following argument(s) require some explanation:
currentStreams
a matrix with dimPerPeriod columns, the rows containing the updated streams for currentPeriod.
currentLogWeights
a vector of log weights corresponding to the streams in the argument matrix currentStreams.
Return value:
TRUE or FALSE reflecting the decision of the resampling scheme implemented by this function.
Note:
The following points are in order:
--
resampling schemes manily depend on currentLogWeights, the other two arguments might come in handy for implementing period or stream specific resampling schemes.
--
if nStreamsPreResamp > nStreams, then this function should always return TRUE.

Optional function: resampFunc

Arguments:
see the sub-section Arguments: for section Optional function: resampCriterionFunc.
Return value:
a named list with the following components:
currentStreams
a matrix of dimension nStreams $x$ dimPerPeriod. The rows of this matrix contain the streams for period currentPeriod + 1 that were resampled from those of the argument currentStreams matrix, which may contain $>=$ nStreams rows.
currentLogWeights
The log weights vector of length nStreams, associated with the streams that were resampled in the returned currentStreams matrix. Note, after the resampling step, usually all the log weights are set to 0.
Note:
the components of the list returned by this function and the arguments to this function have two common names, namely, currentStreams and currentLogWeights. These entities have different meanings, as explained above. For example, the argument matrix currentStreams could possibly have $>=$ nStreams rows, whereas the returned currentStreams has exactly nStreams number of (resampled) streams in its rows.

Optional function: summaryFunc

Arguments:
The following argument(s) require some explanation:
currentStreams
a matrix of dimension nStreams $x$ dimPerPeriod of streams for currentPeriod.
currentLogWeights
a vector of log weights corresponding to the streams in the argument matrix currentStreams.
Return value:
a vector of length of dimSummPerPeriod of summaries for currentPeriod given the currentStreams and the currentLogWeights.

Optional function: MHUpdateFunc

Arguments:
The following argument(s) require some explanation:
nMHSteps
the number of Metropolis Hastings (MH) steps (iterations) to be performed.
currentStreams
a matrix of dimension nStreams $x$ dimPerPeriod of streams for currentPeriod.
lag1Streams
a matrix of dimension nStreams $x$ dimPerPeriod of streams for currentPeriod - 1.
lag1LogWeights
a vector of length nStreams of log weights corresponding to the streams in the argument matrix lag1Streams.
Return value:
a named list with the following components:
currentStreams
a matrix of dimension nStreams $x$ dimPerPeriod. The rows of this matrix contain the streams for period currentPeriod that are (possibly) MH-updated versions of the rows of the argument currentStreams matrix.
acceptanceRates
a vector of length nStreams, representing the acceptance rates of the nMHSteps-many MH steps for each of the streams in the rows of the argument currentStreams matrix.
Note:
a positive value of nMHSteps performs as many MH steps on the rows of the argument currentStreams matrix. This is done to reduce the possible degeneracy after the resampling.

Warning

Using very small values ($<=$ 1e3) for nStreams might not give reliable results.

Details

We introduce the following terms, which will be used in the sections Required Function and Optional Function below:

stream
the state vector also called the particle, the hidden state or the latent variable. Below we will use the terms stream and state vector interchangeably.

dimPerPeriod
the dimension of the space, the state vectors live in.

References

Jun S. Liu (2001). Monte Carlo strategies for scientific computing. Springer. Chapter 3.

Jun S. Liu and Rong Chen (1998). Sequential Monte Carlo methods for dynamical systems. Journal of the American Statistical Association 98(443): 1032-1044.

See Also

particleFilter, auxiliaryParticleFilter

Examples

Run this code
MSObj  <- MarkovSwitchingFuncGenerator(-12345)
smcObj <-
    with(MSObj,
     {
         sequentialMonteCarlo(nStreams         = 5000,
                              nPeriods         = nrow(yy),
                              dimPerPeriod     = ncol(yy),
                              propagateFunc    = propagateFunc,
                              returnStreams    = TRUE,
                              returnLogWeights = TRUE,
                              verboseLevel     = 1)
     })
print(smcObj)
print(names(smcObj))
with(c(smcObj, MSObj),
 {
     par(mfcol = c(2, 1))
     plot(as.ts(yy),
          main     = expression('The data and the underlying regimes'),
          cex.main = 0.8,
          xlab     = 'period',
          ylab     = 'data and the regime means',
          cex.lab  = 0.8)
     lines(as.ts(mu), col = 2, lty = 2)
     plot(as.ts(draws$summary[1, ]),
          main     = expression('The underlying regimes and their estimates'),
          cex.main = 0.8,
          xlab     = 'period',
          ylab     = 'regime means',
          cex.lab  = 0.8)
     lines(as.ts(mu), col = 2, lty = 2)        
 })

MSObj  <- MarkovSwitchingFuncGenerator(-54321)
smcObj <-
    with(MSObj,
     {
         sequentialMonteCarlo(nStreams         = 5000,
                              nPeriods         = nrow(yy),
                              dimPerPeriod     = ncol(yy),
                              propagateFunc    = propagateFunc,
                              returnStreams    = TRUE,
                              returnLogWeights = TRUE,
                              verboseLevel     = 1)
     })
print(smcObj)
print(names(smcObj))
with(c(smcObj, MSObj),
 {
     par(mfcol = c(2, 1))
     plot(as.ts(yy),
          main     = expression('The data and the underlying regimes'),
          cex.main = 0.8,
          xlab     = 'period',
          ylab     = 'data and the regime means',
          cex.lab  = 0.8)
     lines(as.ts(mu), col = 2, lty = 2)
     plot(as.ts(draws$summary[1, ]),
          main     = expression('The underlying regimes and their estimates'),
          cex.main = 0.8,
          xlab     = 'period',
          ylab     = 'regime means',
          cex.lab  = 0.8)
     lines(as.ts(mu), col = 2, lty = 2)        
 })

Run the code above in your browser using DataLab