online_auxiliary_filter.m 17 KB
Newer Older
1
function [xparam,std_param,lb_95,ub_95,median_param] = online_auxiliary_filter(xparam1,DynareDataset,dataset_info,DynareOptions,Model,EstimatedParameters,BayesInfo,bounds,DynareResults)
Stéphane Adjemian's avatar
Stéphane Adjemian committed
2

3
% Liu & West particle filter = auxiliary particle filter including Liu & West filter on parameters.
4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
%
% 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.

Houtan Bastani's avatar
Houtan Bastani committed
24
% Copyright (C) 2013-2017 Dynare Team
25 26 27 28 29 30 31 32 33 34 35 36 37 38 39
%
% 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/>.
40
persistent Y init_flag mf0 mf1 number_of_particles number_of_parameters liu_west_delta liu_west_chol_sigma_bar
41 42 43 44 45
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
second_resample = DynareOptions.particle.resampling.status.systematic ;
47 48 49
variance_update = 1 ;

% initialization of state particles
Stéphane Adjemian's avatar
Stéphane Adjemian committed
50 51
[ys,trend_coeff,exit_flag,info,Model,DynareOptions,BayesInfo,DynareResults,ReducedForm] = ...
    solve_model_for_online_filter(1,xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults) ;
52 53 54 55 56 57 58 59

% 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 ;
60
    sample_size = size(Y,1);
61 62 63
    number_of_observed_variables = length(mf1);
    number_of_structural_innovations = length(ReducedForm.Q);
    liu_west_delta = DynareOptions.particle.liu_west_delta ;
Stéphane Adjemian's avatar
Stéphane Adjemian committed
64
    start_param = xparam1 ;
65 66 67
    init_flag = 1;
end

Stéphane Adjemian's avatar
Stéphane Adjemian committed
68
% Get initial conditions for the state particles
69
StateVectorMean = ReducedForm.StateVectorMean;
70
StateVectorVarianceSquareRoot = chol(ReducedForm.StateVectorVariance)';
71 72 73 74 75
state_variance_rank = size(StateVectorVarianceSquareRoot,2);
StateVectors = bsxfun(@plus,StateVectorVarianceSquareRoot*randn(state_variance_rank,number_of_particles),StateVectorMean);
if pruning
    StateVectors_ = StateVectors;
end
Stéphane Adjemian's avatar
Stéphane Adjemian committed
76 77

% parameters for the Liu & West filter
78 79
small_a = (3*liu_west_delta-1)/(2*liu_west_delta) ;
b_square = 1-small_a*small_a ;
80

Stéphane Adjemian's avatar
Stéphane Adjemian committed
81
% Initialization of parameter particles
82
xparam = zeros(number_of_parameters,number_of_particles) ;
83 84
%stderr = sqrt(bsxfun(@power,bounds.ub-bounds.lb,2)/12)/100 ;
%stderr = sqrt(bsxfun(@power,bounds.ub-bounds.lb,2)/12)/50 ;
85
%stderr = sqrt(bsxfun(@power,bounds.ub-bounds.lb,2)/12)/20 ;
86 87 88 89 90 91 92 93 94 95 96 97 98 99 100
bounds = prior_bounds(BayesInfo,DynareOptions.prior_trunc); %reset bounds as lb and ub must only be operational during mode-finding
prior_draw(BayesInfo,DynareOptions.prior_trunc);
for i=1:number_of_particles
    info = 1;
    while info==1
        %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)) ;
        candidate = prior_draw()';
        if all(candidate(:) >= bounds.lb) && all(candidate(:) <= bounds.ub)
            [ys,trend_coeff,exit_flag,info,Model,DynareOptions,BayesInfo,DynareResults,ReducedForm] = ...
                solve_model_for_online_filter(1,candidate(:),DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults) ;
            if info==0
                xparam(:,i) = candidate(:) ;
            end
        end 
101
    end
102
end 
103 104 105 106 107 108 109 110 111 112 113 114 115
%xparam = bsxfun(@plus,bounds(:,1),bsxfun(@times,(bounds(:,2)-bounds(:,1)),rand(number_of_parameters,number_of_particles))) ;

% 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) ;

Stéphane Adjemian's avatar
Stéphane Adjemian committed
116
%% The Online filter
117
for t=1:sample_size
118 119 120 121 122
    if t>1
        fprintf('\nSubsample with %s first observations.\n\n', int2str(t))
    else
        fprintf('\nSubsample with only the first observation.\n\n', int2str(t))
    end
Stéphane Adjemian's avatar
Stéphane Adjemian committed
123
    % Moments of parameters particles distribution
124 125 126 127
    m_bar = xparam*(weights') ;
    temp = bsxfun(@minus,xparam,m_bar) ;
    sigma_bar = (bsxfun(@times,weights,temp))*(temp') ;
    if variance_update==1
128
        chol_sigma_bar = chol(b_square*sigma_bar)' ;
129 130
    end
    % Prediction (without shocks)
131
    fore_xparam = bsxfun(@plus,(1-small_a).*m_bar,small_a.*xparam) ;
132
    tau_tilde = zeros(1,number_of_particles) ;
133
    for i=1:number_of_particles
Stéphane Adjemian's avatar
Stéphane Adjemian committed
134 135
        % model resolution
        [ys,trend_coeff,exit_flag,info,Model,DynareOptions,BayesInfo,DynareResults,ReducedForm] = ...
136 137
            solve_model_for_online_filter(t,fore_xparam(:,i),DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults) ;
        if info==0
138 139
            steadystate = ReducedForm.steadystate;
            state_variables_steady_state = ReducedForm.state_variables_steady_state;
140
            % Set local state space model (second-order approximation).
141 142 143 144 145 146
            constant = ReducedForm.constant;
            ghx  = ReducedForm.ghx;
            ghu  = ReducedForm.ghu;
            ghxx = ReducedForm.ghxx;
            ghuu = ReducedForm.ghuu;
            ghxu = ReducedForm.ghxu;
147
            % particle likelihood contribution
148 149 150
            yhat = bsxfun(@minus,StateVectors(:,i),state_variables_steady_state);
            if pruning
                yhat_ = bsxfun(@minus,StateVectors_(:,i),state_variables_steady_state);
151
                [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);
152
            else
153
                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);
154
            end
155
            PredictionError = bsxfun(@minus,Y(t,:)',tmp(mf1,:));
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
            % Replace Gaussian density with a Student density with 3 degrees of
            % freedom for fat tails.
            z = sum(PredictionError.*(ReducedForm.H\PredictionError),1) ;
            tau_tilde(i) = weights(i).*(tpdf(z,3*ones(size(z)))+1e-99) ;
            %tau_tilde(i) = weights(i).*exp(-.5*(const_lik+log(det(ReducedForm.H))+sum(PredictionError.*(ReducedForm.H\PredictionError),1))) ;
        end 
    end
    % particles selection
    tau_tilde = tau_tilde/sum(tau_tilde)  ;
    indx = resample(0,tau_tilde',DynareOptions.particle);
    StateVectors = StateVectors(:,indx) ;
    xparam = fore_xparam(:,indx);	
    if pruning
        StateVectors_ = StateVectors_(:,indx) ;
    end
    w_stage1 = weights(indx)./tau_tilde(indx) ;
    % draw in the new distributions
    wtilde = zeros(1,number_of_particles) ;
    for i=1:number_of_particles
        info=1 ;
        while info==1 
            candidate = xparam(:,i) + chol_sigma_bar*randn(number_of_parameters,1) ;
            if all(candidate >= bounds.lb) && all(candidate <= bounds.ub)
	                % model resolution for new parameters particles
                [ys,trend_coeff,exit_flag,info,Model,DynareOptions,BayesInfo,DynareResults,ReducedForm] = ...
                    solve_model_for_online_filter(t,candidate,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults) ;
                if info==0
                    xparam(:,i) = candidate ;
                    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
                    epsilon = chol(ReducedForm.Q)'*randn(number_of_structural_innovations,1) ;
                    % 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,:) ;
                    PredictionError = bsxfun(@minus,Y(t,:)',tmp(mf1,:));
                    wtilde(i) = w_stage1(i)*exp(-.5*(const_lik+log(det(ReducedForm.H))+sum(PredictionError.*(ReducedForm.H\PredictionError),1)));
                end 
            end
209 210
        end
    end
Stéphane Adjemian's avatar
Stéphane Adjemian committed
211
    % normalization
212 213 214 215
    weights = wtilde/sum(wtilde);
    if (variance_update==1) && (neff(weights)<DynareOptions.particle.resampling.threshold*sample_size)
        variance_update = 0 ;
    end
216
    % final resampling (not advised)
Stéphane Adjemian's avatar
Stéphane Adjemian committed
217
    if second_resample==1
218
        indx = resample(0,weights,DynareOptions.particle);
219 220 221 222 223 224 225 226
        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)) ;
Stéphane Adjemian's avatar
Stéphane Adjemian committed
227
        mat_var_cov = (mat_var_cov*mat_var_cov')/(number_of_particles-1) ;
228 229
        std_xparam(:,t) = sqrt(diag(mat_var_cov)) ;
        for i=1:number_of_parameters
Stéphane Adjemian's avatar
Stéphane Adjemian committed
230 231 232 233
            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) ;
234 235 236 237 238 239 240 241
        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
Stéphane Adjemian's avatar
Stéphane Adjemian committed
242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260
            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) ;
                    pass1 = 2 ;
                end
                if cumulated_weights(j)>=0.5 && pass2==1
                    median_xparam(i,t) = temp(j,1) ;
                    pass2 = 2 ;
                end
                if cumulated_weights(j)>=0.975 && pass3==1
                    ub95_xparam(i,t) = temp(j,1) ;
                    pass3 = 2 ;
                end
            end
261 262
        end
    end
263 264 265 266 267 268
    str = sprintf(' Lower Bound (95%%) \t Mean \t\t\t Upper Bound (95%%)');
    for l=1:size(xparam,1)
        str = sprintf('%s\n %5.4f \t\t %7.5f \t\t %5.4f', str, lb95_xparam(l,t), mean_xparam(l,t), ub95_xparam(l,t));
    end
    disp([str])
    disp('')
269 270 271 272 273 274 275 276
end
distrib_param = xparam ;
xparam = mean_xparam(:,sample_size) ;
std_param = std_xparam(:,sample_size) ;
lb_95 = lb95_xparam(:,sample_size) ;
ub_95 = ub95_xparam(:,sample_size) ;
median_param = median_xparam(:,sample_size) ;

Stéphane Adjemian's avatar
Stéphane Adjemian committed
277
%% Plot parameters trajectory
278 279 280
TeX = DynareOptions.TeX;

[nbplt,nr,nc,lr,lc,nstar] = pltorg(number_of_parameters);
281 282 283 284
nr = ceil(sqrt(number_of_parameters)) ;
nc = floor(sqrt(number_of_parameters));
nbplt = 1 ;

285 286

if TeX
287
    fidTeX = fopen([Model.fname '_param_traj.tex'],'w');
288 289 290 291 292 293 294 295 296 297 298 299
    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
Houtan Bastani's avatar
Houtan Bastani committed
300
    hh = dyn_figure(DynareOptions.nodisplay,'Name','Parameters Trajectories');
301
    for k=1:length(xparam)
302
        subplot(nr,nc,k)
303
        [name,texname] = get_the_name(k,TeX,Model,EstimatedParameters,DynareOptions);
304 305 306 307 308 309 310 311 312
        if TeX
            if isempty(NAMES)
                NAMES = name;
                TeXNAMES = texname;
            else
                NAMES = char(NAMES,name);
                TeXNAMES = char(TeXNAMES,texname);
            end
        end
313
        y = [mean_xparam(k,:)' median_xparam(k,:)' lb95_xparam(k,:)' ub95_xparam(k,:)' xparam(k)*ones(sample_size,1)] ;
314 315 316 317 318 319 320
        plot(z,y);
        hold on
        title(name,'interpreter','none')
        hold off
        axis tight
        drawnow
    end
Houtan Bastani's avatar
Houtan Bastani committed
321
    dyn_saveas(hh,[ Model.fname '_param_traj' int2str(plt) ],DynareOptions.nodisplay,DynareOptions.graph_format);
322 323 324
    if TeX
        % TeX eps loader file
        fprintf(fidTeX,'\\begin{figure}[H]\n');
325
        for jj = 1:length(x)
326 327 328 329 330 331 332 333 334 335 336 337 338 339
            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.
Stéphane Adjemian's avatar
Stéphane Adjemian committed
340
kernel_function = 'gaussian';     % Gaussian kernel for Fast Fourier Transform approximation.
341 342 343 344 345
for plt = 1:nbplt,
    if TeX
        NAMES = [];
        TeXNAMES = [];
    end
Houtan Bastani's avatar
Houtan Bastani committed
346
    hh = dyn_figure(DynareOptions.nodisplay,'Name','Parameters Densities');
347
    for k=1:length(xparam)
348
        subplot(nr,nc,k)
349
        [name,texname] = get_the_name(k,TeX,Model,EstimatedParameters,DynareOptions);
350 351 352 353 354 355 356 357 358
        if TeX
            if isempty(NAMES)
                NAMES = name;
                TeXNAMES = texname;
            else
                NAMES = char(NAMES,name);
                TeXNAMES = char(TeXNAMES,texname);
            end
        end
359 360
        optimal_bandwidth = mh_optimal_bandwidth(distrib_param(k,:)',number_of_particles,bandwidth,kernel_function);
        [density(:,1),density(:,2)] = kernel_density_estimate(distrib_param(k,:)',number_of_grid_points,...
361 362 363 364 365 366 367 368
                                                          number_of_particles,optimal_bandwidth,kernel_function);
        plot(density(:,1),density(:,2));
        hold on
        title(name,'interpreter','none')
        hold off
        axis tight
        drawnow
    end
Houtan Bastani's avatar
Houtan Bastani committed
369
    dyn_saveas(hh,[ Model.fname '_param_density' int2str(plt) ],DynareOptions.nodisplay,DynareOptions.graph_format);
370 371 372
    if TeX
        % TeX eps loader file
        fprintf(fidTeX,'\\begin{figure}[H]\n');
373
        for jj = 1:length(x)
374 375 376 377 378 379 380 381 382
            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
Stéphane Adjemian's avatar
Stéphane Adjemian committed
383
end