conditional_particle_filter.m 4.93 KB
Newer Older
1
function [LIK,lik] = conditional_particle_filter(ReducedForm,Y,start,ParticleOptions,ThreadsOptions)
Stéphane Adjemian's avatar
Stéphane Adjemian committed
2
%
3
% Evaluates the likelihood of a non-linear model with a particle filter
Stéphane Adjemian's avatar
Stéphane Adjemian committed
4
5
6
7
8
9
10
% - the proposal is built using the Kalman updating step for each particle.
% - we need draws in the errors distributions
% Whether we use Monte-Carlo draws from a multivariate gaussian distribution
% as in Amisano & Tristani (JEDC 2010).
% Whether we use multidimensional Gaussian sparse grids approximations:
% - a univariate Kronrod-Paterson Gaussian quadrature combined by the Smolyak
% operator (ref: Winschel & Kratzig, 2010).
11
% - a spherical-radial cubature (ref: Arasaratnam & Haykin, 2009a,2009b).
Stéphane Adjemian's avatar
Stéphane Adjemian committed
12
% - a scaled unscented transform cubature (ref: Julier & Uhlmann 1997, van der
13
% Merwe & Wan 2003).
Stéphane Adjemian's avatar
Stéphane Adjemian committed
14
15
16
17
18
19
%
% Pros:
% - Allows using current observable information in the proposal
% - The use of sparse grids Gaussian approximation is much faster than the Monte-Carlo approach
% Cons:
% - The use of the Kalman updating step may biais the proposal distribution since
20
% it has been derived in a linear context and is implemented in a nonlinear
Stéphane Adjemian's avatar
Stéphane Adjemian committed
21
% context. That is why particle resampling is performed.
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
%
% INPUTS
%    reduced_form_model     [structure] Matlab's structure describing the reduced form model.
%                                       reduced_form_model.measurement.H   [double]   (pp x pp) variance matrix of measurement errors.
%                                       reduced_form_model.state.Q         [double]   (qq x qq) variance matrix of state errors.
%                                       reduced_form_model.state.dr        [structure] output of resol.m.
%    Y                      [double]    pp*smpl matrix of (detrended) data, where pp is the maximum number of observed variables.
%    start                  [integer]   scalar, likelihood evaluation starts at 'start'.
%    smolyak_accuracy       [integer]   scalar.
%
% OUTPUTS
%    LIK        [double]    scalar, likelihood
%    lik        [double]    vector, density of observations in each period.
%
% REFERENCES
%
% NOTES
%   The vector "lik" is used to evaluate the jacobian of the likelihood.
% Copyright (C) 2009-2010 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.

% AUTHOR(S) frederic DOT karame AT univ DASH lemans DOT fr
%           stephane DOT adjemian AT univ DASH lemans DOT fr

60
persistent init_flag mf1
Stéphane Adjemian's avatar
Stéphane Adjemian committed
61
62
persistent number_of_particles
persistent sample_size number_of_observed_variables
63
64
65
66
67
68
69
70
71
72
73
74

% Set default
if isempty(start)
    start = 1;
end

% Set persistent variables.
if isempty(init_flag)
    mf1 = ReducedForm.mf1;
    sample_size = size(Y,2);
    number_of_observed_variables = length(mf1);
    init_flag = 1;
75
    number_of_particles = ParticleOptions.number_of_particles ;
76
77
78
79
80
81
82
83
84
end

% Get covariance matrices
Q = ReducedForm.Q;
H = ReducedForm.H;
if isempty(H)
    H = 0;
    H_lower_triangular_cholesky = 0;
else
Stéphane Adjemian's avatar
Stéphane Adjemian committed
85
    H_lower_triangular_cholesky = chol(H)';
86
87
88
89
end

% Get initial condition for the state vector.
StateVectorMean = ReducedForm.StateVectorMean;
90
StateVectorVarianceSquareRoot = chol(ReducedForm.StateVectorVariance)';
91
state_variance_rank = size(StateVectorVarianceSquareRoot,2);
Stéphane Adjemian's avatar
Stéphane Adjemian committed
92
Q_lower_triangular_cholesky = chol(Q)';
93
94
95
96
97

% Set seed for randn().
set_dynare_seed('default');

% Initialization of the likelihood.
98
normconst2 = sqrt( ((2*pi)^number_of_observed_variables)*prod(det(H)) ) ;
99
100
101
102
103
104
lik  = NaN(sample_size,1);
LIK  = NaN;
ks = 0 ;
StateParticles = bsxfun(@plus,StateVectorVarianceSquareRoot*randn(state_variance_rank,number_of_particles),StateVectorMean);
SampleWeights = ones(1,number_of_particles)/number_of_particles ;
for t=1:sample_size
Stéphane Adjemian's avatar
Stéphane Adjemian committed
105
106
107
    for i=1:number_of_particles
        [StateParticles(:,i),SampleWeights(i)] = ...
            conditional_filter_proposal(ReducedForm,Y(:,t),StateParticles(:,i),SampleWeights(i),Q_lower_triangular_cholesky,H_lower_triangular_cholesky,H,ParticleOptions,ThreadsOptions,normconst2) ;
108
109
    end
    SumSampleWeights = sum(SampleWeights) ;
Stéphane Adjemian's avatar
Stéphane Adjemian committed
110
111
    lik(t) = log(SumSampleWeights) ;
    SampleWeights = SampleWeights./SumSampleWeights ;
112
    if (ParticleOptions.resampling.status.generic && neff(SampleWeights)<ParticleOptions.resampling.threshold*sample_size) || ParticleOptions.resampling.status.systematic
113
        ks = ks + 1 ;
114
        StateParticles = resample(StateParticles',SampleWeights',ParticleOptions)';
115
116
117
118
        SampleWeights = ones(1,number_of_particles)/number_of_particles ;
    end
end

119
LIK = -sum(lik(start:end));