The auxiliary particle filtering algorithm
Description
Function for doing auxiliary particle filtering given the state
equation (via generateNextStreamsFunc
), the stream
representative generation rule (via generateStreamRepsFunc
),
and the observation equation density (via
logObsDensFunc
).
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 15 16  auxiliaryParticleFilter(nStreams,
nPeriods,
dimPerPeriod,
generateStreamRepsFunc,
generateNextStreamsFunc,
logObsDensFunc,
resampCriterionFunc = NULL,
resampFunc = NULL,
summaryFunc = NULL,
nMHSteps = 0,
MHUpdateFunc = NULL,
nStreamsPreResamp = NULL,
returnStreams = FALSE,
returnLogWeights = FALSE,
verboseLevel = 0,
...)

Arguments
nStreams 

nPeriods 

dimPerPeriod 

generateStreamRepsFunc 

generateNextStreamsFunc 

logObsDensFunc 

resampCriterionFunc 

resampFunc 

summaryFunc 

nMHSteps 

MHUpdateFunc 

nStreamsPreResamp 

returnStreams 

returnLogWeights 

verboseLevel 

... 
optional arguments to be passed to

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: 
nStreams 
the 
nPeriods 
the 
dimPerPeriod 
the 
nStreamsPreResamp 
the 
nMHSteps 
the 
filterType 
type of the filter: “auxiliaryParticleFilter”. 
time 
the time taken by the run. 
Required function: generateStreamRepsFunc
 Arguments:
The following argument(s) require some explanation:
lag1Streams
a matrix of dimension
nStreams
xdimPerPeriod
of streams forcurrentPeriod  1
.lag1LogWeights
a vector of length
nStreams
of log weights corresponding to the streams in the argument matrixlag1Streams
.streamIndices
a vector of length
nStreams
for which the stream representatives (μ_t^k of Pitt and Shephard, 1999) forcurrentPeriod
are to be generated. See the subsection Note: below.
 Return value:
a matrix of dimension
nStreamIndices
xdimPerPeriod
. The rows of this matrix contain the stream representative for periodcurrentPeriod
, given the state vectors to be found in thestreamIndices
rows of the argumentlag1Streams
matrix. HerenStreamIndices
is the length of the argumentstreamIndices
.
 Note:
The following points are in order:
 –
this function should distinguish the cases
currentPeriod == 1
andcurrentPeriod > 1
inside of it. –
for details on the stream representatives (i.e., μ_t^k), see of Pitt and Shephard, 1999. The quantity μ_t^k could be the mean, the mode, a draw or some other likely value associated with the state density for period
currentPeriod
(i.e., f(α_t \mid α_{t  1})). –
this function is called by setting
streamIndices
to1:nStreams
, i.e., stream representatives for all the streams in the argumentlag1Streams
matrix is generated.
Optional function: generateNextStreamsFunc
 Arguments:
The following argument(s) require some explanation:
lag1Streams
a matrix of dimension
nStreams
xdimPerPeriod
of streams forcurrentPeriod  1
.lag1LogWeights
a vector of length
nStreams
of log weights corresponding to the streams in the argument matrixlag1Streams
.streamIndices
a vector of length >=
nStreams
which are to be updated fromcurrentPeriod  1
tocurrentPeriod
.streamReps
a matrix of dimension
nStreams
xdimPerPeriod
of the stream representatives forcurrentPeriod
.startingStreams
a matrix of dimension
nStreams
xdimPerPeriod
to be used forcurrentPeriod = 1
. If this isNULL
, then the function should provide a way to generate streams forcurrentPeriod = 1
.
 Return value:
a matrix of dimension
nStreamIndices
xdimPerPeriod
. The rows of this matrix contain the state vectors for periodcurrentPeriod
, given the state vectors to be found in thestreamIndices
rows of the argumentlag1Streams
matrix. HerenStreamIndices
is the length of the argumentstreamIndices
.
 Note:
The following points are in order:
 –
this function should distinguish the cases
currentPeriod == 1
andcurrentPeriod > 1
inside of it. –
this function is called by setting
streamIndices
such thatnStreamIndices
takes either of the two valuesnStreams
ornStreamsPreResamp
in different ocassions.
Optional function: logObsDensFunc
 Arguments:
The following argument(s) require some explanation:
currentStreams
a matrix with
dimPerPeriod
columns, the rows containing the streams forcurrentPeriod
.
 Return value:
a vector of length
nCurrentStreams
, wherenCurrentStreams
refers to the number of rows of thecurrentStreams
matrix argument. This vector contains the observation equation density values forcurrentPeriod
in the log scale, evaluated at the rows ofcurrentStreams
.
 Note:
nCurrentStreams
might be >=nStreams
.
Optional function: resampCriterionFunc
 Arguments:
The following argument(s) require some explanation:
currentStreams
a matrix with
dimPerPeriod
columns, the rows containing the updated streams forcurrentPeriod
.currentLogWeights
a vector of log weights corresponding to the streams in the argument matrix
currentStreams
.
 Return value:
TRUE
orFALSE
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 returnTRUE
.
Optional function: resampFunc
 Arguments:
see the subsection Arguments: for section Optional function: resampCriterionFunc.
 Return value:
a named list with the following components:
currentStreams
a matrix of dimension
nStreams
xdimPerPeriod
. The rows of this matrix contain the streams for periodcurrentPeriod + 1
that were resampled from those of the argumentcurrentStreams
matrix, which may contain >=nStreams
rows.currentLogWeights
The log weights vector of length
nStreams
, associated with the streams that were resampled in the returnedcurrentStreams
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
andcurrentLogWeights
. These entities have different meanings, as explained above. For example, the argument matrixcurrentStreams
could possibly have >=nStreams
rows, whereas the returnedcurrentStreams
has exactlynStreams
number of (resampled) streams in its rows.
Optional function: summaryFunc
 Arguments:
The following argument(s) require some explanation:
currentStreams
a matrix of dimension
nStreams
xdimPerPeriod
of streams forcurrentPeriod
.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 forcurrentPeriod
given thecurrentStreams
and thecurrentLogWeights
.
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
xdimPerPeriod
of streams forcurrentPeriod
.lag1Streams
a matrix of dimension
nStreams
xdimPerPeriod
of streams forcurrentPeriod  1
.lag1LogWeights
a vector of length
nStreams
of log weights corresponding to the streams in the argument matrixlag1Streams
.
 Return value:
a named list with the following components:
currentStreams
a matrix of dimension
nStreams
xdimPerPeriod
. The rows of this matrix contain the streams for periodcurrentPeriod
that are (possibly) MHupdated versions of the rows of the argumentcurrentStreams
matrix.acceptanceRates
a vector of length
nStreams
, representing the acceptance rates of thenMHSteps
MH steps for each of the streams in the rows of the argumentcurrentStreams
matrix.
 Note:
a positive value of
nMHSteps
performs as many MH steps on the rows of the argumentcurrentStreams
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 forcurrentPeriod
streams using those ofcurrentPeriod  1
. The user needs to implement this function ifnMHSteps
> 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
xdimSummPerPeriod
.propUniqueStreamIds
a vector of length
nPeriods
. The values are either proportions of unique streams accpeted (at each period) if resampling was done orNA
.streams
an array of dimension
nStreams
xdimPerPeriod
xnPeriods
. This is returned ifreturnStreams = TRUE
.logWeights
a matrix of dimension
nStreams
xnPeriods
. This is returned ifreturnLogWeights = TRUE
.acceptanceRates
a matrix of dimension
nStreams
xnPeriods
. This is returned ifnMHSteps > 0
.
Author(s)
Gopi Goswami goswami@stat.harvard.edu
References
Michael K. Pitt and Meil Shephard (1999). Filtering via Simulation: Auxiloary Particle Filters. Journal of the American Statistical Association 94(446): 590599.
See Also
particleFilter
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 66 67 68 69  MSObj < MarkovSwitchingFuncGenerator(2468)
smcObj <
with(MSObj,
{
auxiliaryParticleFilter(nStreams = 5000,
nPeriods = nrow(yy),
dimPerPeriod = ncol(yy),
generateStreamRepsFunc = generateStreamRepsFunc,
generateNextStreamsFunc = generateNextStreamsFunc,
logObsDensFunc = logObsDensFunc,
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(8642)
smcObj <
with(MSObj,
{
auxiliaryParticleFilter(nStreams = 5000,
nPeriods = nrow(yy),
dimPerPeriod = ncol(yy),
generateStreamRepsFunc = generateStreamRepsFunc,
generateNextStreamsFunc = generateNextStreamsFunc,
logObsDensFunc = logObsDensFunc,
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)
})

Want to suggest features or report bugs for rdrr.io? Use the GitHub issue tracker. Vote for new features on Trello.