Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • giovanma/dynare
  • giorgiomas/dynare
  • Vermandel/dynare
  • Dynare/dynare
  • normann/dynare
  • MichelJuillard/dynare
  • wmutschl/dynare
  • FerhatMihoubi/dynare
  • sebastien/dynare
  • lnsongxf/dynare
  • rattoma/dynare
  • CIMERS/dynare
  • FredericKarame/dynare
  • SumuduK/dynare
  • MinjeJeon/dynare
  • camilomrch/dynare
  • DoraK/dynare
  • avtishin/dynare
  • selma/dynare
  • claudio_olguin/dynare
  • jeffjiang07/dynare
  • EthanSystem/dynare
  • stepan-a/dynare
  • wjgatt/dynare
  • JohannesPfeifer/dynare
  • gboehl/dynare
26 results
Show changes
Commits on Source (22)
Showing
with 308 additions and 261 deletions
......@@ -8,8 +8,8 @@ Software requirements
=====================
Packaged versions of Dynare are available for Windows (10 and 11), several
GNU/Linux distributions (Debian, Ubuntu, Linux Mint, Arch Linux), macOS (13
Ventura), and FreeBSD. Dynare should work on other systems, but some
GNU/Linux distributions (Debian, Ubuntu, Linux Mint, Arch Linux), macOS (15
“Sequoia”), and FreeBSD. Dynare should work on other systems, but some
compilation steps are necessary in that case.
In order to run Dynare, you need one of the following:
......
......@@ -15,7 +15,7 @@ function write(M_)
% REMARKS
% - The trends are assumed to be multiplicative.
% Copyright © 2019-2023 Dynare Team
% Copyright © 2019-2024 Dynare Team
%
% This file is part of Dynare.
%
......@@ -84,37 +84,40 @@ fprintf(fid, 'y = z(1:%u);\n\n', M_.endo_nbr);
fprintf(fid, 'g = z(%u:%u);\n', M_.endo_nbr+1, 2*M_.endo_nbr);
% Define the point where the dynamic model is to be evaluated.
fprintf(fid, 'Y = zeros(%u, 1);\n', 2*(n0+n1+n2));
% In period t, then in period t+1
fprintf(fid, 'Y0 = NaN(%u, 1);\n', 3*M_.endo_nbr);
fprintf(fid, 'Y1 = NaN(%u, 1);\n', 3*M_.endo_nbr);
for i=1:length(I0) % period t equations, lagged variables.
if I0(i)
fprintf(fid, 'Y(%u) = y(%u);\n', I0(i), i);
fprintf(fid, 'Y0(%u) = y(%u);\n', i+(purely_forward_model*M_.endo_nbr), i);
end
end
for i=1:length(I1) % period t equations, current variables.
if I1(i)
fprintf(fid, 'Y(%u) = y(%u)*g(%u);\n', I1(i), i, i);
fprintf(fid, 'Y0(%u) = y(%u)*g(%u);\n', i+M_.endo_nbr+(purely_forward_model*M_.endo_nbr), i, i);
end
end
for i=1:length(I2) % period t equations, leaded variables.
if I2(i)
fprintf(fid, 'Y(%u) = y(%u)*g(%u)*g(%u);\n', I2(i), i, i, i);
fprintf(fid, 'Y0(%u) = y(%u)*g(%u)*g(%u);\n', i+2*M_.endo_nbr, i, i, i);
end
end
for i=1:length(I0) % period t+1 equations lagged variables.
if I0(i)
fprintf(fid, 'Y(%u) = y(%u)*g(%u);\n', n0+n1+n2+I0(i), i, i);
fprintf(fid, 'Y1(%u) = y(%u)*g(%u);\n', i+(purely_forward_model*M_.endo_nbr), i, i);
end
end
for i=1:length(I1) % period t+1 equations current variables.
if I1(i)
fprintf(fid, 'Y(%u) = y(%u)*g(%u)*g(%u);\n', n0+n1+n2+I1(i), i, i, i);
fprintf(fid, 'Y1(%u) = y(%u)*g(%u)*g(%u);\n', i+M_.endo_nbr+(purely_forward_model*M_.endo_nbr), i, i, i);
end
end
for i=1:length(I2) % period t+1 equations leaded variables.
if I2(i)
fprintf(fid, 'Y(%u) = y(%u)*g(%u)*g(%u)*g(%u);\n', n0+n1+n2+I2(i), i, i, i, i);
fprintf(fid, 'Y1(%u) = y(%u)*g(%u)*g(%u)*g(%u);\n', i+2*M_.endo_nbr, i, i, i, i);
end
end
fprintf(fid, '\n');
% Define the vector of parameters.
......@@ -131,19 +134,36 @@ fprintf(fid, 'F = NaN(%u, 1);\n', 2*M_.endo_nbr);
fprintf(fid, 'x = zeros(1, %u);\n\n', M_.exo_nbr);
% Evaluate the residuals and jacobian of the dynamic model in periods t and t+1.
fprintf(fid, '[F(1:%u), T0_order, T0] = %s.sparse.dynamic_resid(Y0, x, p, y);\n', M_.endo_nbr, M_.fname);
fprintf(fid, '[F(%u:%u), T1_order, T1] = %s.sparse.dynamic_resid(Y1, x, p, y);\n', M_.endo_nbr+1, 2*M_.endo_nbr, M_.fname);
fprintf(fid, 'if nargout>1\n');
fprintf(fid, ' sparse_rowval = [');
fprintf(fid, '%u ', M_.dynamic_g1_sparse_rowval);
fprintf(fid, '];\n');
fprintf(fid, ' sparse_colval = [');
fprintf(fid, '%u ', M_.dynamic_g1_sparse_colval);
fprintf(fid, '];\n');
fprintf(fid, ' sparse_colptr = [');
fprintf(fid, '%u ', M_.dynamic_g1_sparse_colptr);
fprintf(fid, '];\n');
fprintf(fid, ' J0 = %s.sparse.dynamic_g1(Y0, x, p, y, sparse_rowval, sparse_colval, sparse_colptr, T0_order, T0);\n', M_.fname);
fprintf(fid, ' J1 = %s.sparse.dynamic_g1(Y1, x, p, y, sparse_rowval, sparse_colval, sparse_colptr, T1_order, T1);\n', M_.fname);
% Transform back the Jacobians J0 and J1 in the legacy format (non-sparse)
% NB: it is probably possible to simplify the rest of this file, but maintaining
% decent performance does not seem straightforward.
lli = find(M_.lead_lag_incidence');
if purely_forward_model
lli = lli + M_.endo_nbr;
end
fprintf(fid, ' lli = [');
fprintf(fid, '%u ', lli);
fprintf(fid, '];\n');
fprintf(fid, ' J = zeros(%u, %u);\n', 2*M_.endo_nbr, n0+n1+n2+M_.endo_nbr);
fprintf(fid, ' [F(1:%u), tmp] = %s.dynamic(Y(1:%u), x, p, y, 1);\n', M_.endo_nbr, M_.fname, n0+n1+n2);
fprintf(fid, ' J(1:%u,1:%u) = tmp(:,1:%u);\n', M_.endo_nbr, n0+n1+n2, n0+n1+n2);
fprintf(fid, ' [F(%u:%u), tmp] = %s.dynamic(Y(1+%u:%u), x, p, y, 1);\n', M_.endo_nbr+1, 2*M_.endo_nbr, M_.fname, n0+n1+n2, 2*(n0+n1+n2));
fprintf(fid, ' J(%u:%u,1:%u) = tmp(:,1:%u);\n', M_.endo_nbr+1, 2*M_.endo_nbr, n0+n1+n2, n0+n1+n2);
fprintf(fid, 'else\n');
fprintf(fid, ' F(1:%u) = %s.dynamic(Y(1:%u), x, p, y, 1);\n', M_.endo_nbr, M_.fname, n0+n1+n2);
fprintf(fid, ' F(%u:%u) = %s.dynamic(Y(1+%u:%u), x, p, y, 1);\n', M_.endo_nbr+1, 2*M_.endo_nbr, M_.fname, n0+n1+n2, 2*(n0+n1+n2));
fprintf(fid, 'end\n\n');
fprintf(fid, ' J(1:%u,1:%u) = full(J0(:,lli));\n', M_.endo_nbr, n0+n1+n2);
fprintf(fid, ' J(%u:%u,1:%u) = full(J1(:,lli));\n', M_.endo_nbr+1, 2*M_.endo_nbr, n0+n1+n2);
% Compute the jacobian if required.
fprintf(fid, 'if nargout>1\n');
fprintf(fid, ' JAC = zeros(%u,%u);\n', 2*M_.endo_nbr, 2*M_.endo_nbr);
% Compute the derivatives of the first block of equations (period t)
......@@ -279,7 +299,7 @@ else
end
% Compute the derivatives of the second block of equations (period t+1)
% with respect to the endogenous variables.
% with respect to the growth factors.
if purely_backward_model || purely_forward_model
for i=M_.eq_nbr+1:2*M_.eq_nbr
for j=1:M_.endo_nbr
......@@ -328,4 +348,4 @@ end
fprintf(fid, ' JAC(:,%u:%u) = J(:,%u:%u);\n', M_.endo_nbr+1, 2*M_.endo_nbr, n0+n1+n2+1, n0+n1+n2+M_.endo_nbr);
fprintf(fid,'end\n');
fclose(fid);
\ No newline at end of file
fclose(fid);
......@@ -13,7 +13,7 @@ function run(json)
% SPECIAL REQUIREMENTS
% none
% Copyright © 2019-2023 Dynare Team
% Copyright © 2019-2024 Dynare Team
%
% This file is part of Dynare.
%
......@@ -103,7 +103,6 @@ if ~isempty(jm.anticipated_transitory_shocks)
'periods', s.start_date:s.end_date, ...
'value', s.value)];
end
M_.exo_det_length = 0;
end
%% Make unanticipated shock map
......
......@@ -13,7 +13,7 @@ function read(json)
% SPECIAL REQUIREMENTS
% none
% Copyright © 2019-2020 Dynare Team
% Copyright © 2019-2024 Dynare Team
%
% This file is part of Dynare.
%
......@@ -36,7 +36,6 @@ global M_ options_ oo_
jm = loadjson_(json, 'SimplifyCell', 1);
data2json=struct();
M_.exo_det_length = 0;
for nshocks = 1:length(jm.stochasticshocksdescription)
covartype=jm.stochasticshocksdescription{nshocks}.shockattributevalue;
thisshock=(jm.stochasticshocksdescription{nshocks}.shockindex)+1;
......
......@@ -13,7 +13,7 @@ function [h_minus_1, h, h_plus_1, h_exo, resid] = get_deriv(M_, ys_)
% - h_exo [N by N_exo] derivative matrix with respect to exogenous variables
% - resid [N by 1] vector of residuals
% Copyright © 2021 Dynare Team
% Copyright © 2021-2024 Dynare Team
%
% This file is part of Dynare.
%
......@@ -30,40 +30,16 @@ function [h_minus_1, h, h_plus_1, h_exo, resid] = get_deriv(M_, ys_)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
x = zeros(M_.maximum_lag + M_.maximum_lead + 1,M_.exo_nbr);
dyn_endo_ss = repmat(ys_, 3, 1);
x = zeros(M_.exo_nbr);
iyv = M_.lead_lag_incidence';
iyr0 = find(iyv(:)) ;
z=repmat(ys_,1,M_.maximum_lag + M_.maximum_lead + 1);
[resid, T_order, T] = feval([M_.fname '.sparse.dynamic_resid'], dyn_endo_ss, x, M_.params, ys_);
[resid,g1]=feval([M_.fname,'.dynamic'],z(iyr0),x, M_.params, ys_, M_.maximum_exo_lag+1);
g1 = feval([M_.fname '.sparse.dynamic_g1'], dyn_endo_ss, x, M_.params, ys_, ...
M_.dynamic_g1_sparse_rowval, M_.dynamic_g1_sparse_colval, ...
M_.dynamic_g1_sparse_colptr, T_order, T);
% Initialize matrices
h_minus_1=zeros(M_.endo_nbr);
h = h_minus_1;
h_plus_1 = h_minus_1;
% build h_minus_1
if M_.maximum_lag
lag_columns=find(iyv(:,1));
n_lag_columns=length(lag_columns);
h_minus_1(:,lag_columns) = g1(:,1:n_lag_columns);
else
n_lag_columns=0;
end
% build h
contemporaneous_columns=find(iyv(:,M_.maximum_lag+1));
n_contemporaneous_columns = length(contemporaneous_columns);
h(:,contemporaneous_columns) = g1(:,1+n_lag_columns:n_lag_columns+n_contemporaneous_columns);
%build h_plus_1
if M_.maximum_lead
lead_columns=find(iyv(:,end));
n_lead_columns = length(lead_columns);
h_plus_1(:,lead_columns) = g1(:,n_lag_columns+n_contemporaneous_columns+1:n_lag_columns+n_contemporaneous_columns+n_lead_columns);
else
n_lead_columns=0;
end
h_exo =g1(:,n_lag_columns+n_contemporaneous_columns+n_lead_columns+1:end);
h_minus_1 = full(g1(:,1:M_.endo_nbr));
h = full(g1(:,M_.endo_nbr + (1:M_.endo_nbr)));
h_plus_1 = full(g1(:,2*M_.endo_nbr + (1:M_.endo_nbr)));
h_exo = full(g1(:,3*M_.endo_nbr+1:end));
......@@ -12,7 +12,7 @@ function options_ = default_option_values(M_)
% SPECIAL REQUIREMENTS
% none
% Copyright © 2018-2023 Dynare Team
% Copyright © 2018-2024 Dynare Team
%
% This file is part of Dynare.
%
......@@ -288,9 +288,6 @@ particle.resampling.number_of_partitions = 200;
particle.mixture_state_variables = 5 ;
particle.mixture_structural_shocks = 1 ;
particle.mixture_measurement_shocks = 1 ;
% Online approach
particle.liu_west_delta = 0.99 ;
particle.liu_west_max_resampling_tries = 5000;
% Options for setting the weights in conditional particle filters.
particle.cpf_weights_method.amisanotristani = true;
particle.cpf_weights_method.murrayjonesparslow = false;
......@@ -372,8 +369,8 @@ options_.initial_period = NaN; %dates(1,1);
options_.no_init_estimation_check_first_obs=false;
options_.dataset.file = [];
options_.dataset.series = [];
options_.dataset.firstobs = dates();
options_.dataset.lastobs = dates();
options_.dataset.first_obs = dates();
options_.dataset.last_obs = dates();
options_.dataset.nobs = NaN;
options_.dataset.xls_sheet = [];
options_.dataset.xls_range = [];
......@@ -514,6 +511,11 @@ options_.posterior_sampler_options.dsmh.particles = 20000 ;
options_.posterior_sampler_options.dsmh.alpha0 = 0.2 ;
options_.posterior_sampler_options.dsmh.alpha1 = 0.3 ;
options_.posterior_sampler_options.dsmh.tau = 10 ;
% Liu and West online Sampler
options_.posterior_sampler_options.online.particles= 5000 ;
options_.posterior_sampler_options.online.liu_west_delta = 0.99 ;
options_.posterior_sampler_options.online.liu_west_max_resampling_tries = 5000 ;
options_.posterior_sampler_options.online.second_resampling = false;
options_.trace_plot_ma = 200;
options_.mh_autocorrelation_function_size = 30;
......@@ -657,6 +659,12 @@ particleswarm.UseParallel = false;
particleswarm.UseVectorized = false;
options_.particleswarm = particleswarm;
% One step
one_step.gradient = [] ;
one_step.hessian = [] ;
one_step.outer_product = 0 ;
options_.one_step= one_step ;
% prior analysis
options_.prior_mc = 20000;
......
......@@ -53,6 +53,7 @@ p = {'/../contrib/ms-sbvar/TZcode/MatlabFiles/' ; ...
'/ep/' ; ...
'/estimation/'; ...
'/estimation/smc/'; ...
'/estimation/online/'; ...
'/estimation/resampler/'; ...
'/kalman/' ; ...
'/kalman/likelihood' ; ...
......
......@@ -51,6 +51,18 @@ if issmc(options_)
else
draws = posterior.tlogpostkernel;
end
elseif isonline(options_)
% Load draws from the posterior distribution
posterior = load(sprintf('%s/online/parameters_particles_final.mat', dname));
if column>0 || strcmp(column,'all')
if strcmp(column,'all')
draws = transpose(posterior.param);
else
draws = transpose(posterior.param(column,:));
end
else
draws=NaN(size(posterior.param,2),1);
end
elseif isdime(options_)
posterior = load(sprintf('%s%s%s%schains.mat', dname, filesep(), 'dime', filesep()));
tune = posterior.tune;
......
......@@ -36,6 +36,12 @@ if issmc(options_)
mean = sum(posterior.particles, 2)/length(posterior.tlogpostkernel);
% Compute the posterior covariance
variance = (posterior.particles-mean)*(posterior.particles-mean)'/length(posterior.tlogpostkernel);
elseif isonline(options_)
posterior = load(sprintf('%s/online/parameters_particles_final.mat', M_.dname));
% Compute the posterior mean
mean = sum(posterior.param, 2)/size(posterior.param, 2);
% Compute the posterior covariance
variance = (posterior.param-mean)*(posterior.param-mean)'/size(posterior.param, 2);
elseif isdime(options_)
posterior = load(sprintf('%s%s%s%schains.mat', M_.dname, filesep(), 'dime', filesep()));
tune = posterior.tune;
......
......@@ -55,9 +55,11 @@ skipline(2)
disp('ESTIMATION RESULTS')
skipline()
if ishssmc(options_)
if ishssmc(options_)
dprintf('Log data density is %f.', oo_.MarginalDensity.hssmc);
hpd_draws = round((1-options_.mh_conf_sig)*num_draws);
elseif isonline(options_)
hpd_draws = num_draws;
elseif isdime(options_)
hpd_draws = round((1-options_.mh_conf_sig)*num_draws);
else
......@@ -87,7 +89,7 @@ if np
disp(tit2)
ip = nvx+nvn+ncx+ncn+1;
for i=1:np
if options_.mh_replic || (options_.load_mh_file && ~options_.load_results_after_load_mh) || ishssmc(options_)
if options_.mh_replic || (options_.load_mh_file && ~options_.load_results_after_load_mh) || ishssmc(options_) || isonline(options_)
draws = getalldraws(ip);
[post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(draws, options_.mh_conf_sig);
name = bayestopt_.name{ip};
......
......@@ -216,10 +216,12 @@ else
localVars.ifil2=ifil2;
% which files have to be copied to run remotely
NamFileInput(1,:) = {'',[M_.fname '.static.m']};
NamFileInput(2,:) = {'',[M_.fname '.dynamic.m']};
NamFileInput(1,:) = {'',[M_.fname '.sparse.static_resid.m']};
NamFileInput(2,:) = {'',[M_.fname '.sparse.static_g1.m']};
NamFileInput(3,:) = {'',[M_.fname '.sparse.dynamic_resid.m']};
NamFileInput(4,:) = {'',[M_.fname '.sparse.dynamic_g1.m']};
if M_.set_auxiliary_variables
NamFileInput(3,:) = {'',[M_.fname '.set_auxiliary_variables.m']};
NamFileInput(5,:) = {'',[M_.fname '.set_auxiliary_variables.m']};
end
if options_.steadystate_flag
if options_.steadystate_flag == 1
......@@ -436,4 +438,4 @@ if ~options_.nograph && ~options_.no_graph.posterior
% END parallel code!
end
fprintf('%s: Posterior IRFs, done!\n',dispString);
\ No newline at end of file
fprintf('%s: Posterior IRFs, done!\n',dispString);
......@@ -45,6 +45,8 @@ if ~issmc(options_)
else
if ishssmc(options_)
[MetropolisFolder, info] = CheckPath('hssmc',M_.dname);
elseif isonline(options_)
[MetropolisFolder, info] = CheckPath('online',M_.dname);
elseif isdime(options_)
[MetropolisFolder, info] = CheckPath('dime',M_.dname);
else
......@@ -76,6 +78,10 @@ else
% Load draws from the posterior distribution
pfiles = dir(sprintf('%s/hssmc/particles-*.mat', M_.dname));
mhdate = get_date_of_a_file(sprintf('%s/hssmc/particles-%u-%u.mat', M_.dname, length(pfiles), length(pfiles)));
elseif isonline(options_)
% Load draws from the posterior distribution
pfiles = dir(sprintf('%s/online/parameters_particles_final.mat', M_.dname));
mhdate = get_date_of_a_file(sprintf('%s/online/parameters_particles_final.mat', M_.dname, length(pfiles), length(pfiles)));
elseif isdime(options_)
mhdate = get_date_of_a_file(sprintf('%s%s%s%schains.mat', M_.dname, filesep(), 'dime', filesep()));
end
......
......@@ -374,125 +374,150 @@ if init
posterior_sampler_options.WR=[];
end
case 'hssmc'
case 'hssmc'
% default options
posterior_sampler_options = add_fields_(posterior_sampler_options, options_.posterior_sampler_options.hssmc);
% default options
posterior_sampler_options = add_fields_(posterior_sampler_options, options_.posterior_sampler_options.hssmc);
% user defined options
if ~isempty(options_.posterior_sampler_options.sampling_opt)
options_list = read_key_value_string(options_.posterior_sampler_options.sampling_opt);
for i=1:rows(options_list)
switch options_list{i,1}
% user defined options
if ~isempty(options_.posterior_sampler_options.sampling_opt)
options_list = read_key_value_string(options_.posterior_sampler_options.sampling_opt);
for i=1:rows(options_list)
switch options_list{i,1}
case 'target'
posterior_sampler_options.target = options_list{i,2};
posterior_sampler_options.target = options_list{i,2};
case 'steps'
posterior_sampler_options.steps = options_list{i,2};
posterior_sampler_options.steps = options_list{i,2};
case 'scale'
posterior_sampler_options.scale = options_list{i,2};
posterior_sampler_options.scale = options_list{i,2};
case 'particles'
posterior_sampler_options.particles = options_list{i,2};
posterior_sampler_options.particles = options_list{i,2};
case 'lambda'
posterior_sampler_options.lambda = options_list{i,2};
posterior_sampler_options.lambda = options_list{i,2};
otherwise
warning(['hssmc: Unknown option (' options_list{i,1} ')!'])
end
end
end
options_.mode_compute = 0;
options_.cova_compute = 0;
options_.mh_replic = 0;
options_.mh_posterior_mode_estimation = false;
case 'dime_mcmc'
% default options
posterior_sampler_options = add_fields_(posterior_sampler_options, options_.posterior_sampler_options.dime);
% user defined options
if ~isempty(options_.posterior_sampler_options.sampling_opt)
options_list = read_key_value_string(options_.posterior_sampler_options.sampling_opt);
for i=1:rows(options_list)
switch options_list{i,1}
warning(['hssmc: Unknown option (' options_list{i,1} ')!'])
end
end
end
options_.mode_compute = 0;
options_.cova_compute = 0;
options_.mh_replic = 0;
options_.mh_posterior_mode_estimation = false;
case 'dime_mcmc'
% default options
posterior_sampler_options = add_fields_(posterior_sampler_options, options_.posterior_sampler_options.dime);
if ~isempty(options_.posterior_sampler_options.sampling_opt)
options_list = read_key_value_string(options_.posterior_sampler_options.sampling_opt);
for i=1:rows(options_list)
switch options_list{i,1}
case 'nchain'
posterior_sampler_options.nchain = options_list{i,2};
posterior_sampler_options.nchain = options_list{i,2};
case 'niter'
posterior_sampler_options.niter = options_list{i,2};
posterior_sampler_options.niter = options_list{i,2};
case 'parallel'
posterior_sampler_options.parallel = options_list{i,2};
posterior_sampler_options.parallel = options_list{i,2};
case 'tune'
posterior_sampler_options.tune = options_list{i,2};
posterior_sampler_options.tune = options_list{i,2};
case 'aimh_prob'
posterior_sampler_options.aimh_prob = options_list{i,2};
posterior_sampler_options.aimh_prob = options_list{i,2};
case 'sigma'
posterior_sampler_options.sigma = options_list{i,2};
posterior_sampler_options.sigma = options_list{i,2};
case 'gamma'
posterior_sampler_options.gamma = options_list{i,2};
posterior_sampler_options.gamma = options_list{i,2};
case 'df_proposal_dist'
posterior_sampler_options.df_proposal_dist = options_list{i,2};
posterior_sampler_options.df_proposal_dist = options_list{i,2};
case 'rho'
posterior_sampler_options.rho = options_list{i,2};
posterior_sampler_options.rho = options_list{i,2};
otherwise
warning(['dime: Unknown option (' options_list{i,1} ')!'])
end
end
end
if ~isfield(posterior_sampler_options,'tune')
posterior_sampler_options.tune = ceil(50000 / posterior_sampler_options.nchain);
dprintf('dime: setting number of ensemble iterations to keep (`tune`) to %d.', posterior_sampler_options.tune)
end
options_.mode_compute = 0;
options_.cova_compute = 0;
options_.mh_replic = 0;
options_.mh_posterior_mode_estimation = false;
warning(['dime: Unknown option (' options_list{i,1} ')!'])
end
end
end
if ~isfield(posterior_sampler_options,'tune')
posterior_sampler_options.tune = ceil(50000 / posterior_sampler_options.nchain);
dprintf('dime: setting number of ensemble iterations to keep (`tune`) to %d.', posterior_sampler_options.tune)
end
case 'online'
% default options
posterior_sampler_options = add_fields_(posterior_sampler_options, options_.posterior_sampler_options.online);
% user defined options
if ~isempty(options_.posterior_sampler_options.sampling_opt)
options_list = read_key_value_string(options_.posterior_sampler_options.sampling_opt);
for i=1:rows(options_list)
switch options_list{i,1}
case 'particles'
posterior_sampler_options.particles = options_list{i,2};
case 'liu_west_delta'
posterior_sampler_options.liu_west_delta = options_list{i,2};
case 'liu_west_max_resampling_tries'
posterior_sampler_options.liu_west_max_resampling_tries = options_list{i,2};
case 'second_resampling'
posterior_sampler_options.second_resampling = options_list{i,2};
otherwise
warning(['online: Unknown option (' options_list{i,1} ')!'])
end
end
end
options_.mode_compute = 0;
options_.cova_compute = 0;
options_.mh_replic = 0;
options_.mh_posterior_mode_estimation = false;
case 'dsmh'
% default options
posterior_sampler_options = add_fields_(posterior_sampler_options, options_.posterior_sampler_options.dsmh);
% user defined options
if ~isempty(options_.posterior_sampler_options.sampling_opt)
options_list = read_key_value_string(options_.posterior_sampler_options.sampling_opt);
for i=1:rows(options_list)
switch options_list{i,1}
case 'proposal_distribution'
if ~(strcmpi(options_list{i,2}, 'rand_multivariate_student') || ...
strcmpi(options_list{i,2}, 'rand_multivariate_normal'))
error(['initial_estimation_checks:: the proposal_distribution option to estimation takes either ' ...
'rand_multivariate_student or rand_multivariate_normal as options']);
else
posterior_sampler_options.proposal_distribution=options_list{i,2};
end
case 'student_degrees_of_freedom'
if options_list{i,2} <= 0
error('initial_estimation_checks:: the student_degrees_of_freedom takes a positive integer argument');
else
posterior_sampler_options.student_degrees_of_freedom=options_list{i,2};
end
case 'save_tmp_file'
posterior_sampler_options.save_tmp_file = options_list{i,2};
case 'number_of_particles'
posterior_sampler_options.particles = options_list{i,2};
otherwise
warning(['rwmh_sampler: Unknown option (' options_list{i,1} ')!'])
end
end
end
options_.mode_compute = 0;
options_.cova_compute = 0;
options_.mh_replic = 0;
options_.mh_posterior_mode_estimation = true;
% default options
posterior_sampler_options = add_fields_(posterior_sampler_options, options_.posterior_sampler_options.dsmh);
% user defined options
if ~isempty(options_.posterior_sampler_options.sampling_opt)
options_list = read_key_value_string(options_.posterior_sampler_options.sampling_opt);
for i=1:rows(options_list)
switch options_list{i,1}
case 'proposal_distribution'
if ~(strcmpi(options_list{i,2}, 'rand_multivariate_student') || ...
strcmpi(options_list{i,2}, 'rand_multivariate_normal'))
error(['initial_estimation_checks:: the proposal_distribution option to estimation takes either ' ...
'rand_multivariate_student or rand_multivariate_normal as options']);
else
posterior_sampler_options.proposal_distribution=options_list{i,2};
end
case 'student_degrees_of_freedom'
if options_list{i,2} <= 0
error('initial_estimation_checks:: the student_degrees_of_freedom takes a positive integer argument');
else
posterior_sampler_options.student_degrees_of_freedom=options_list{i,2};
end
case 'save_tmp_file'
posterior_sampler_options.save_tmp_file = options_list{i,2};
case 'number_of_particles'
posterior_sampler_options.particles = options_list{i,2};
otherwise
warning(['rwmh_sampler: Unknown option (' options_list{i,1} ')!'])
end
end
end
options_.mode_compute = 0;
options_.cova_compute = 0;
options_.mh_replic = 0;
options_.mh_posterior_mode_estimation = true;
otherwise
error('check_posterior_sampler_options:: Unknown posterior_sampling_method option %s ',posterior_sampler_options.posterior_sampling_method);
end
otherwise
error('check_posterior_sampler_options:: Unknown posterior_sampling_method option %s ',posterior_sampler_options.posterior_sampling_method);
end
return
return
end
end % if init
% here are all samplers requiring a proposal distribution
if ~strcmp(posterior_sampler_options.posterior_sampling_method,'slice')
......
......@@ -238,11 +238,7 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation &
optimizer_vec = [options_.mode_compute;num2cell(options_.additional_optimizer_steps)];
for optim_iter = 1:length(optimizer_vec)
current_optimizer = optimizer_vec{optim_iter};
if current_optimizer~=11
[xparam1, fval, ~, hh, options_, Scale, new_rat_hess_info] = dynare_minimize_objective(objective_function,xparam1,current_optimizer,options_,[bounds.lb bounds.ub],bayestopt_.name,bayestopt_,hh,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_.dr, oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
else
[xparam1, fval, ~, hh, options_, Scale, new_rat_hess_info] = dynare_minimize_objective(objective_function,xparam1,current_optimizer,options_,[bounds.lb bounds.ub],bayestopt_.name,bayestopt_,hh,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_);
end
[xparam1, fval, ~, hh, options_, Scale, new_rat_hess_info] = dynare_minimize_objective(objective_function,xparam1,current_optimizer,options_,[bounds.lb bounds.ub],bayestopt_.name,bayestopt_,hh,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_.dr, oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
fprintf('\nFinal value of minus the log posterior (or likelihood):%f \n', fval);
if isnumeric(current_optimizer)
......@@ -272,7 +268,7 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation &
hh = hessian(penalized_objective_function, xparam1, options_.gstep, objective_function, fval, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, bounds,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
else
hh = hessian(objective_function, xparam1, options_.gstep, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, bounds,oo_.dr,oo_.steady_state,oo_.exo_steady_state,oo_.exo_det_steady_state);
end
end
hh = reshape(hh, nx, nx);
elseif isnumeric(current_optimizer) && isequal(current_optimizer,5)
% other numerical hessian options available with optimizer
......@@ -399,7 +395,11 @@ end
% Run SMC sampler.
%
if ishssmc(options_)
if isonline(options_)
[posterior_sampler_options, options_, bayestopt_] = check_posterior_sampler_options([], M_.fname, M_.dname, options_, bounds, bayestopt_);
options_.posterior_sampler_options.current_options = posterior_sampler_options;
online_auxiliary_filter(xparam1, dataset_, options_, M_, estim_params_, bayestopt_, oo_);
elseif ishssmc(options_)
[posterior_sampler_options, options_, bayestopt_] = check_posterior_sampler_options([], M_.fname, M_.dname, options_, bounds, bayestopt_);
options_.posterior_sampler_options.current_options = posterior_sampler_options;
oo_.MarginalDensity.hssmc = hssmc(objective_function, bounds, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, oo_);
......@@ -472,7 +472,7 @@ if issmc(options_) || (any(bayestopt_.pshape>0) && options_.mh_replic) || (any(
oo_.lprob = trace_plot_dime(options_, M_);
end
% Estimation of the marginal density from the Mh draws:
if ishssmc(options_) || isdime(options_) || options_.mh_replic || (options_.load_mh_file && ~options_.load_results_after_load_mh)
if ishssmc(options_) || isonline(options_) || isdime(options_) || options_.mh_replic || (options_.load_mh_file && ~options_.load_results_after_load_mh)
if ~issmc(options_)
[~, oo_] = marginal_density(M_, options_, estim_params_, oo_, bayestopt_);
end
......
......@@ -36,6 +36,8 @@ else
folder_name='hssmc';
elseif isdime(options_)
folder_name='dime';
elseif isonline(options_)
folder_name='online';
else
error('get_posterior_folder_name:: case should not happen. Please contact the developers')
end
......
function bool = isonline(options_)
% Copyright © 2024 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 <https://www.gnu.org/licenses/>.
bool = false;
if isfield(options_, 'posterior_sampler_options')
if strcmp(options_.posterior_sampler_options.posterior_sampling_method, 'online')
bool = true;
end
end
function [pmean, pmode, pmedian, pstdev, p025, p975, covariance] = online_auxiliary_filter(xparam1, dataset_, options_, M_, estim_params_, bayestopt_, oo_)
% [pmean, pmode, pmedian, pstdev, p025, p975, covariance] = online_auxiliary_filter(xparam1, dataset_, options_, M_, estim_params_, bayestopt_, oo_)
function online_auxiliary_filter(xparam1, dataset_, options_, M_, estim_params_, bayestopt_, oo_)
% online_auxiliary_filter(xparam1, dataset_, options_, M_, estim_params_, bayestopt_, oo_)
% Liu & West particle filter = auxiliary particle filter including Liu & West filter on parameters.
%
% INPUTS
......@@ -12,16 +12,7 @@ function [pmean, pmode, pmedian, pstdev, p025, p975, covariance] = online_auxili
% - bayestopt_ [struct] Prior definition.
% - oo_ [struct] Results.
%
% OUTPUTS
% - pmean [double] n×1 vector, mean of the particles at the end of the sample (for the parameters).
% - pmode [double] n×1 vector, mode of the particles at the end of the sample (for the parameters).
% - pmedian [double] n×1 vector, median of the particles at the end of the sample (for the parameters).
% - pstdev [double] n×1 vector, st. dev. of the particles at the end of the sample (for the parameters).
% - p025 [double] n×1 vector, 2.5 percent of the particles are below p025(i) for i=1,…,n.
% - p975 [double] n×1 vector, 97.5 percent of the particles are below p975(i) for i=1,…,n.
% - covariance [double] n×n matrix, covariance of the particles at the end of the sample.
% Copyright © 2013-2023 Dynare Team
% Copyright © 2013-2024 Dynare Team
%
% This file is part of Dynare.
%
......@@ -41,8 +32,11 @@ function [pmean, pmode, pmedian, pstdev, p025, p975, covariance] = online_auxili
% Set seed for randn().
options_ = set_dynare_seed_local_options(options_,'default');
pruning = options_.particle.pruning;
second_resample = options_.particle.resampling.status.systematic;
variance_update = true;
online_opt = options_.posterior_sampler_options.current_options;
% Set location for the simulated particles.
SimulationFolder = CheckPath('online', M_.dname);
bounds = prior_bounds(bayestopt_, options_.prior_trunc); % Reset bounds as lb and ub must only be operational during mode-finding
......@@ -52,13 +46,13 @@ bounds = prior_bounds(bayestopt_, options_.prior_trunc); % Reset bounds as lb an
order = options_.order;
mf0 = ReducedForm.mf0;
mf1 = ReducedForm.mf1;
number_of_particles = options_.particle.number_of_particles;
number_of_particles = online_opt.particles;
number_of_parameters = size(xparam1,1);
Y = dataset_.data;
sample_size = size(Y,1);
number_of_observed_variables = length(mf1);
number_of_structural_innovations = length(ReducedForm.Q);
liu_west_delta = options_.particle.liu_west_delta;
liu_west_delta = online_opt.liu_west_delta;
% Get initial conditions for the state particles
StateVectorMean = ReducedForm.StateVectorMean;
......@@ -99,7 +93,6 @@ 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);
mode_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);
......@@ -194,7 +187,6 @@ for t=1:sample_size
% Replace Gaussian density with a Student density with 3 degrees of freedom for fat tails.
z = sum(PredictionError.*(ReducedForm.H\PredictionError), 1) ;
ddl = 3 ;
%tau_tilde(i) = weights(i).*(tpdf(z,3*ones(size(z)))+1e-99) ;
tau_tilde(i) = weights(i)*(exp(gammaln((ddl + 1) / 2) - gammaln(ddl/2))/(sqrt(ddl*pi)*(1 + (z^2)/ddl)^((ddl + 1)/2))+1e-99) ;
else
tau_tilde(i) = 0 ;
......@@ -202,7 +194,7 @@ for t=1:sample_size
end
% particles selection
tau_tilde = tau_tilde/sum(tau_tilde);
indx = resample(0, tau_tilde', options_.particle);
indx = kitagawa(tau_tilde');
StateVectors = StateVectors(:,indx);
xparam = fore_xparam(:,indx);
if pruning
......@@ -214,7 +206,7 @@ for t=1:sample_size
for i=1:number_of_particles
info = 12042009;
counter=0;
while info(1) && counter <options_.particle.liu_west_max_resampling_tries
while info(1) && counter <online_opt.liu_west_max_resampling_tries
counter=counter+1;
candidate = xparam(:,i) + chol_sigma_bar*randn(number_of_parameters, 1);
if all(candidate>=bounds.lb) && all(candidate<=bounds.ub)
......@@ -298,12 +290,12 @@ for t=1:sample_size
wtilde(i) = w_stage1(i)*exp(-.5*(const_lik+log(det(ReducedForm.H))+sum(PredictionError.*(ReducedForm.H\PredictionError), 1)));
end
end
if counter==options_.particle.liu_west_max_resampling_tries
fprintf('\nLiu & West particle filter: I haven''t been able to solve the model in %u tries.\n',options_.particle.liu_west_max_resampling_tries)
fprintf('Liu & West particle filter: The last error message was: %s\n',get_error_message(info))
fprintf('Liu & West particle filter: You can try to increase liu_west_max_resampling_tries, but most\n')
fprintf('Liu & West particle filter: likely there is an issue with the model.\n')
error('Liu & West particle filter: unable to solve the model.')
if counter==online_opt.liu_west_max_resampling_tries
fprintf('\nLiu & West online filter: I haven''t been able to solve the model in %u tries.\n',online_opt.liu_west_max_resampling_tries)
fprintf('Liu & West online filter: The last error message was: %s\n',get_error_message(info))
fprintf('Liu & West online filter: You can try to increase liu_west_max_resampling_tries, but most\n')
fprintf('Liu & West online filter: likely there is an issue with the model.\n')
error('Liu & West online filter: unable to solve the model.')
end
end
end
......@@ -313,10 +305,8 @@ for t=1:sample_size
variance_update = false;
end
% final resampling (not advised)
if second_resample
[~, idmode] = max(weights);
mode_xparam(:,t) = xparam(:,idmode);
indx = resample(0, weights,options_.particle);
if online_opt.second_resampling
indx = kitagawa(weights);
StateVectors = StateVectors(:,indx) ;
if pruning
StateVectors_ = StateVectors_(:,indx);
......@@ -333,10 +323,11 @@ for t=1:sample_size
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
[~, idmode] = max(weights);
mode_xparam(:,t) = xparam(:,idmode);
if t==sample_size
param = xparam ;
save(sprintf('%s%sparameters_particles_final.mat', SimulationFolder, filesep()), 'param');
end
else
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)');
......@@ -362,22 +353,21 @@ for t=1:sample_size
end
end
end
if t==sample_size
param = xparam(:,kitagawa(weights));
save(sprintf('%s%sparameters_particles_final.mat', SimulationFolder, filesep()), 'param');
end
end
save(sprintf('%s%sparameters_particles-%u.mat', SimulationFolder, filesep(), t), 'xparam');
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('')
end
pmean = xparam(:,sample_size);
pmode = mode_xparam(:,sample_size);
pstdev = std_xparam(:,sample_size) ;
p025 = lb95_xparam(:,sample_size) ;
p975 = ub95_xparam(:,sample_size) ;
pmedian = median_xparam(:,sample_size) ;
covariance = mat_var_cov;
end
%% Plot parameters trajectory
TeX = options_.TeX;
......@@ -395,7 +385,7 @@ end
for plt = 1:nbplt
hh_fig = dyn_figure(options_.nodisplay,'Name','Parameters Trajectories');
for k=1:length(pmean)
for k=1:number_of_parameters
subplot(nr,nc,k)
[name,texname] = get_the_name(k,TeX,M_,estim_params_,options_.varobs);
% Draw the surface for an interval containing 95% of the particles.
......@@ -426,38 +416,3 @@ for plt = 1:nbplt
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
hh_fig = dyn_figure(options_.nodisplay,'Name','Parameters Densities');
for k=1:length(pmean)
subplot(nr,nc,k)
[name,texname] = get_the_name(k,TeX,M_,estim_params_,options_.varobs);
optimal_bandwidth = mh_optimal_bandwidth(xparam(k,:)',number_of_particles,bandwidth,kernel_function);
[density(:,1),density(:,2)] = kernel_density_estimate(xparam(k,:)', number_of_grid_points, ...
number_of_particles, optimal_bandwidth, kernel_function);
plot(density(:,1), density(:,2));
hold on
if TeX
title(texname,'interpreter','latex')
else
title(name,'interpreter','none')
end
hold off
axis tight
drawnow
end
dyn_saveas(hh_fig,[ M_.fname '_param_density' int2str(plt) ],options_.nodisplay,options_.graph_format);
if TeX && any(strcmp('eps',cellstr(options_.graph_format)))
% TeX eps loader file
fprintf(fidTeX, '\\begin{figure}[H]\n');
fprintf(fidTeX,'\\centering \n');
fprintf(fidTeX,'\\includegraphics[width=%2.2f\\textwidth]{%_param_density%s}\n',min(k/nc,1),M_.fname,int2str(plt));
fprintf(fidTeX,'\\caption{Parameter densities based on the Liu/West particle filter.}');
fprintf(fidTeX,'\\label{Fig:ParameterDensities:%s}\n',int2str(plt));
fprintf(fidTeX,'\\end{figure}\n');
fprintf(fidTeX,' \n');
end
end
......@@ -37,7 +37,7 @@ function posterior_sampler(TargetFun,ProposalFun,xparam1,sampler_options,mh_boun
% Then the comments write here can be used for all the other pairs of
% parallel functions and also for management functions.
% Copyright © 2006-2023 Dynare Team
% Copyright © 2006-2024 Dynare Team
%
% This file is part of Dynare.
%
......@@ -128,10 +128,12 @@ else
% Global variables for parallel routines.
globalVars = struct();
% which files have to be copied to run remotely
NamFileInput(1,:) = {'',[ModelName '.static.m']};
NamFileInput(2,:) = {'',[ModelName '.dynamic.m']};
NamFileInput(1,:) = {'',[ModelName '.sparse.static_resid.m']};
NamFileInput(2,:) = {'',[ModelName '.sparse.static_g1.m']};
NamFileInput(3,:) = {'',[ModelName '.sparse.dynamic_resid.m']};
NamFileInput(4,:) = {'',[ModelName '.sparse.dynamic_g1.m']};
if M_.set_auxiliary_variables
NamFileInput(3,:) = {'',[M_.fname '.set_auxiliary_variables.m']};
NamFileInput(5,:) = {'',[M_.fname '.set_auxiliary_variables.m']};
end
if options_.steadystate_flag
if options_.steadystate_flag == 1
......
......@@ -271,10 +271,12 @@ else
localVars.ifil = ifil;
globalVars = [];
% which files have to be copied to run remotely
NamFileInput(1,:) = {'',[M_.fname '.static.m']};
NamFileInput(2,:) = {'',[M_.fname '.dynamic.m']};
NamFileInput(1,:) = {'',[M_.fname '.sparse.static_resid.m']};
NamFileInput(2,:) = {'',[M_.fname '.sparse.static_g1.m']};
NamFileInput(3,:) = {'',[M_.fname '.sparse.dynamic_resid.m']};
NamFileInput(4,:) = {'',[M_.fname '.sparse.dynamic_g1.m']};
if M_.set_auxiliary_variables
NamFileInput(3,:) = {'',[M_.fname '.set_auxiliary_variables.m']};
NamFileInput(5,:) = {'',[M_.fname '.set_auxiliary_variables.m']};
end
if options_.steadystate_flag
if options_.steadystate_flag == 1
......
......@@ -57,6 +57,11 @@ else
posterior = load(sprintf('%s/hssmc/particles-%u-%u.mat', M_.dname, length(pfiles), length(pfiles)));
NumberOfDraws = size(posterior.particles,2);
NumberOfDrawsPerChain=NumberOfDraws;
elseif isonline(options_)
% Load draws from the posterior distribution
posterior = load(sprintf('%s/online/parameters_particles_final.mat', M_.dname));
NumberOfDraws = size(posterior.param,2);
NumberOfDrawsPerChain=NumberOfDraws;
elseif isdime(options_)
posterior = load(sprintf('%s%s%s%schains.mat', M_.dname, filesep(), 'dime', filesep()));
tune = posterior.tune;
......