From ef22c716d3e23f749b82cb0e73392ddcf93c42b4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fr=C3=A9d=C3=A9ric=20Karam=C3=A9?= Date: Wed, 13 Jun 2018 18:43:04 +0200 Subject: [PATCH] Initial particles are drawn in the prior distribution + bug fixes. --- src/online_auxiliary_filter.m | 176 ++++++++++++++++++---------------- 1 file changed, 95 insertions(+), 81 deletions(-) diff --git a/src/online_auxiliary_filter.m b/src/online_auxiliary_filter.m index a14548559..8598a1318 100644 --- a/src/online_auxiliary_filter.m +++ b/src/online_auxiliary_filter.m @@ -43,7 +43,7 @@ persistent start_param sample_size number_of_observed_variables number_of_struct % Set seed for randn(). set_dynare_seed('default') ; pruning = DynareOptions.particle.pruning; -second_resample = 0 ; +second_resample = DynareOptions.particle.resampling.status.systematic ; variance_update = 1 ; % initialization of state particles @@ -75,25 +75,31 @@ if pruning end % parameters for the Liu & West filter -h_square = (3*liu_west_delta-1)/(2*liu_west_delta) ; -h_square = 1-h_square*h_square ; -small_a = sqrt(1-h_square) ; +small_a = (3*liu_west_delta-1)/(2*liu_west_delta) ; +b_square = 1-small_a*small_a ; % Initialization of parameter particles xparam = zeros(number_of_parameters,number_of_particles) ; -stderr = sqrt(bsxfun(@power,bounds.ub-bounds.lb,2)/12)/100 ; -stderr = sqrt(bsxfun(@power,bounds.ub-bounds.lb,2)/12)/50 ; +%stderr = sqrt(bsxfun(@power,bounds.ub-bounds.lb,2)/12)/100 ; +%stderr = sqrt(bsxfun(@power,bounds.ub-bounds.lb,2)/12)/50 ; %stderr = sqrt(bsxfun(@power,bounds.ub-bounds.lb,2)/12)/20 ; -i = 1 ; -while i<=number_of_particles - %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)) ; - if all(candidate(:) >= bounds.lb) && all(candidate(:) <= bounds.ub) - xparam(:,i) = candidate(:) ; - i = i+1 ; +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 end -end - +end %xparam = bsxfun(@plus,bounds(:,1),bsxfun(@times,(bounds(:,2)-bounds(:,1)),rand(number_of_parameters,number_of_particles))) ; % Initialization of the weights of particles. @@ -119,81 +125,87 @@ for t=1:sample_size temp = bsxfun(@minus,xparam,m_bar) ; sigma_bar = (bsxfun(@times,weights,temp))*(temp') ; if variance_update==1 - chol_sigma_bar = chol(h_square*sigma_bar)' ; + chol_sigma_bar = chol(b_square*sigma_bar)' ; end % Prediction (without shocks) + fore_xparam = bsxfun(@plus,(1-small_a).*m_bar,small_a.*xparam) ; tau_tilde = zeros(1,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 - PredictionError = bsxfun(@minus,Y(t,:)',tmp(mf1,:)); - % 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 - % particles selection - tau_tilde = tau_tilde/sum(tau_tilde) ; - indx = resample(0,tau_tilde',DynareOptions.particle); - StateVectors = StateVectors(:,indx) ; - if pruning - StateVectors_ = StateVectors_(:,indx) ; - end - xparam = bsxfun(@plus,(1-small_a).*m_bar,small_a.*xparam(:,indx)) ; - w_stage1 = weights(indx)./tau_tilde(indx) ; - % draw in the new distributions - wtilde = zeros(1,number_of_particles) ; - i = 1 ; - while i<=number_of_particles - candidate = xparam(:,i) + chol_sigma_bar*randn(number_of_parameters,1) ; - if all(candidate >= bounds.lb) && all(candidate <= bounds.ub) - 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) ; + solve_model_for_online_filter(t,fore_xparam(:,i),DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults) ; + if info==0 steadystate = ReducedForm.steadystate; state_variables_steady_state = ReducedForm.state_variables_steady_state; - % Set local state space model (second order approximation). + % 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 + % 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,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,yhat_,steadystate,DynareOptions.threads.local_state_space_iteration_2); - StateVectors_(:,i) = tmp_(mf0,:) ; + [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,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,DynareOptions.threads.local_state_space_iteration_2); + 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 - 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))); - i = i+1 ; + % 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 end end % normalization @@ -266,6 +278,10 @@ median_param = median_xparam(:,sample_size) ; TeX = DynareOptions.TeX; [nbplt,nr,nc,lr,lc,nstar] = pltorg(number_of_parameters); +nr = ceil(sqrt(number_of_parameters)) ; +nc = floor(sqrt(number_of_parameters)); +nbplt = 1 ; + if TeX fidTeX = fopen([Model.fname '_param_traj.tex'],'w'); @@ -282,10 +298,9 @@ for plt = 1:nbplt, TeXNAMES = []; end hh = dyn_figure(DynareOptions.nodisplay,'Name','Parameters Trajectories'); - for k=1:min(nstar,length(xparam)-(plt-1)*nstar) + for k=1:length(xparam) subplot(nr,nc,k) - kk = (plt-1)*nstar+k; - [name,texname] = get_the_name(kk,TeX,Model,EstimatedParameters,DynareOptions); + [name,texname] = get_the_name(k,TeX,Model,EstimatedParameters,DynareOptions); if TeX if isempty(NAMES) NAMES = name; @@ -295,7 +310,7 @@ for plt = 1:nbplt, 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)] ; + y = [mean_xparam(k,:)' median_xparam(k,:)' lb95_xparam(k,:)' ub95_xparam(k,:)' xparam(k)*ones(sample_size,1)] ; plot(z,y); hold on title(name,'interpreter','none') @@ -307,7 +322,7 @@ for plt = 1:nbplt, if TeX % TeX eps loader file fprintf(fidTeX,'\\begin{figure}[H]\n'); - for jj = 1:min(nstar,length(x)-(plt-1)*nstar) + for jj = 1:length(x) fprintf(fidTeX,'\\psfrag{%s}[1][][0.5][0]{%s}\n',deblank(NAMES(jj,:)),deblank(TeXNAMES(jj,:))); end fprintf(fidTeX,'\\centering \n'); @@ -329,10 +344,9 @@ for plt = 1:nbplt, TeXNAMES = []; end hh = dyn_figure(DynareOptions.nodisplay,'Name','Parameters Densities'); - for k=1:min(nstar,length(xparam)-(plt-1)*nstar) + for k=1:length(xparam) subplot(nr,nc,k) - kk = (plt-1)*nstar+k; - [name,texname] = get_the_name(kk,TeX,Model,EstimatedParameters,DynareOptions); + [name,texname] = get_the_name(k,TeX,Model,EstimatedParameters,DynareOptions); if TeX if isempty(NAMES) NAMES = name; @@ -342,8 +356,8 @@ for plt = 1:nbplt, TeXNAMES = char(TeXNAMES,texname); end end - 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,... + 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,... number_of_particles,optimal_bandwidth,kernel_function); plot(density(:,1),density(:,2)); hold on @@ -356,7 +370,7 @@ for plt = 1:nbplt, if TeX % TeX eps loader file fprintf(fidTeX,'\\begin{figure}[H]\n'); - for jj = 1:min(nstar,length(x)-(plt-1)*nstar) + for jj = 1:length(x) fprintf(fidTeX,'\\psfrag{%s}[1][][0.5][0]{%s}\n',deblank(NAMES(jj,:)),deblank(TeXNAMES(jj,:))); end fprintf(fidTeX,'\\centering \n');