online_auxiliary_filter.m 16.5 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
function [xparam,std_param,lb_95,ub_95,median_param] = online_auxiliary_filter(xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults)

% Carvalho & Lopes particle filter = auxiliary particle filter including Liu & West filter on parameters.
%
% INPUTS
%    ReducedForm     [structure] Matlab's structure describing the reduced form model.
%                                       ReducedForm.measurement.H   [double]   (pp x pp) variance matrix of measurement errors.
%                                       ReducedForm.state.Q         [double]   (qq x qq) variance matrix of state errors.
%                                       ReducedForm.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'.
%    mf                     [integer]   pp*1 vector of indices.
%    number_of_particles    [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) 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 Y init_flag mf0 mf1 bounds number_of_particles number_of_parameters liu_west_delta liu_west_chol_sigma_bar
persistent start_param sample_size number_of_observed_variables number_of_structural_innovations

% Set seed for randn().
set_dynare_seed('default') ;
pruning = DynareOptions.particle.pruning;
46
47
second_resample = 1 ;
variance_update = 1 ;
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63

% initialization of state particles
[ys,trend_coeff,exit_flag,info,Model,DynareOptions,BayesInfo,DynareResults,ReducedForm] = ... 
            solve_model_for_online_filter(1,xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults) ;

% Set persistent variables.
if isempty(init_flag)
    mf0 = ReducedForm.mf0;
    mf1 = ReducedForm.mf1;
    number_of_particles = DynareOptions.particle.number_of_particles;
    number_of_parameters = size(xparam1,1) ;
    Y = DynareDataset.data ;
    sample_size = size(Y,2);
    number_of_observed_variables = length(mf1);
    number_of_structural_innovations = length(ReducedForm.Q);
    liu_west_delta = DynareOptions.particle.liu_west_delta ;
Frédéric Karamé's avatar
Frédéric Karamé committed
64
65
66
67
    %liu_west_chol_sigma_bar = DynareOptions.particle.liu_west_chol_sigma_bar*eye(number_of_parameters) ;
    start_param = xparam1 ; 
    %liu_west_chol_sigma_bar = sqrt(bsxfun(@times,eye(number_of_parameters),BayesInfo.p2)) ; 
    %start_param = BayesInfo.p1 ;
68
69
70
71
72
73
74
75
76
77
78
79
80
    bounds = [BayesInfo.lb BayesInfo.ub] ;
    init_flag = 1;
end

% Get initial conditions for the state particles 
StateVectorMean = ReducedForm.StateVectorMean;
StateVectorVarianceSquareRoot = reduced_rank_cholesky(ReducedForm.StateVectorVariance)';
state_variance_rank = size(StateVectorVarianceSquareRoot,2);
StateVectors = bsxfun(@plus,StateVectorVarianceSquareRoot*randn(state_variance_rank,number_of_particles),StateVectorMean);
if pruning
    StateVectors_ = StateVectors;
end
        
81
% parameters for the Liu & West filter 
82
83
84
h_square = (3*liu_west_delta-1)/(2*liu_west_delta) ;
h_square = 1-h_square*h_square ;
small_a = sqrt(1-h_square) ;
85
86

% Initialization of parameter particles 
87
xparam = zeros(number_of_parameters,number_of_particles) ;
Frédéric Karamé's avatar
Frédéric Karamé committed
88
89
stderr = sqrt(bsxfun(@power,bounds(:,2)+bounds(:,1),2)/12)/100 ;
stderr = sqrt(bsxfun(@power,bounds(:,2)+bounds(:,1),2)/12)/50 ;
90
91
i = 1 ;
while i<=number_of_particles
Frédéric Karamé's avatar
Frédéric Karamé committed
92
93
    %candidate = start_param + .001*liu_west_chol_sigma_bar*randn(number_of_parameters,1) ;
    candidate = start_param + bsxfun(@times,stderr,randn(number_of_parameters,1)) ;
94
95
96
97
98
99
    if all(candidate(:) >= bounds(:,1)) && all(candidate(:) <= bounds(:,2))
        xparam(:,i) = candidate(:) ;
        i = i+1 ;
    end
end

100
101
%xparam = bsxfun(@plus,bounds(:,1),bsxfun(@times,(bounds(:,2)-bounds(:,1)),rand(number_of_parameters,number_of_particles))) ;

102
103
104
105
106
107
108
109
110
111
112
113
114
% Initialization of the weights of particles.
weights = ones(1,number_of_particles)/number_of_particles ;

% Initialization of the likelihood.
const_lik = log(2*pi)*number_of_observed_variables;
mean_xparam = zeros(number_of_parameters,sample_size) ;
median_xparam = zeros(number_of_parameters,sample_size) ;
std_xparam = zeros(number_of_parameters,sample_size) ;
lb95_xparam = zeros(number_of_parameters,sample_size) ;
ub95_xparam = zeros(number_of_parameters,sample_size) ;

%% The Online filter 
for t=1:sample_size
Frédéric Karamé's avatar
Frédéric Karamé committed
115
    disp(t)
116
117
118
119
    % Moments of parameters particles distribution 
    m_bar = xparam*(weights') ;
    temp = bsxfun(@minus,xparam,m_bar) ;
    sigma_bar = (bsxfun(@times,weights,temp))*(temp') ;
120
121
122
    if variance_update==1
        chol_sigma_bar = chol(h_square*sigma_bar)' ;
    end
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
    % Prediction (without shocks)
    ObservedVariables = zeros(number_of_observed_variables,number_of_particles) ;
    for i=1:number_of_particles
        % model resolution 
        [ys,trend_coeff,exit_flag,info,Model,DynareOptions,BayesInfo,DynareResults,ReducedForm] = ... 
            solve_model_for_online_filter(t,xparam(:,i),DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults) ;
        steadystate = ReducedForm.steadystate;
        state_variables_steady_state = ReducedForm.state_variables_steady_state;
        % Set local state space model (second-order approximation).
        constant = ReducedForm.constant;
        ghx  = ReducedForm.ghx;
        ghu  = ReducedForm.ghu;
        ghxx = ReducedForm.ghxx;
        ghuu = ReducedForm.ghuu;
        ghxu = ReducedForm.ghxu;
        % particle likelihood contribution  
        yhat = bsxfun(@minus,StateVectors(:,i),state_variables_steady_state);
        if pruning
            yhat_ = bsxfun(@minus,StateVectors_(:,i),state_variables_steady_state);
            [tmp, tmp_] = local_state_space_iteration_2(yhat,zeros(number_of_structural_innovations,1),ghx,ghu,constant,ghxx,ghuu,ghxu,yhat_,steadystate,DynareOptions.threads.local_state_space_iteration_2);
        else
            tmp = local_state_space_iteration_2(yhat,zeros(number_of_structural_innovations,1),ghx,ghu,constant,ghxx,ghuu,ghxu,DynareOptions.threads.local_state_space_iteration_2);
        end
        ObservedVariables(:,i) = tmp(mf1,:) ;
    end
    PredictedObservedMean = sum(bsxfun(@times,weights,ObservedVariables),2) ;
    PredictionError = bsxfun(@minus,Y(:,t),ObservedVariables);
    dPredictedObservedMean = bsxfun(@minus,ObservedVariables,PredictedObservedMean);
    PredictedObservedVariance = bsxfun(@times,weights,dPredictedObservedMean)*dPredictedObservedMean' + ReducedForm.H ;
    wtilde = exp(-.5*(const_lik+log(det(PredictedObservedVariance))+sum(PredictionError.*(PredictedObservedVariance\PredictionError),1))) ;
    % unormalized weights and observation likelihood contribution 
    tau_tilde = weights.*wtilde ;
    sum_tau_tilde = sum(tau_tilde) ;
    % particles selection 
    tau_tilde = tau_tilde/sum_tau_tilde ;
    indx = index_resample(0,tau_tilde',DynareOptions);
    StateVectors = StateVectors(:,indx) ;
    if pruning
        StateVectors_ = StateVectors_(:,indx) ;
    end
    xparam = bsxfun(@plus,(1-small_a).*m_bar,small_a.*xparam) ;
    xparam = xparam(:,indx) ;
    wtilde = wtilde(indx) ;
    % draw in the new distributions 
    i = 1 ;
    while i<=number_of_particles
169
        candidate = xparam(:,i) + chol_sigma_bar*randn(number_of_parameters,1) ;
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
        if all(candidate >= bounds(:,1)) && all(candidate <= bounds(:,2))
            xparam(:,i) = candidate ;
            % model resolution for new parameters particles
            [ys,trend_coeff,exit_flag,info,Model,DynareOptions,BayesInfo,DynareResults,ReducedForm] = ...
                solve_model_for_online_filter(t,xparam(:,i),DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults) ;
            steadystate = ReducedForm.steadystate;
            state_variables_steady_state = ReducedForm.state_variables_steady_state;
            % Set local state space model (second order approximation).
            constant = ReducedForm.constant;
            ghx  = ReducedForm.ghx;
            ghu  = ReducedForm.ghu;
            ghxx = ReducedForm.ghxx;
            ghuu = ReducedForm.ghuu;
            ghxu = ReducedForm.ghxu;
            % Get covariance matrices and structural shocks 
185
            epsilon = chol(ReducedForm.Q)'*randn(number_of_structural_innovations,1) ;
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
            % compute particles likelihood contribution 
            yhat = bsxfun(@minus,StateVectors(:,i),state_variables_steady_state);
            if pruning
                yhat_ = bsxfun(@minus,StateVectors_(:,i),state_variables_steady_state);
                [tmp, tmp_] = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,yhat_,steadystate,DynareOptions.threads.local_state_space_iteration_2);
                StateVectors_(:,i) = tmp_(mf0,:) ;
            else
                tmp = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,DynareOptions.threads.local_state_space_iteration_2);
            end
            StateVectors(:,i) = tmp(mf0,:) ;
            ObservedVariables(:,i) = tmp(mf1,:) ;
            i = i+1 ;
        end
    end
    PredictedObservedMean = sum(bsxfun(@times,weights,ObservedVariables),2) ;
    PredictionError = bsxfun(@minus,Y(:,t),ObservedVariables);
    dPredictedObservedMean = bsxfun(@minus,ObservedVariables,PredictedObservedMean);
    PredictedObservedVariance = bsxfun(@times,weights,dPredictedObservedMean)*dPredictedObservedMean' + ReducedForm.H ;
    lnw = exp(-.5*(const_lik+log(det(PredictedObservedVariance))+sum(PredictionError.*(PredictedObservedVariance\PredictionError),1)));
    % importance ratio 
    wtilde = lnw./wtilde ;
    % normalization 
    weights = wtilde/sum(wtilde);
209
210
    if (variance_update==1) && (neff(weights)<DynareOptions.particle.resampling.neff_threshold*sample_size)
        variance_update = 0 ;
211
    end
Frédéric Karamé's avatar
Frédéric Karamé committed
212
    % final resampling (advised)
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
    if second_resample==1 
        indx = index_resample(0,weights,DynareOptions);
        StateVectors = StateVectors(:,indx) ;
        if pruning
            StateVectors_ = StateVectors_(:,indx) ;
        end
        xparam = xparam(:,indx) ;
        weights = ones(1,number_of_particles)/number_of_particles ;
        mean_xparam(:,t) = mean(xparam,2);
        mat_var_cov = bsxfun(@minus,xparam,mean_xparam(:,t)) ;
        mat_var_cov = (mat_var_cov*mat_var_cov')/(number_of_particles-1) ;    
        std_xparam(:,t) = sqrt(diag(mat_var_cov)) ;
        for i=1:number_of_parameters
           temp = sortrows(xparam(i,:)') ;
           lb95_xparam(i,t) = temp(0.025*number_of_particles) ;
           median_xparam(i,t) = temp(0.5*number_of_particles) ;
           ub95_xparam(i,t) = temp(0.975*number_of_particles) ;
        end
    end
    if second_resample==0
        mean_xparam(:,t) = xparam*(weights') ;
        mat_var_cov = bsxfun(@minus,xparam,mean_xparam(:,t)) ;
        mat_var_cov = mat_var_cov*(bsxfun(@times,mat_var_cov,weights)') ;
        std_xparam(:,t) = sqrt(diag(mat_var_cov)) ;
        for i=1:number_of_parameters
           temp = sortrows([xparam(i,:)' weights'],1) ;
           cumulated_weights = cumsum(temp(:,2)) ;
           pass1=1 ;
           pass2=1 ;
           pass3=1 ;
           for j=1:number_of_particles
               if cumulated_weights(j)>0.025 && pass1==1 
                   lb95_xparam(i,t) = (temp(j-1,1)+temp(j,1))/2 ;
                   pass1 = 2 ;
               end
               if cumulated_weights(j)>0.5 && pass2==1
                   median_xparam(i,t) = (temp(j-1,1)+temp(j,1))/2 ;
                   pass2 = 2 ;
               end
               if cumulated_weights(j)>0.975 && pass3==1
                   ub95_xparam(i,t) = (temp(j-1,1)+temp(j,1))/2 ;
                   pass3 = 2 ;
               end
           end
        end
258
    end
259
    disp([lb95_xparam(:,t) mean_xparam(:,t) ub95_xparam(:,t)])
260
end
Frédéric Karamé's avatar
Frédéric Karamé committed
261
262
distrib_param = xparam ;
xparam = mean_xparam(:,sample_size) ;
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
std_param = std_xparam(:,sample_size) ;
lb_95 = lb95_xparam(:,sample_size) ;
ub_95 = ub95_xparam(:,sample_size) ;
median_param = median_xparam(:,sample_size) ;

%% Plot parameters trajectory  
TeX = DynareOptions.TeX;

[nbplt,nr,nc,lr,lc,nstar] = pltorg(number_of_parameters);

if TeX
    fidTeX = fopen([Model.fname '_param_traj.TeX'],'w');
    fprintf(fidTeX,'%% TeX eps-loader file generated by online_auxiliary_filter.m (Dynare).\n');
    fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
    fprintf(fidTeX,' \n');
end

z = 1:1:sample_size ;

for plt = 1:nbplt,
    if TeX
        NAMES = [];
        TeXNAMES = [];
    end
    hh = dyn_figure(DynareOptions,'Name','Parameters Trajectories');
    for k=1:min(nstar,length(xparam)-(plt-1)*nstar)
        subplot(nr,nc,k)
        kk = (plt-1)*nstar+k;
        [name,texname] = get_the_name(kk,TeX,Model,EstimatedParameters,DynareOptions);
        if TeX
            if isempty(NAMES)
                NAMES = name;
                TeXNAMES = texname;
            else
                NAMES = char(NAMES,name);
                TeXNAMES = char(TeXNAMES,texname);
            end
        end
        y = [mean_xparam(kk,:)' median_xparam(kk,:)' lb95_xparam(kk,:)' ub95_xparam(kk,:)' xparam(kk)*ones(sample_size,1)] ;
        plot(z,y);
        hold on
        title(name,'interpreter','none')
        hold off
        axis tight
        drawnow
    end
    dyn_saveas(hh,[ Model.fname '_param_traj' int2str(plt) ],DynareOptions);
    if TeX
        % TeX eps loader file
        fprintf(fidTeX,'\\begin{figure}[H]\n');
        for jj = 1:min(nstar,length(x)-(plt-1)*nstar)
            fprintf(fidTeX,'\\psfrag{%s}[1][][0.5][0]{%s}\n',deblank(NAMES(jj,:)),deblank(TeXNAMES(jj,:)));
        end
        fprintf(fidTeX,'\\centering \n');
        fprintf(fidTeX,'\\includegraphics[scale=0.5]{%s_ParamTraj%s}\n',Model.fname,int2str(plt));
        fprintf(fidTeX,'\\caption{Parameters trajectories.}');
        fprintf(fidTeX,'\\label{Fig:ParametersPlots:%s}\n',int2str(plt));
        fprintf(fidTeX,'\\end{figure}\n');
        fprintf(fidTeX,' \n');
    end
end

%% Plot Parameter Densities
number_of_grid_points = 2^9;      % 2^9 = 512 !... Must be a power of two.
bandwidth = 0;                    % Rule of thumb optimal bandwidth parameter.
kernel_function = 'gaussian';     % Gaussian kernel for Fast Fourier Transform approximation.  
for plt = 1:nbplt,
    if TeX
        NAMES = [];
        TeXNAMES = [];
    end
    hh = dyn_figure(DynareOptions,'Name','Parameters Densities');
    for k=1:min(nstar,length(xparam)-(plt-1)*nstar)
        subplot(nr,nc,k)
        kk = (plt-1)*nstar+k;
        [name,texname] = get_the_name(kk,TeX,Model,EstimatedParameters,DynareOptions);
        if TeX
            if isempty(NAMES)
                NAMES = name;
                TeXNAMES = texname;
            else
                NAMES = char(NAMES,name);
                TeXNAMES = char(TeXNAMES,texname);
            end
        end
Frédéric Karamé's avatar
Frédéric Karamé committed
348
349
        optimal_bandwidth = mh_optimal_bandwidth(distrib_param(kk,:)',number_of_particles,bandwidth,kernel_function);
        [density(:,1),density(:,2)] = kernel_density_estimate(distrib_param(kk,:)',number_of_grid_points,...
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
                                                          number_of_particles,optimal_bandwidth,kernel_function);
        plot(density(:,1),density(:,2));
        hold on
        title(name,'interpreter','none')
        hold off
        axis tight
        drawnow
    end
    dyn_saveas(hh,[ Model.fname '_param_density' int2str(plt) ],DynareOptions);
    if TeX
        % TeX eps loader file
        fprintf(fidTeX,'\\begin{figure}[H]\n');
        for jj = 1:min(nstar,length(x)-(plt-1)*nstar)
            fprintf(fidTeX,'\\psfrag{%s}[1][][0.5][0]{%s}\n',deblank(NAMES(jj,:)),deblank(TeXNAMES(jj,:)));
        end
        fprintf(fidTeX,'\\centering \n');
        fprintf(fidTeX,'\\includegraphics[scale=0.5]{%s_ParametersDensities%s}\n',Model.fname,int2str(plt));
        fprintf(fidTeX,'\\caption{ParametersDensities.}');
        fprintf(fidTeX,'\\label{Fig:ParametersDensities:%s}\n',int2str(plt));
        fprintf(fidTeX,'\\end{figure}\n');
        fprintf(fidTeX,' \n');
    end
end    
Frédéric Karamé's avatar
Frédéric Karamé committed
373