Commit 8820c63f authored by Frédéric Karamé's avatar Frédéric Karamé
Browse files

Bug fixes.

parent 5db18090
...@@ -37,7 +37,7 @@ function [LIK,lik] = conditional_particle_filter(ReducedForm,Y,start,DynareOptio ...@@ -37,7 +37,7 @@ function [LIK,lik] = conditional_particle_filter(ReducedForm,Y,start,DynareOptio
% %
% NOTES % NOTES
% The vector "lik" is used to evaluate the jacobian of the likelihood. % The vector "lik" is used to evaluate the jacobian of the likelihood.
% Copyright (C) 2009-2013 Dynare Team % Copyright (C) 2009-2010 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
...@@ -116,7 +116,7 @@ for t=1:sample_size ...@@ -116,7 +116,7 @@ for t=1:sample_size
if (strcmp(DynareOptions.particle.resampling.status,'generic') && neff(SampleWeights)<DynareOptions.particle.resampling.neff_threshold*sample_size ) || ... if (strcmp(DynareOptions.particle.resampling.status,'generic') && neff(SampleWeights)<DynareOptions.particle.resampling.neff_threshold*sample_size ) || ...
strcmp(DynareOptions.particle.resampling.status,'systematic') strcmp(DynareOptions.particle.resampling.status,'systematic')
ks = ks + 1 ; ks = ks + 1 ;
StateParticles = resample(StateParticles',SampleWeights,DynareOptions)'; StateParticles = resample(StateParticles',SampleWeights',DynareOptions)';
SampleWeights = ones(1,number_of_particles)/number_of_particles ; SampleWeights = ones(1,number_of_particles)/number_of_particles ;
end end
end end
......
...@@ -19,7 +19,7 @@ function IncrementalWeights = gaussian_densities(obs,mut_t,sqr_Pss_t_t,st_t_1,sq ...@@ -19,7 +19,7 @@ function IncrementalWeights = gaussian_densities(obs,mut_t,sqr_Pss_t_t,st_t_1,sq
% %
% NOTES % NOTES
% The vector "lik" is used to evaluate the jacobian of the likelihood. % The vector "lik" is used to evaluate the jacobian of the likelihood.
% Copyright (C) 2009-2012 Dynare Team % Copyright (C) 2009-2010 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
......
...@@ -59,7 +59,6 @@ function new_particles = multivariate_smooth_resampling(particles,weights) ...@@ -59,7 +59,6 @@ function new_particles = multivariate_smooth_resampling(particles,weights)
% stephane DOT adjemian AT univ DASH lemans DOT fr % stephane DOT adjemian AT univ DASH lemans DOT fr
number_of_particles = length(weights); number_of_particles = length(weights);
weights = weights' ;
number_of_states = size(particles,2); number_of_states = size(particles,2);
[P,D] = eig(particles'*(bsxfun(@times,1/number_of_particles,particles))) ; [P,D] = eig(particles'*(bsxfun(@times,1/number_of_particles,particles))) ;
D = diag(D) ; D = diag(D) ;
......
...@@ -41,7 +41,6 @@ persistent Y init_flag mf0 mf1 bounds number_of_particles number_of_parameters l ...@@ -41,7 +41,6 @@ persistent Y init_flag mf0 mf1 bounds number_of_particles number_of_parameters l
persistent start_param sample_size number_of_observed_variables number_of_structural_innovations persistent start_param sample_size number_of_observed_variables number_of_structural_innovations
% Set seed for randn(). % Set seed for randn().
%set_dynare_seed('mt19937ar',1234) ;
set_dynare_seed('default') ; set_dynare_seed('default') ;
pruning = DynareOptions.particle.pruning; pruning = DynareOptions.particle.pruning;
second_resample = 1 ; second_resample = 1 ;
...@@ -62,11 +61,10 @@ if isempty(init_flag) ...@@ -62,11 +61,10 @@ if isempty(init_flag)
number_of_observed_variables = length(mf1); number_of_observed_variables = length(mf1);
number_of_structural_innovations = length(ReducedForm.Q); number_of_structural_innovations = length(ReducedForm.Q);
liu_west_delta = DynareOptions.particle.liu_west_delta ; liu_west_delta = DynareOptions.particle.liu_west_delta ;
liu_west_chol_sigma_bar = DynareOptions.particle.liu_west_chol_sigma_bar*eye(number_of_parameters) ; %liu_west_chol_sigma_bar = DynareOptions.particle.liu_west_chol_sigma_bar*eye(number_of_parameters) ;
%start_param = xparam1 ; start_param = xparam1 ;
% Conditions initiales %liu_west_chol_sigma_bar = sqrt(bsxfun(@times,eye(number_of_parameters),BayesInfo.p2)) ;
%liu_west_chol_sigma_bar = bsxfun(@times,eye(number_of_parameters),BayesInfo.p2) ; %start_param = BayesInfo.p1 ;
start_param = BayesInfo.p1 ;
bounds = [BayesInfo.lb BayesInfo.ub] ; bounds = [BayesInfo.lb BayesInfo.ub] ;
init_flag = 1; init_flag = 1;
end end
...@@ -87,11 +85,12 @@ small_a = sqrt(1-h_square) ; ...@@ -87,11 +85,12 @@ small_a = sqrt(1-h_square) ;
% Initialization of parameter particles % Initialization of parameter particles
xparam = zeros(number_of_parameters,number_of_particles) ; xparam = zeros(number_of_parameters,number_of_particles) ;
stderr = sqrt(bsxfun(@power,bounds(:,2)+bounds(:,1),2)/12)/1000 ; stderr = sqrt(bsxfun(@power,bounds(:,2)+bounds(:,1),2)/12)/100 ;
stderr = sqrt(bsxfun(@power,bounds(:,2)+bounds(:,1),2)/12)/50 ;
i = 1 ; i = 1 ;
while i<=number_of_particles while i<=number_of_particles
candidate = start_param + 10*liu_west_chol_sigma_bar*randn(number_of_parameters,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 = start_param + bsxfun(@times,stderr,randn(number_of_parameters,1)) ;
if all(candidate(:) >= bounds(:,1)) && all(candidate(:) <= bounds(:,2)) if all(candidate(:) >= bounds(:,1)) && all(candidate(:) <= bounds(:,2))
xparam(:,i) = candidate(:) ; xparam(:,i) = candidate(:) ;
i = i+1 ; i = i+1 ;
...@@ -113,7 +112,7 @@ ub95_xparam = zeros(number_of_parameters,sample_size) ; ...@@ -113,7 +112,7 @@ ub95_xparam = zeros(number_of_parameters,sample_size) ;
%% The Online filter %% The Online filter
for t=1:sample_size for t=1:sample_size
disp(t) disp(t)
% Moments of parameters particles distribution % Moments of parameters particles distribution
m_bar = xparam*(weights') ; m_bar = xparam*(weights') ;
temp = bsxfun(@minus,xparam,m_bar) ; temp = bsxfun(@minus,xparam,m_bar) ;
...@@ -210,7 +209,7 @@ for t=1:sample_size ...@@ -210,7 +209,7 @@ for t=1:sample_size
if (variance_update==1) && (neff(weights)<DynareOptions.particle.resampling.neff_threshold*sample_size) if (variance_update==1) && (neff(weights)<DynareOptions.particle.resampling.neff_threshold*sample_size)
variance_update = 0 ; variance_update = 0 ;
end end
% final resampling (advised) % final resampling (advised)
if second_resample==1 if second_resample==1
indx = index_resample(0,weights,DynareOptions); indx = index_resample(0,weights,DynareOptions);
StateVectors = StateVectors(:,indx) ; StateVectors = StateVectors(:,indx) ;
...@@ -259,7 +258,8 @@ for t=1:sample_size ...@@ -259,7 +258,8 @@ for t=1:sample_size
end end
disp([lb95_xparam(:,t) mean_xparam(:,t) ub95_xparam(:,t)]) disp([lb95_xparam(:,t) mean_xparam(:,t) ub95_xparam(:,t)])
end end
xparam = mean_xparam(:,sample_size) ; distrib_param = xparam ;
xparam = mean_xparam(:,sample_size) ;
std_param = std_xparam(:,sample_size) ; std_param = std_xparam(:,sample_size) ;
lb_95 = lb95_xparam(:,sample_size) ; lb_95 = lb95_xparam(:,sample_size) ;
ub_95 = ub95_xparam(:,sample_size) ; ub_95 = ub95_xparam(:,sample_size) ;
...@@ -345,8 +345,8 @@ for plt = 1:nbplt, ...@@ -345,8 +345,8 @@ for plt = 1:nbplt,
TeXNAMES = char(TeXNAMES,texname); TeXNAMES = char(TeXNAMES,texname);
end end
end end
optimal_bandwidth = mh_optimal_bandwidth(xparam(kk,:)',number_of_particles,bandwidth,kernel_function); optimal_bandwidth = mh_optimal_bandwidth(distrib_param(kk,:)',number_of_particles,bandwidth,kernel_function);
[density(:,1),density(:,2)] = kernel_density_estimate(xparam(kk,:)',number_of_grid_points,... [density(:,1),density(:,2)] = kernel_density_estimate(distrib_param(kk,:)',number_of_grid_points,...
number_of_particles,optimal_bandwidth,kernel_function); number_of_particles,optimal_bandwidth,kernel_function);
plot(density(:,1),density(:,2)); plot(density(:,1),density(:,2));
hold on hold on
...@@ -370,3 +370,4 @@ for plt = 1:nbplt, ...@@ -370,3 +370,4 @@ for plt = 1:nbplt,
fprintf(fidTeX,' \n'); fprintf(fidTeX,' \n');
end end
end end
\ No newline at end of file
...@@ -54,6 +54,7 @@ function resampled_particles = resample(particles,weights,DynareOptions) ...@@ -54,6 +54,7 @@ function resampled_particles = resample(particles,weights,DynareOptions)
% AUTHOR(S) frederic DOT karame AT univ DASH evry DOT fr % AUTHOR(S) frederic DOT karame AT univ DASH evry DOT fr
% stephane DOT adjemian AT univ DASH lemans DOT fr % stephane DOT adjemian AT univ DASH lemans DOT fr
switch DynareOptions.particle.resampling.method1 switch DynareOptions.particle.resampling.method1
case 'residual' case 'residual'
if strcmpi(DynareOptions.particle.resampling.method2,'kitagawa') if strcmpi(DynareOptions.particle.resampling.method2,'kitagawa')
......
...@@ -161,7 +161,7 @@ for t=1:sample_size ...@@ -161,7 +161,7 @@ for t=1:sample_size
StateVectors = temp(:,1:number_of_state_variables)' ; StateVectors = temp(:,1:number_of_state_variables)' ;
StateVectors_ = temp(:,number_of_state_variables+1:2*number_of_state_variables)'; StateVectors_ = temp(:,number_of_state_variables+1:2*number_of_state_variables)';
else else
StateVectors = resample(tmp(mf0,:)',weights,DynareOptions)'; StateVectors = resample(tmp(mf0,:)',weights',DynareOptions)';
end end
weights = ones(1,number_of_particles)/number_of_particles; weights = ones(1,number_of_particles)/number_of_particles;
elseif strcmp(DynareOptions.particle.resampling.status,'none') elseif strcmp(DynareOptions.particle.resampling.status,'none')
......
...@@ -168,8 +168,7 @@ end ...@@ -168,8 +168,7 @@ end
offset = EstimatedParameters.nvx; offset = EstimatedParameters.nvx;
if EstimatedParameters.nvn if EstimatedParameters.nvn
for i=1:EstimatedParameters.nvn for i=1:EstimatedParameters.nvn
k = EstimatedParameters.var_endo(i,1); H(i,i) = xparam1(i+offset)*xparam1(i+offset);
H(k,k) = xparam1(i+offset)*xparam1(i+offset);
end end
offset = offset+EstimatedParameters.nvn; offset = offset+EstimatedParameters.nvn;
else else
......
...@@ -35,12 +35,12 @@ steady; ...@@ -35,12 +35,12 @@ steady;
//disp(oo_.mean) ; //disp(oo_.mean) ;
estimated_params; estimated_params;
alp, uniform_pdf,,, 0.0001, 1; alp, uniform_pdf,,, 0.0001, 0.5;
bet, uniform_pdf,,, 0.75, 0.999; bet, uniform_pdf,,, 0.0001, 0.99;
tet, uniform_pdf,,, 0.0001, 1; tet, uniform_pdf,,, 0.0001, 1;
tau, uniform_pdf,,, 0.0001, 100; tau, uniform_pdf,,, 0.0001, 100;
delt, uniform_pdf,,, 0.0001, 0.05; delt, uniform_pdf,,, 0.0001, 0.05;
rho, uniform_pdf,,, 0.0001, 0.999; rho, uniform_pdf,,, 0.8, 0.99;
stderr e_a, uniform_pdf,,, 0.00001, 0.1; stderr e_a, uniform_pdf,,, 0.00001, 0.1;
stderr y, uniform_pdf,,, 0.00001, 0.1; stderr y, uniform_pdf,,, 0.00001, 0.1;
stderr l, uniform_pdf,,, 0.00001, 0.1; stderr l, uniform_pdf,,, 0.00001, 0.1;
...@@ -49,11 +49,11 @@ end; ...@@ -49,11 +49,11 @@ end;
estimated_params_init; estimated_params_init;
alp, 0.4; alp, 0.4;
bet, 0.99; bet, 0.97;
tet, 0.357 ; tet, 0.357 ;
tau, 50; tau, 50;
delt, 0.02; delt, 0.02;
rho, 0.95; rho, 0.9 ;
stderr e_a, .035; stderr e_a, .035;
stderr y, .0175;//.00158; stderr y, .0175;//.00158;
stderr l, .00312;//.0011; stderr l, .00312;//.0011;
...@@ -66,28 +66,42 @@ varobs y l i ; ...@@ -66,28 +66,42 @@ varobs y l i ;
//options_.gstep(2) = .1; //options_.gstep(2) = .1;
options_.particle.status = 1; options_.particle.status = 1;
options_.particle.algorithm = 'sequential_importance_particle_filter';
options_.particle.initialization = 1; options_.particle.initialization = 1;
options_.particle.pruning = 1; options_.particle.pruning = 0;
options_.particle.number_of_particles = 2000; options_.particle.number_of_particles = 1000 ;
options_.particle.resampling.status = 'systematic'; options_.particle.resampling.status = 'systematic';
options_.particle.resampling.neff_threshold = .1;
//options_.particle.resampling.method1 = 'traditional' ;
//options_.particle.resampling.method1 = 'residual' ;
options_.particle.resampling.method1 = 'smooth' ;
options_.particle.reampling.method2 = 'kitagawa' ;
//options_.particle.resampling.method2 = 'stratified' ;
options_.particle.resampling.number_of_partitions = 1;
options_.particle.resampling.neff_threshold = .5;
set_dynare_threads('local_state_space_iteration_2',3); set_dynare_threads('local_state_space_iteration_2',3);
options_.particle.algorithm = 'sequential_importance_particle_filter'; options_.particle.algorithm = 'sequential_importance_particle_filter';
//options_.particle.algorithm = 'auxiliary_particle_filter'; //options_.particle.algorithm = 'auxiliary_particle_filter';
//options_.particle.algorithm = 'gaussian_mixture_filter'; //options_.particle.algorithm = 'gaussian_mixture_filter';
//options_.particle.algorithm = 'each_gaussian_filter'; //options_.particle.algorithm = 'conditional_particle_filter';
//options_.particle.algorithm = 'gaussian_filter'; //options_.particle.algorithm = 'gaussian_filter';
//options_.particle.IS_approximation_method = 'quadrature' ; //options_.particle.IS_approximation_method = 'quadrature' ;
//options_.particle.IS_approximation_method = 'cubature' ; options_.particle.IS_approximation_method = 'cubature' ;
//options_.particle.IS_approximation_method = 'unscented' ; //options_.particle.IS_approximation_method = 'unscented' ;
//options_.particle.approximation_method = 'quadrature' ; //options_.particle.approximation_method = 'quadrature' ;
//options_.particle.approximation_method = 'cubature' ; options_.particle.approximation_method = 'cubature' ;
//options_.particle.approximation_method = 'unscented' ; //options_.particle.approximation_method = 'unscented' ;
//options_.particle.approximation_method = 'MonteCarlo' ; //options_.particle.approximation_method = 'MonteCarlo' ;
estimation(datafile=data_risky_perturb2,nograph,order=2,nobs=100,mh_replic=0,mode_compute=7,mode_check); //options_.mh_posterior_mode_estimation=1 ;
// online
options_.particle.liu_west_delta = 0.9 ;
options_.mode_check_node_number = 250 ;
estimation(datafile=data_risky_perturb3,nograph,order=2,nobs=100,mh_replic=0,mode_compute=7,mode_check);
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment