# sequentialMonteCarlo: The sequential Monte Carlo (SMC) algorithm In SMC: 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

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14``` ```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`.

## 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.

## 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.

## Note

The effect of leaving the default value `NULL` for some of the arguments above are as follows:

`resampCriterionFunc`

the builtin resampling criterion, namely, resample when square of the coefficient of variation of the weights >= 1, is used.

`resampFunc`

the builtin resampling function, which resamples streams with probability proportional to their weights, is used.

`summaryFunc`

the builtin summary function, which returns the weighted average of each of the `dimPerPeriod` dimensions, is used.

`MHUpdateFunc`

unlike, `particleFilter`, there is no builtin Metropolis Hastings updating function, which generates proposals for `currentPeriod` streams using those of ```currentPeriod - 1```. The user needs to implement this function if `nMHSteps` > 0.

`nStreamsPreResamp`

it is set to `nStreams`.

Also, the following point is worth noting:

`resampCriterionFunc`, `resampFunc`, `summaryFunc`

are only necessary when user wants to try out new resampling schemes or enhanced summary generation procedures, as part of their research. The default builtins take care of the typical problems.

This function returns a list with component called `draw`. The detailed description of this component, as promised in section Value, is as follows. It is a list itself with the following components:

`summary`

a matrix of dimension `nPeriods` x `dimSummPerPeriod`.

`propUniqueStreamIds`

a vector of length `nPeriods`. The values are either proportions of unique stream ids accpeted (at each period) if resampling was done or `NA`.

`streams`

an array of dimension `nStreams` x `dimPerPeriod` x `nPeriods`. This is returned if `returnStreams = TRUE`.

`logWeights`

a matrix of dimension `nStreams` x `nPeriods`. This is returned if `returnLogWeights = TRUE`.

`acceptanceRates`

a matrix of dimension `nStreams` x `nPeriods`. This is returned if `nMHSteps > 0`.

## Author(s)

Gopi Goswami goswami@stat.harvard.edu

## 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.

`particleFilter`, `auxiliaryParticleFilter`

## Examples

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65``` ```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) }) ```

### Example output

```##
## Sequential Monte Carlo Package (SMC)
##
## Functionality: sequential Monte Carlo (SMC) or sequential importance
## sampling (SIS) or hidden Markov models (HMM), particle filter (PF)
## and auxiliary particle filter (APF)
##
## Use: "help(package = SMC)" at the R prompt for more info
##
## Copyright (C) 2006-2019 Gopi Goswami
##
##    Created by: Gopi Goswami <goswami@stat.harvard.edu>
## Maintained by: Gopi Goswami <grgoswami@gmail.com>
##

BEGIN: SMC
..........
[Time to finish (est):        1 secs, this period:       10]....................
[Time to finish (est):        1 secs, this period:       30]....................
[Time to finish (est):        1 secs, this period:       50]....................
[Time to finish (est):        1 secs, this period:       70]....................
[Time to finish (est):        1 secs, this period:       90]..........
[Total time:          0 secs (usr),     0 secs (sys), this period:        100]
E N D: SMC

The settings for this sequentialMonteCarlo run:

nStreams:          5000
nStreamsPreResamp: 5000
nPeriods:           100
dimPerPeriod:         1

The summary of the resampling proportions:
[Note: resampling was done for 11 periods]

Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's
0.0920  0.2510  0.3152  0.3104  0.4170  0.4556      89

[1] "draws"             "nStreams"          "nPeriods"
[4] "dimPerPeriod"      "nStreamsPreResamp" "nMHSteps"
[7] "summaryFunc"       "time"              "propagateFunc"
[10] "filterType"

BEGIN: SMC
..........
[Time to finish (est):        1 secs, this period:       10]....................
[Time to finish (est):        1 secs, this period:       30]....................
[Time to finish (est):        1 secs, this period:       50]....................
[Time to finish (est):        1 secs, this period:       70]....................
[Time to finish (est):        1 secs, this period:       90]..........
[Total time:          0 secs (usr),     0 secs (sys), this period:        100]
E N D: SMC

The settings for this sequentialMonteCarlo run:

nStreams:          5000
nStreamsPreResamp: 5000
nPeriods:           100
dimPerPeriod:         1

The summary of the resampling proportions:
[Note: resampling was done for 16 periods]

Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's
0.0636  0.1475  0.2528  0.2329  0.3118  0.4294      84

[1] "draws"             "nStreams"          "nPeriods"
[4] "dimPerPeriod"      "nStreamsPreResamp" "nMHSteps"
[7] "summaryFunc"       "time"              "propagateFunc"
[10] "filterType"
```

SMC documentation built on May 1, 2019, 9:12 p.m.