gaussian_mixture_filter.m 10.8 KB
Newer Older
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
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
function [LIK,lik] = gaussian_mixture_filter(ReducedForm,Y,start,DynareOptions)
% Evaluates the likelihood of a non-linear model approximating the state
% variables distributions with gaussian mixtures. Gaussian Mixture allows reproducing
% a wide variety of generalized distributions (when multimodal for instance).
% Each gaussian distribution is obtained whether
%   - with a Smolyak quadrature à la Kronrod & Paterson (Heiss & Winschel 2010, Winschel & Kratzig 2010).
%   - with a radial-spherical cubature
%   - with scaled unscented sigma-points
% A Sparse grid Kalman Filter is implemented on each component of the mixture,
% which confers it a weight about current information.
% Information on the current observables is then embodied in the proposal
% distribution in which we draw particles, which allows
%   - reaching a greater precision relatively to a standard particle filter,
%   - reducing the number of particles needed,
%   - still being faster.
%
%
% 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'.
%
% OUTPUTS
%    LIK        [double]    scalar, likelihood
%    lik        [double]    vector, density of observations in each period.
%
% REFERENCES
%
% Van der Meerwe & Wan, Gaussian Mixture Sigma-Point Particle Filters for Sequential
% Probabilistic Inference in Dynamic State-Space Models.
% Heiss & Winschel, 2010, Journal of Applied Economics.
% Winschel & Kratzig, 2010, Econometrica.
%
% NOTES
%   The vector "lik" is used to evaluate the jacobian of the likelihood.
% Copyright (C) 2009-2013 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/>.


persistent init_flag mf0 mf1 Gprime Gsecond
persistent nodes weights weights_c I J G number_of_particles
persistent sample_size number_of_state_variables number_of_observed_variables number_of_structural_innovations

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

% Set persistent variables.
if isempty(init_flag)
  mf0 = ReducedForm.mf0;
  mf1 = ReducedForm.mf1;
  sample_size = size(Y,2);
  number_of_state_variables = length(mf0);
  number_of_observed_variables = length(mf1);
  number_of_structural_innovations = length(ReducedForm.Q);
  G = DynareOptions.particle.mixture_state_variables;           % number of GM components in state
  I = DynareOptions.particle.mixture_structural_shocks ;        % number of GM components in structural noise
  J = DynareOptions.particle.mixture_measurement_shocks ;       % number of GM components in observation noise
  Gprime = G*I ;
  Gsecond = G*I*J ;
  number_of_particles = DynareOptions.particle.number_of_particles;
  init_flag = 1;
end

SampleWeights = ones(Gsecond,1)/Gsecond ;

% compute gaussian quadrature nodes and weights on states and shocks
if isempty(nodes)
    if DynareOptions.particle.distribution_approximation.cubature
        [nodes,weights] = spherical_radial_sigma_points(number_of_state_variables);
        weights_c = weights;
    elseif DynareOptions.particle.distribution_approximation.unscented
        [nodes,weights,weights_c] = unscented_sigma_points(number_of_state_variables,DynareOptions);
    else
        if ~DynareOptions.particle.distribution_approximation.montecarlo
            error('Estimation: This approximation for the proposal is not implemented or unknown!')
        end
    end
end

if DynareOptions.particle.distribution_approximation.montecarlo
    set_dynare_seed('default');
    SampleWeights = 1/number_of_particles ;
end

% Get covariance matrices
Q = ReducedForm.Q;
H = ReducedForm.H;
if isempty(H)
    H = 0;
    H_lower_triangular_cholesky = 0;
else
    H_lower_triangular_cholesky = reduced_rank_cholesky(H)';
end
Q_lower_triangular_cholesky = reduced_rank_cholesky(Q)';

% Initialize all matrices
StateWeights = ones(1,G)/G ;
StateMu = ReducedForm.StateVectorMean*ones(1,G) ;
StateSqrtP = zeros(number_of_state_variables,number_of_state_variables,G) ;
for g=1:G
  StateSqrtP(:,:,g) = reduced_rank_cholesky(ReducedForm.StateVectorVariance)' ;
end

StructuralShocksWeights = ones(1,I)/I ;
StructuralShocksMu = zeros(number_of_structural_innovations,I) ;
StructuralShocksSqrtP = zeros(number_of_structural_innovations,number_of_structural_innovations,I) ;
for i=1:I
  StructuralShocksSqrtP(:,:,i) = Q_lower_triangular_cholesky ;
end

ObservationShocksWeights  = ones(1,J)/J ;
ObservationShocksMu = zeros(number_of_observed_variables,J) ;
ObservationShocksSqrtP = zeros(number_of_observed_variables,number_of_observed_variables,J) ;
for j=1:J
  ObservationShocksSqrtP(:,:,j) = H_lower_triangular_cholesky ;
end

StateWeightsPrior = zeros(1,Gprime) ;
StateMuPrior = zeros(number_of_state_variables,Gprime) ;
StateSqrtPPrior = zeros(number_of_state_variables,number_of_state_variables,Gprime) ;

StateWeightsPost = zeros(1,Gsecond) ;
StateMuPost = zeros(number_of_state_variables,Gsecond) ;
StateSqrtPPost = zeros(number_of_state_variables,number_of_state_variables,Gsecond) ;

%estimate = zeros(sample_size,number_of_state_variables,3) ;
const_lik = (2*pi)^(.5*number_of_observed_variables) ;

ks = 0 ;
lik  = NaN(sample_size,1);
LIK  = NaN;
for t=1:sample_size
    % Build the proposal joint quadratures of Gaussian on states, structural
    % shocks and observation shocks based on each combination of mixtures
    for i=1:I
        for j=1:J
            for g=1:G ;
                a = g + (j-1)*G ;
                b = a + (i-1)*Gprime ;
                [StateMuPrior(:,a),StateSqrtPPrior(:,:,a),StateWeightsPrior(1,a),...
                 StateMuPost(:,b),StateSqrtPPost(:,:,b),StateWeightsPost(1,b)] =...
                 gaussian_mixture_filter_bank(ReducedForm,Y(:,t),StateMu(:,g),StateSqrtP(:,:,g),StateWeights(1,g),...
                                                                 StructuralShocksMu(:,i),StructuralShocksSqrtP(:,:,i),StructuralShocksWeights(1,i),...
                                                                 ObservationShocksMu(:,j),ObservationShocksSqrtP(:,:,j),ObservationShocksWeights(1,j),...
                                                                 H,H_lower_triangular_cholesky,const_lik,DynareOptions) ;
            end
        end
    end

    % Normalize weights
    StateWeightsPrior = StateWeightsPrior/sum(StateWeightsPrior,2) ;
    StateWeightsPost = StateWeightsPost/sum(StateWeightsPost,2) ;

    if DynareOptions.particle.distribution_approximation.cubature || DynareOptions.particle.distribution_approximation.unscented
        for i=1:Gsecond
            StateParticles = bsxfun(@plus,StateMuPost(:,i),StateSqrtPPost(:,:,i)*nodes') ;
            IncrementalWeights = gaussian_mixture_densities(Y(:,t),StateMuPrior,StateSqrtPPrior,StateWeightsPrior,...
                                                                   StateMuPost,StateSqrtPPost,StateWeightsPost,...
                                                                   StateParticles,H,const_lik,weights,weights_c,ReducedForm,DynareOptions) ;
            SampleWeights(i) = sum(StateWeightsPost(i)*weights.*IncrementalWeights) ;
        end
        SumSampleWeights = sum(SampleWeights) ;
        lik(t) = log(SumSampleWeights) ;
        SampleWeights = SampleWeights./SumSampleWeights ;
        [ras,SortedRandomIndx] = sort(rand(1,Gsecond));
        SortedRandomIndx = SortedRandomIndx(1:G);
        indx = index_resample(0,SampleWeights,DynareOptions) ;
        indx = indx(SortedRandomIndx) ;
        StateMu = StateMuPost(:,indx);
        StateSqrtP = StateSqrtPPost(:,:,indx);
        StateWeights = ones(1,G)/G ;
    else
        % Sample particle in the proposal distribution, ie the posterior state GM
        StateParticles = importance_sampling(StateMuPost,StateSqrtPPost,StateWeightsPost',number_of_particles) ;
        % Compute prior, proposal and likelihood of particles
        IncrementalWeights = gaussian_mixture_densities(Y(:,t),StateMuPrior,StateSqrtPPrior,StateWeightsPrior,...
                                                               StateMuPost,StateSqrtPPost,StateWeightsPost,...
                                                               StateParticles,H,const_lik,1/number_of_particles,...
                                                               1/number_of_particles,ReducedForm,DynareOptions) ;
        % calculate importance weights of particles
        SampleWeights = SampleWeights.*IncrementalWeights ;
        SumSampleWeights = sum(SampleWeights,1) ;
        SampleWeights = SampleWeights./SumSampleWeights ;
        lik(t) = log(SumSampleWeights) ;
        % First possible state point estimates
        %estimate(t,:,1) = SampleWeights*StateParticles' ;
        % Resampling if needed of required
        Neff = 1/sum(bsxfun(@power,SampleWeights,2)) ;
        if (DynareOptions.particle.resampling.status.generic && Neff<.5*sample_size) || DynareOptions.particle.resampling.status.systematic
            ks = ks + 1 ;
            StateParticles = resample(StateParticles',SampleWeights,DynareOptions)' ;
            StateVectorMean = mean(StateParticles,2) ;
            StateVectorVarianceSquareRoot = reduced_rank_cholesky( (StateParticles*StateParticles')/number_of_particles - StateVectorMean*(StateVectorMean') )';
            SampleWeights = 1/number_of_particles ;
        elseif DynareOptions.particle.resampling.status.none
            StateVectorMean = StateParticles*sampleWeights ;
            temp = sqrt(SampleWeights').*StateParticles ;
            StateVectorVarianceSquareRoot = reduced_rank_cholesky( temp*temp' - StateVectorMean*(StateVectorMean') )';
        end
        % Use the information from particles to update the gaussian mixture on state variables
        [StateMu,StateSqrtP,StateWeights] = fit_gaussian_mixture(StateParticles,StateMu,StateSqrtP,StateWeights,0.001,10,1) ;
        %estimate(t,:,3) = StateWeights*StateMu' ;
    end
end

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