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
Select Git revision
Loading items

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
  • chskcau/dynare-doc-fixes
27 results
Select Git revision
Loading items
Show changes
Show whitespace changes
Inline Side-by-side
Showing
with 1168 additions and 571 deletions
matlab/compute_mh_covariance_matrix.m
View file @ 28eba237
......@@ -32,7 +32,7 @@ function [posterior_mean,posterior_covariance,posterior_mode,posterior_kernel_at
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
global M_ options_ estim_params_
global M_ options_ estim_params_ bayestopt_
n = estim_params_.np + ...
......@@ -76,5 +76,6 @@ end
xparam1 = posterior_mode';
hh = inv(posterior_covariance);
fval = posterior_kernel_at_the_mode;
parameter_names = bayestopt_.name;
save([M_.fname '_mh_mode.mat'],'xparam1','hh','fval');
\ No newline at end of file
save([M_.fname '_mh_mode.mat'],'xparam1','hh','fval','parameter_names');
\ No newline at end of file
matlab/cycle_reduction.m
View file @ 28eba237
......@@ -116,7 +116,7 @@ end
%$ t = zeros(3,1);
%$
%$ % Set the dimension of the problem to be solved.
%$ n = 2000;
%$ n = 500;
%$
%$ % Set the equation to be solved
%$ A = eye(n);
......
matlab/det_cond_forecast.m
View file @ 28eba237
function det_cond_forecast(constrained_paths, constrained_vars, options_cond_fcst, constrained_perfect_foresight)
% Computes conditional forecasts for a deterministic model.
function data_set = det_cond_forecast(varargin)
% Computes conditional forecasts using the extended path method.
%
% INPUTS
% o constrained_paths [double] m*p array, where m is the number of constrained endogenous variables and p is the number of constrained periods.
% o constrained_vars [char] m*x array holding the names of the controlled endogenous variables.
% o options_cond_fcst [structure] containing the options. The fields are:
% + replic [integer] scalar, number of monte carlo simulations.
% + parameter_set [char] values of the estimated parameters:
% "posterior_mode",
% "posterior_mean",
% "posterior_median",
% "prior_mode" or
% "prior mean".
% [double] np*1 array, values of the estimated parameters.
% + controlled_varexo [char] m*x array, list of controlled exogenous variables.
% + conf_sig [double] scalar in [0,1], probability mass covered by the confidence bands.
% o constrained_perfect_foresight [double] m*1 array indicating if the endogenous variables path is perfectly foresight (1) or is a surprise (0)
% o plan [structure] A structure describing the different shocks and the endogenous varibales, the date of the shocks and the path of the shock.
% The plan structure is created by the functions init_plan, basic_plan and flip_plan
% o [dataset] [dseries] A dserie containing the initial values of the shocks and the endogenous variables (usually the dseries generated by the smoother).
% o [starting_date] [dates] The first date of the forecast.
%
%
% OUTPUTS
% None.
% dataset [dseries] Returns a dseries containing the forecasted endgenous variables and shocks
%
% SPECIAL REQUIREMENTS
% This routine has to be called after an estimation statement or an estimated_params block.
%
% REMARKS
% [1] Results are stored in a structure which is saved in a mat file called conditional_forecasts.mat.
% [2] Use the function plot_icforecast to plot the results.
% Copyright (C) 2013 Dynare Team
% Copyright (C) 2013-2014 Dynare Team
%
% This file is part of Dynare.
%
......@@ -46,29 +30,251 @@ function det_cond_forecast(constrained_paths, constrained_vars, options_cond_fcs
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
global options_ oo_ M_
pp = 2;
initial_conditions = oo_.steady_state;
%We have to get an initial guess for the conditional forecast
% and terminal conditions for the non-stationary variables, we
% use the first order approximation of the rational expectation solution.
dr = struct();
oo_.dr=set_state_space(dr,M_,options_);
options_.order = 1;
[dr,Info,M_,options_,oo_] = resol(0,M_,options_,oo_);
b_surprise = 0;
b_pf = 0;
surprise = 0;
pf = 0;
is_shock = [];
is_constraint = [];
if length(varargin) > 3
% regular way to call
constrained_paths = varargin{1};
max_periods_simulation = size(constrained_paths, 2);
constrained_vars = varargin{2};
options_cond_fcst = varargin{3};
constrained_perfect_foresight = varargin{4};
constraint_index = cell(max_periods_simulation,1);
nvars = length(constrained_vars);
for i = 1:max_periods_simulation
constraint_index{i} = 1:nvars;
end;
direct_mode = 0;
shocks_present = 0;
controlled_varexo = options_cond_fcst.controlled_varexo;
nvarexo = size(controlled_varexo, 1);
options_cond_fcst.controlled_varexo = zeros(nvarexo,1);
exo_names = cell(M_.exo_nbr,1);
for i = 1:M_.exo_nbr
exo_names{i} = deblank(M_.exo_names(i,:));
end
for i = 1:nvarexo
j = find(strcmp(controlled_varexo(i,:), exo_names));
if ~isempty(j)
options_cond_fcst.controlled_varexo(i) = j;
else
error(['Unknown exogenous variable ' controlled_varexo(i,:)]);
end
end
else
% alternative way to call
plan = varargin{1};
if length(varargin) >= 2
dset = varargin{2};
if ~isa(dset,'dseries')
error('the second argmuent should be a dseries');
end
if length(varargin) >= 3
range = varargin{3};
if ~isa(range,'dates')
error('the third argmuent should be a dates');
end
%if (range(range.ndat) > dset.time(dset.nobs) )
if (range(range.ndat) > dset.dates(dset.nobs) )
s1 = strings(dset.dates(dset.nobs));
s2 = strings(range(range.ndat));
error(['the dseries ' inputname(2) ' finish at time ' s1{1} ' before the last period of forecast ' s2{1}]);
end
sym_dset = dset(dates(-range(1)):dates(range(range.ndat)));
oo_.exo_simul = repmat(oo_.exo_steady_state',max(range.ndat + 1, options_.periods),1);
oo_.endo_simul = repmat(oo_.steady_state, 1, max(range.ndat + 1, options_.periods));
for i = 1:sym_dset.vobs
iy = find(strcmp(strtrim(sym_dset.name{i}), strtrim(plan.endo_names)));
if ~isempty(iy)
oo_.endo_simul(iy,1:sym_dset.nobs) = sym_dset.data(:, i);
initial_conditions(iy) = sym_dset.data(1, i);
else
ix = find(strcmp(strtrim(sym_dset.name{i}), strtrim(plan.exo_names)));
if ~isempty(ix)
oo_.exo_simul(1, ix) = sym_dset.data(1, i)';
else
%warning(['The variable ' sym_dset.name{i} ' in the dataset ' inputname(2) ' is not a endogenous neither an exogenous variable!!']);
end
end
end
for i = 1:length(M_.aux_vars)
if M_.aux_vars(i).type == 1 %lag variable
iy = find(strcmp(deblank(M_.endo_names(M_.aux_vars(i).orig_index,:)), sym_dset.name));
if ~isempty(iy)
oo_.endo_simul(M_.aux_vars(i).endo_index, 1:sym_dset.nobs) = dset(dates(range(1) + (M_.aux_vars(i).orig_lead_lag - 1)):dates(range(range.ndat) + M_.aux_vars(i).orig_lead_lag)).data(:,iy);
initial_conditions(M_.aux_vars(i).endo_index) = dset(dates(range(1) + (M_.aux_vars(i).orig_lead_lag - 1))).data(:,iy);
else
warning(['The variable auxiliary ' M_.endo_names(M_.aux_vars(i).endo_index, :) ' associated to the variable ' M_.endo_names(M_.aux_vars(i).orig_index,:) ' do not appear in the dataset']);
end
else
oo_.endo_simul(M_.aux_vars(i).endo_index, 1:sym_dset.nobs) = repmat(oo_.steady_state(M_.aux_vars(i).endo_index), 1, range.ndat + 1);
end
end
else
error('impossible case');
end;
else
oo_.exo_simul = repmat(oo_.exo_steady_state',options_.periods+2,1);
oo_.endo_simul = repmat(oo_.steady_state, 1, options_.periods+2);
end
direct_mode = 1;
constrained_paths = plan.constrained_paths_;
constrained_vars = plan.constrained_vars_;
options_cond_fcst = plan.options_cond_fcst_;
constrained_perfect_foresight = plan.constrained_perfect_foresight_;
constrained_periods = plan.constrained_date_;
if ~isempty(plan.shock_paths_)
shock_paths = plan.shock_paths_;
shock_vars = plan.shock_vars_;
shock_perfect_foresight = plan.shock_perfect_foresight_;
shock_periods = plan.shock_date_;
shocks_present = 1;
else
shocks_present = 0;
end
total_periods = plan.date;
end;
if ~isfield(options_cond_fcst,'periods') || isempty(options_cond_fcst.periods)
options_cond_fcst.periods = 100;
end
options_.periods = 10;
if direct_mode == 1
n_periods = length(constrained_periods);
is_constraint = zeros(size(total_periods), n_periods);
constrained_paths_cell = constrained_paths;
clear constrained_paths;
constrained_paths = zeros(n_periods, size(total_periods));
max_periods_simulation = 0;
for i = 1:n_periods
period_i = constrained_periods{i};
%period_i
tp = total_periods(1);
if size(period_i) > 1
init_periods = period_i(1);
tp_end = period_i(end);
else
init_periods = period_i;
tp_end = period_i;
end;
tp0 = tp;
while tp < init_periods
tp = tp + 1;
end
j = 0;
while tp <= tp_end
is_constraint(tp - tp0 + 1, i) = 1;
constrained_paths(i, tp - tp0 + 1) = constrained_paths_cell{i}(j + 1);
tp = tp + 1;
j = j + 1;
end;
if tp - tp0 > max_periods_simulation
max_periods_simulation = tp - tp0;
end;
end
n_nnz = length(sum(is_constraint,2));
if n_nnz > 0
constraint_index = cell(n_nnz,1);
for i= 1:n_nnz
constraint_index{i} = find(is_constraint(i,:));
end;
end;
if shocks_present
n_periods = length(shock_periods);
shock_paths_cell = shock_paths;
clear shock_paths;
shock_paths = zeros(n_periods, size(total_periods));
is_shock = zeros(size(total_periods), n_periods);
for i = 1:n_periods
period_i = shock_periods{i};
%period_i
tp = total_periods(1);
if size(period_i) > 1
init_periods = period_i(1);
tp_end = period_i(end);
else
init_periods = period_i;
tp_end = period_i;
end;
tp0 = tp;
while tp < init_periods
tp = tp + 1;
end
j = 0;
while tp <= tp_end
is_shock(tp - tp0 + 1, i) = 1;
shock_paths(i, tp - tp0 + 1) = shock_paths_cell{i}(j + 1);
tp = tp + 1;
j = j + 1;
end;
if tp - tp0 > max_periods_simulation
max_periods_simulation = tp - tp0;
end;
end;
n_nnz = length(sum(is_shock,2));
if n_nnz > 0
shock_index = cell(n_nnz, 1);
for i= 1:n_nnz
shock_index{i} = find(is_shock(i,:));
end;
end
end
else
is_constraint = ones(size(constrained_paths));
end
maximum_lag = M_.maximum_lag;
maximum_lead = M_.maximum_lead;
ys = oo_.steady_state;
ny = size(ys,1);
xs = [oo_.exo_steady_state ; oo_.exo_det_steady_state];
nx = size(xs,1);
constrained_periods = size(constrained_paths,2);
constrained_periods = max_periods_simulation;
n_endo_constrained = size(constrained_vars,1);
if isfield(options_cond_fcst,'controlled_varexo')
n_control_exo = size(options_cond_fcst.controlled_varexo, 1);
if n_control_exo ~= n_endo_constrained
error(['det_cond_forecast:: the number of exogenous controlled variables (' int2str(n_control_exo) ') has to be equal to the number of constrained endogenous variabes (' int2str(n_endo_constrained) ')'])
error(['det_cond_forecast: the number of exogenous controlled variables (' int2str(n_control_exo) ') has to be equal to the number of constrained endogenous variabes (' int2str(n_endo_constrained) ')'])
end;
else
error('det_cond_forecast:: to run a deterministic conditional forecast you have to specified the exogenous variables controlled using the option controlled_varex in forecast command');
error('det_cond_forecast: to run a deterministic conditional forecast you have to specified the exogenous variables controlled using the option controlled_varexo in forecast command');
end;
% if n_endo_constrained == 0
% options_.ep.use_bytecode = options_.bytecode;
% data_set = extended_path(initial_conditions, max_periods_simulation);
% end
if length(varargin) >= 1
controlled_varexo = options_cond_fcst.controlled_varexo;
else
exo_names = M_.exo_names;
controlled_varexo = zeros(1,n_control_exo);
for i = 1:nx
......@@ -78,6 +284,7 @@ for i = 1:nx
end
end
end
end
%todo check if zero => error message
......@@ -99,21 +306,26 @@ else
surprise = length(indx_endo_solve_surprise);
end;
eps = options_.solve_tolf;
maxit = options_.simul.maxit;
initial_conditions = oo_.steady_state;
past_val = 0;
save_options_periods = options_.periods;
options_.periods = options_cond_fcst.periods;
save_options_dynatol_f = options_.dynatol.f;
options_.dynatol.f = 1e-8;
eps1 = 1e-4;
eps1 = 1e-7;%1e-4;
exo = zeros(maximum_lag + options_cond_fcst.periods, nx);
endo = zeros(maximum_lag + options_cond_fcst.periods, ny);
endo(1,:) = oo_.steady_state';
verbosity = options_.verbosity;
options_.verbosity = 0;
% if all the endogenous paths are perfectly anticipated we do not need to
% implement extended path
% implement the extended path
if pf && ~surprise
time_index_constraint = maximum_lag + 1:maximum_lag + constrained_periods;
second_system_size = pf * constrained_periods;
......@@ -129,6 +341,7 @@ if pf && ~surprise
make_y_;
it = 1;
convg = 0;
normra = 1e+50;
while ~convg && it <= maxit
disp('---------------------------------------------------------------------------------------------');
disp(['iteration ' int2str(it)]);
......@@ -137,12 +350,23 @@ if pf && ~surprise
while not_achieved
simul();
result = sum(sum(isfinite(oo_.endo_simul(:,time_index_constraint)))) == ny * constrained_periods;
if (~oo_.deterministic_simulation.status || ~result) && it > 1
if result
y = oo_.endo_simul(constrained_vars, time_index_constraint);
ys = oo_.endo_simul(indx_endo);
col_count = 1;
for j = 1:length(constrained_vars)
y = oo_.endo_simul(constrained_vars(j), time_index_constraint);
r(col_count:col_count+constrained_periods - 1) = (y - constrained_paths(j,1:constrained_periods))';
col_count = col_count + constrained_periods;
end;
normr = norm(r, 1);
end;
if (~oo_.deterministic_simulation.status || ~result || normr > normra) && it > 1
not_achieved = 1;
alpha = alpha / 2;
col_count = 1;
for j = controlled_varexo
oo_.exo_simul(time_index_constraint,j) = (old_exo(:,j) + alpha * D_exo(col_count: col_count + constrained_periods - 1));
for j = controlled_varexo'
oo_.exo_simul(time_index_constraint,j) = (old_exo(:,j) + alpha * D_exo(col_count: (col_count + constrained_periods - 1)));
col_count = col_count + constrained_periods;
end;
disp(['Divergence in Newton: reducing the path length alpha=' num2str(alpha)]);
......@@ -152,26 +376,101 @@ if pf && ~surprise
end;
end;
y = oo_.endo_simul(constrained_vars(j), time_index_constraint);
ys = oo_.endo_simul(indx_endo);
col_count = 1;
for j = 1:length(constrained_vars)
y = oo_.endo_simul(constrained_vars(j), time_index_constraint);
r(col_count:col_count+constrained_periods - 1) = (y - constrained_paths(j,1:constrained_periods))';
col_count = col_count + constrained_periods;
per = 30;
z = oo_.endo_simul(:, 1 : per + 2 );
zx = oo_.exo_simul(1: per + 2,:);
g1 = spalloc(M_.endo_nbr * (per ), M_.endo_nbr * (per ), 3* M_.endo_nbr * per );
g1_x = spalloc(M_.endo_nbr * (per ), M_.exo_nbr, M_.endo_nbr * (per )* M_.exo_nbr );
lag_indx = find(M_.lead_lag_incidence(1,:));
cur_indx = M_.endo_nbr + find(M_.lead_lag_incidence(2,:));
lead_indx = 2 * M_.endo_nbr + find(M_.lead_lag_incidence(3,:));
cum_l1 = 0;
cum_index_d_y_x = [];
indx_x = [];
for k = 1 : per
if k == 1
if (isfield(M_,'block_structure'))
data1 = M_.block_structure.block;
Size = length(M_.block_structure.block);
else
data1 = M_;
Size = 1;
end;
col_count = 1;
for j = controlled_varexo
for time = time_index_constraint
saved = oo_.exo_simul(time,j);
oo_.exo_simul(time,j) = oo_.exo_simul(time,j) + eps1;
simul();
J(:,col_count) = (oo_.endo_simul(indx_endo) - ys) / eps1;
oo_.exo_simul(time,j) = saved;
col_count = col_count + 1;
data1 = M_;
if (options_.bytecode)
[chck, zz, data1]= bytecode('dynamic','evaluate', z, zx, M_.params, oo_.steady_state, k, data1);
else
[zz, g1b] = feval([M_.fname '_dynamic'], z', zx, M_.params, oo_.steady_state, k);
data1.g1_x = g1b(:,end - M_.exo_nbr + 1:end);
data1.g1 = g1b(:,1 : end - M_.exo_nbr);
chck = 0;
end;
mexErrCheck('bytecode', chck);
end;
normr = norm(r, 1);
if k == 1
g1(1:M_.endo_nbr,-M_.endo_nbr + [cur_indx lead_indx]) = data1.g1(:,M_.nspred + 1:end);
elseif k == per
g1(M_.endo_nbr * (k - 1) + 1 :M_.endo_nbr * k,M_.endo_nbr * (k -2) + [lag_indx cur_indx]) = data1.g1(:,1:M_.nspred + M_.endo_nbr);
else
g1(M_.endo_nbr * (k - 1) + 1 :M_.endo_nbr * k, M_.endo_nbr * (k -2) + [lag_indx cur_indx lead_indx]) = data1.g1;
end;
l2 = 1;
pf_c = 1;
if k <= constrained_periods
for l = constraint_index{k}
l1 = controlled_varexo(l);
g1_x(M_.endo_nbr * (k - 1) + 1:M_.endo_nbr * k,1 + cum_l1) = data1.g1_x(:,l1);
if k == 1
indx_x(l2) = l ;
l2 = l2 + 1;
for ii = 2:constrained_periods
indx_x(l2) = length(controlled_varexo) + pf * (ii - 2) + constraint_index{k + ii - 1}(pf_c);
l2 = l2 + 1;
end;
pf_c = pf_c + 1;
cum_index_d_y_x = [cum_index_d_y_x; constrained_vars(l)];
else
cum_index_d_y_x = [cum_index_d_y_x; constrained_vars(l) + (k - 1) * M_.endo_nbr];
end
cum_l1 = cum_l1 + length(l1);
end;
end;
end;
d_y_x = - g1 \ g1_x;
cum_l1 = 0;
count_col = 1;
cum_index_J = 1:length(cum_index_d_y_x(indx_x));
J= zeros(length(cum_index_J));
for j1 = 1:length(controlled_varexo)
cum_l1 = 0;
for k = 1:(constrained_periods)
l1 = constraint_index{k};
l1 = find(constrained_perfect_foresight(l1) | (k == 1));
if constraint_index{k}( j1)
J(cum_index_J,count_col) = d_y_x(cum_index_d_y_x(indx_x),indx_x(count_col));
count_col = count_col + 1;
end
cum_l1 = cum_l1 + length(l1);
end
cum_l1 = cum_l1 + length(constrained_vars(j1));
end;
% col_count = 1;
% for j = controlled_varexo'
% for time = time_index_constraint
% saved = oo_.exo_simul(time,j);
% oo_.exo_simul(time,j) = oo_.exo_simul(time,j) + eps1;
% simul();
% J1(:,col_count) = (oo_.endo_simul(indx_endo) - ys) / eps1;
% oo_.exo_simul(time,j) = saved;
% col_count = col_count + 1;
% end;
% end;
% J1
% sdfmlksdf;
disp(['iteration ' int2str(it) ' error = ' num2str(normr)]);
......@@ -183,36 +482,108 @@ if pf && ~surprise
old_exo = oo_.exo_simul(time_index_constraint,:);
D_exo = - J \ r;
col_count = 1;
for j = controlled_varexo
oo_.exo_simul(time_index_constraint,j) = (oo_.exo_simul(time_index_constraint,j) + D_exo(col_count: col_count + constrained_periods - 1));
col_count = col_count + constrained_periods;
%constrained_periods
for j = controlled_varexo'
oo_.exo_simul(time_index_constraint,j) = oo_.exo_simul(time_index_constraint,j) + D_exo(col_count: (col_count + constrained_periods - 1));
col_count = col_count + constrained_periods - 1;
end;
end;
it = it + 1;
normra = normr;
end;
endo(maximum_lag + 1 :maximum_lag + options_cond_fcst.periods ,:) = oo_.endo_simul(:,maximum_lag + 1:maximum_lag + options_cond_fcst.periods )';
endo = oo_.endo_simul';
exo = oo_.exo_simul;
else
for t = 1:constrained_periods
disp('=============================================================================================');
disp(['t=' int2str(t) ' constrained_paths(:,t)=' num2str(constrained_paths(:,t)') ]);
disp('=============================================================================================');
if direct_mode && ~isempty(is_constraint)
[pos_constrained_pf, junk] = find(constrained_perfect_foresight .* is_constraint(t, :)');
indx_endo_solve_pf = constrained_vars(pos_constrained_pf);
if isempty(indx_endo_solve_pf)
pf = 0;
else
pf = length(indx_endo_solve_pf);
end;
[pos_constrained_surprise, junk] = find((1-constrained_perfect_foresight) .* is_constraint(t, :)');
indx_endo_solve_surprise = constrained_vars(pos_constrained_surprise);
if isempty(indx_endo_solve_surprise)
surprise = 0;
else
surprise = length(indx_endo_solve_surprise);
end;
end;
if direct_mode && ~isempty(is_shock)
[pos_shock_pf, junk] = find(shock_perfect_foresight .* is_shock(t, :)');
indx_endo_solve_pf = shock_vars(pos_shock_pf);
if isempty(indx_endo_solve_pf)
b_pf = 0;
else
b_pf = length(indx_endo_solve_pf);
end;
[pos_shock_surprise, junk] = find((1-shock_perfect_foresight) .* is_shock(t, :)');
indx_endo_solve_surprise = shock_vars(pos_shock_surprise);
if isempty(indx_endo_solve_surprise)
b_surprise = 0;
else
b_surprise = length(indx_endo_solve_surprise);
end;
end;
disp('===============================================================================================');
disp(['t=' int2str(t) ' conditional (surprise=' int2str(surprise) ' perfect foresight=' int2str(pf) ') unconditional (surprise=' int2str(b_surprise) ' perfect foresight=' int2str(b_pf) ')']);
disp('===============================================================================================');
if t == 1
make_ex_;
if maximum_lag > 0
exo_init = oo_.exo_simul;
else
exo_init = zeros(size(oo_.exo_simul));
end
make_y_;
end;
%exo_init
oo_.exo_simul = exo_init;
oo_.endo_simul(:,1) = initial_conditions;
convg = 0;
it = 1;
time_index_constraint = maximum_lag + 1:maximum_lag + constrained_periods - t + 1;
if direct_mode
nb_shocks = length(plan.shock_vars_);
if nb_shocks > 0
shock_index_t = shock_index{t};
shock_vars_t = shock_vars(shock_index_t);
for shock_indx = shock_index_t
if shock_perfect_foresight(shock_indx)
oo_.exo_simul(time_index_constraint,shock_vars_t(shock_indx)) = shock_paths(shock_indx,t:constrained_periods);
else
oo_.exo_simul(maximum_lag + 1,shock_vars_t(shock_indx)) = shock_paths(shock_indx,t);
end
end
end
end
%Compute the initial path using the first order solution around the
% steady-state
oo_.endo_simul(:, 1) = oo_.endo_simul(:, 1) - oo_.steady_state;
for jj = 2 : (options_.periods + 2)
oo_.endo_simul(:,jj) = oo_.dr.ghx(oo_.dr.inv_order_var,:) * oo_.endo_simul(oo_.dr.state_var, jj-1) + oo_.dr.ghu * oo_.exo_simul(jj,:)';
end
for jj = 1 : (options_.periods + 2)
oo_.endo_simul(:,jj) = oo_.steady_state + oo_.endo_simul(:,jj);
end
constraint_index_t = constraint_index{t};
controlled_varexo_t = controlled_varexo(constraint_index_t);
constrained_vars_t = constrained_vars(constraint_index_t);
second_system_size = surprise + pf * (constrained_periods - t + 1);
indx_endo = zeros(second_system_size,1);
col_count = 1;
for j = 1:length(constrained_vars)
for j = 1:length(constrained_vars_t)
if constrained_perfect_foresight(j)
indx_endo(col_count : col_count + constrained_periods - t) = constrained_vars(j) + (time_index_constraint - 1) * ny;
col_count = col_count + constrained_periods - t + 1;
......@@ -222,12 +593,12 @@ else
end;
end;
J = zeros(second_system_size,second_system_size);
r = zeros(second_system_size,1);
convg = 0;
it = 1;
while ~convg && it <= maxit
disp('---------------------------------------------------------------------------------------------');
disp('-----------------------------------------------------------------------------------------------');
disp(['iteration ' int2str(it) ' time ' int2str(t)]);
not_achieved = 1;
alpha = 1;
......@@ -238,8 +609,9 @@ else
not_achieved = 1;
alpha = alpha / 2;
col_count = 1;
for j = controlled_varexo
if constrained_perfect_foresight(j)
for j1 = constraint_index_t
j = controlled_varexo(j1);
if constrained_perfect_foresight(j1)
oo_.exo_simul(time_index_constraint,j) = (old_exo(time_index_constraint,j) + alpha * D_exo(col_count: col_count + constrained_periods - t));
col_count = col_count + constrained_periods - t + 1;
else
......@@ -260,7 +632,7 @@ else
ys = oo_.endo_simul(indx_endo);
col_count = 1;
for j = 1:length(constrained_vars)
for j = constraint_index_t
if constrained_perfect_foresight(j)
y = oo_.endo_simul(constrained_vars(j), time_index_constraint);
r(col_count:col_count+constrained_periods - t) = (y - constrained_paths(j,t:constrained_periods))';
......@@ -273,27 +645,133 @@ else
end
disp('computation of derivatives w.r. to exogenous shocks');
col_count = 1;
for j = 1:length(controlled_varexo)
j_pos = controlled_varexo(j);
if constrained_perfect_foresight(j)
for time = time_index_constraint
saved = oo_.exo_simul(time,j_pos);
oo_.exo_simul(time,j_pos) = oo_.exo_simul(time,j_pos) + eps1;
simul();
J(:,col_count) = (oo_.endo_simul(indx_endo) - ys) / eps1;
oo_.exo_simul(time,j_pos) = saved;
col_count = col_count + 1;
per = 30;
z = oo_.endo_simul(:, 1 : per + 2 );
zx = oo_.exo_simul(1: per + 2,:);
g1 = spalloc(M_.endo_nbr * (per ), M_.endo_nbr * (per ), 3* M_.endo_nbr * per );
g1_x = spalloc(M_.endo_nbr * (per ), M_.exo_nbr, M_.endo_nbr * (per )* M_.exo_nbr );
lag_indx = find(M_.lead_lag_incidence(1,:));
cur_indx = M_.endo_nbr + find(M_.lead_lag_incidence(2,:));
lead_indx = 2 * M_.endo_nbr + find(M_.lead_lag_incidence(3,:));
cum_l1 = 0;
%indx_x = zeros(length(constraint_index_t)* constrained_periods, 1);
cum_index_d_y_x = [];
indx_x = [];
for k = 1 : per
if k == 1
if (isfield(M_,'block_structure'))
data1 = M_.block_structure.block;
Size = length(M_.block_structure.block);
else
data1 = M_;
Size = 1;
end;
data1 = M_;
if (options_.bytecode)
[chck, zz, data1]= bytecode('dynamic','evaluate', z, zx, M_.params, oo_.steady_state, k, data1);
else
saved = oo_.exo_simul(maximum_lag+1,j_pos);
oo_.exo_simul(maximum_lag+1,j_pos) = oo_.exo_simul(maximum_lag+1,j_pos) + eps1;
simul();
J(:,col_count) = (oo_.endo_simul(indx_endo) - ys) / eps1;
oo_.exo_simul(maximum_lag+1,j_pos) = saved;
col_count = col_count + 1;
[zz, g1b] = feval([M_.fname '_dynamic'], z', zx, M_.params, oo_.steady_state, k);
data1.g1_x = g1b(:,end - M_.exo_nbr + 1:end);
data1.g1 = g1b(:,1 : end - M_.exo_nbr);
chck = 0;
end;
mexErrCheck('bytecode', chck);
end;
if k == 1
g1(1:M_.endo_nbr,-M_.endo_nbr + [cur_indx lead_indx]) = data1.g1(:,M_.nspred + 1:end);
elseif k == per
g1(M_.endo_nbr * (k - 1) + 1 :M_.endo_nbr * k,M_.endo_nbr * (k -2) + [lag_indx cur_indx]) = data1.g1(:,1:M_.nspred + M_.endo_nbr);
else
g1(M_.endo_nbr * (k - 1) + 1 :M_.endo_nbr * k, M_.endo_nbr * (k -2) + [lag_indx cur_indx lead_indx]) = data1.g1;
end;
l2 = 1;
pf_c = 1;
if t + k - 1 <= constrained_periods
for l = constraint_index{t + k - 1}
l1 = controlled_varexo(l);
if (k == 1) || ((k > 1) && constrained_perfect_foresight(l))
g1_x(M_.endo_nbr * (k - 1) + 1:M_.endo_nbr * k,1 + cum_l1) = data1.g1_x(:,l1);
if k == 1
if constrained_perfect_foresight(l)
indx_x(l2) = l ;
l2 = l2 + 1;
for ii = 2:constrained_periods - t + 1
indx_x(l2) = length(controlled_varexo) + pf * (ii - 2) + constraint_index{k + ii - 1}(pf_c);
l2 = l2 + 1;
end;
pf_c = pf_c + 1;
else
indx_x(l2) = l;
l2 = l2 + 1;
end;
cum_index_d_y_x = [cum_index_d_y_x; constrained_vars_t(l)];
else
cum_index_d_y_x = [cum_index_d_y_x; constrained_vars_t(l) + (k - 1) * M_.endo_nbr];
end
cum_l1 = cum_l1 + length(l1);
end;
end;
end;
end;
d_y_x = - g1 \ g1_x;
cum_l1 = 0;
count_col = 1;
cum_index_J = 1:length(cum_index_d_y_x(indx_x));
J= zeros(length(cum_index_J));
for j1 = constraint_index_t
if constrained_perfect_foresight(j1)
cum_l1 = 0;
for k = 1:(constrained_periods - t + 1)
l1 = constraint_index{k};
l1 = find(constrained_perfect_foresight(l1) | (k == 1));
if constraint_index{k}( j1)
J(cum_index_J,count_col) = d_y_x(cum_index_d_y_x(indx_x),indx_x(count_col));
count_col = count_col + 1;
end
cum_l1 = cum_l1 + length(l1);
end
else
J(cum_index_J,count_col) = d_y_x(cum_index_d_y_x(indx_x),indx_x(count_col));
count_col = count_col + 1;
end
cum_l1 = cum_l1 + length(constrained_vars_t(j1));
end;
% % Numerical computation of the derivatives in the second systme
% J1 = [];
% col_count = 1;
% for j = constraint_index_t
% j_pos = controlled_varexo(j);
% if constrained_perfect_foresight(j)
% for time = time_index_constraint
% saved = oo_.exo_simul(time,j_pos);
% oo_.exo_simul(time,j_pos) = oo_.exo_simul(time,j_pos) + eps1;
% simul();
% J1(:,col_count) = (oo_.endo_simul(indx_endo) - ys) / eps1;
% oo_.exo_simul(time,j_pos) = saved;
% col_count = col_count + 1;
% end;
% else
% saved = oo_.exo_simul(maximum_lag+1,j_pos);
% oo_.exo_simul(maximum_lag+1,j_pos) = oo_.exo_simul(maximum_lag+1,j_pos) + eps1;
% simul();
% % indx_endo
% J1(:,col_count) = (oo_.endo_simul(indx_endo) - ys) / eps1;
% % J(:,col_count) = (oo_.endo_simul((pp - 1) * M_.endo_nbr + 1: pp * M_.endo_nbr) - ys) / eps1;
% oo_.exo_simul(maximum_lag+1,j_pos) = saved;
% col_count = col_count + 1;
% end;
% end;
% disp('J1');
% disp(full(J1));
% sdfmlk;
normr = norm(r, 1);
......@@ -304,13 +782,22 @@ else
disp('convergence achieved');
else
% Newton update on exogenous shocks
try
D_exo = - J \ r;
catch
[V, D] = eig(full(J));
ev = diag(D);
[ev abs(ev)]
z_root = find(abs(ev) < 1e-4);
z_root
disp(V(:,z_root));
end;
old_exo = oo_.exo_simul;
col_count = 1;
for j = 1:length(controlled_varexo)
for j = constraint_index_t
j_pos=controlled_varexo(j);
if constrained_perfect_foresight(j)
oo_.exo_simul(time_index_constraint,j_pos) = (oo_.exo_simul(time_index_constraint,j_pos) + D_exo(col_count: col_count + constrained_periods - t));
oo_.exo_simul(time_index_constraint,j_pos) = (oo_.exo_simul(time_index_constraint,j_pos) + D_exo(col_count:(col_count + constrained_periods - t) ));
col_count = col_count + constrained_periods - t + 1;
else
oo_.exo_simul(maximum_lag + 1,j_pos) = oo_.exo_simul(maximum_lag + 1,j_pos) + D_exo(col_count);
......@@ -323,11 +810,11 @@ else
if ~convg
error(['convergence not achived at time ' int2str(t) ' after ' int2str(it) ' iterations']);
end;
for j = 1:length(controlled_varexo)
for j = constraint_index_t
j_pos=controlled_varexo(j);
if constrained_perfect_foresight(j)
% in case of mixed surprise and perfect foresight
% endogenous path at each date all the exogenous paths have to be
% in case of mixed surprise and perfect foresight on the
% endogenous path, at each date all the exogenous paths have to be
% stored. The paths are stacked in exo.
for time = time_index_constraint;
exo(past_val + time,j_pos) = oo_.exo_simul(time,j_pos);
......@@ -343,15 +830,17 @@ else
endo(maximum_lag + t :maximum_lag + options_cond_fcst.periods ,:) = ycc(:,maximum_lag + 1:maximum_lag + options_cond_fcst.periods - constrained_periods + 1)';
end;
initial_conditions = yc;
if maximum_lag > 0
exo_init(1,:) = exo(maximum_lag + t,:);
end;
end;
end;
options_.periods = save_options_periods;
options_.dynatol.f = save_options_dynatol_f;
options_.initval_file = save_options_initval_file;
% disp('endo');
% disp(endo);
% disp('exo');
% disp(exo);
options_.verbosity = verbosity;
oo_.endo_simul = endo';
oo_.exo_simul = exo;
if direct_mode
data_set = dseries([endo(2:end,1:M_.orig_endo_nbr) exo(2:end,:)], total_periods(1), {plan.endo_names{:} plan.exo_names{:}}, {plan.endo_names{:} plan.exo_names{:}});
end;
matlab/disp_moments.m
View file @ 28eba237
......@@ -80,6 +80,10 @@ if options_.nocorr == 0
end
end
if options_.noprint == 0 && length(options_.conditional_variance_decomposition)
fprintf('\nSTOCH_SIMUL: conditional_variance_decomposition requires theoretical moments, i.e. periods=0.\n')
end
ar = options_.ar;
if ar > 0
autocorr = [];
......
matlab/disp_th_moments.m
View file @ 28eba237
......@@ -52,7 +52,7 @@ oo_.var = oo_.gamma_y{1};
if ~options_.noprint %options_.nomoments == 0
if options_.order == 2
title='APROXIMATED THEORETICAL MOMENTS';
title='APPROXIMATED THEORETICAL MOMENTS';
else
title='THEORETICAL MOMENTS';
end
......
matlab/display_estimation_results_table.m
View file @ 28eba237
......@@ -183,13 +183,12 @@ if any(bayestopt_.pshape > 0) && options_.TeX %% Bayesian estimation (posterior
fprintf(fidTeX,' \n');
fprintf(fidTeX,'\\begin{center}\n');
fprintf(fidTeX,'\\begin{longtable}{l|lcccc} \n');
fprintf(fidTeX,'\\caption{Results from posterior maximization (parameters)}\n ');
fprintf(fidTeX,'\\caption{Results from posterior maximization (parameters)}\\\\\n ');
fprintf(fidTeX,'\\label{Table:Posterior:1}\\\\\n');
fprintf(fidTeX,'\\hline\\hline \\\\ \n');
fprintf(fidTeX,' & Prior distribution & Prior mean & Prior s.d. & Posterior mode & s.d. \\\\ \n');
fprintf(fidTeX,'\\hline \\endfirsthead \n');
fprintf(fidTeX,'\\caption{(continued)}\n ');
fprintf(fidTeX,'\\label{Table:Posterior:1}\\\\\n');
fprintf(fidTeX,'\\caption{(continued)}\\\\\n ');
fprintf(fidTeX,'\\hline\\hline \\\\ \n');
fprintf(fidTeX,' & Prior distribution & Prior mean & Prior s.d. & Posterior mode & s.d. \\\\ \n');
fprintf(fidTeX,'\\hline \\endhead \n');
......@@ -221,12 +220,12 @@ if any(bayestopt_.pshape > 0) && options_.TeX %% Bayesian estimation (posterior
fprintf(fidTeX,' \n');
fprintf(fidTeX,'\\begin{center}\n');
fprintf(fidTeX,'\\begin{longtable}{l|lcccc} \n');
fprintf(fidTeX,'\\caption{Results from posterior maximization (standard deviation of structural shocks)}\n ');
fprintf(fidTeX,'\\caption{Results from posterior maximization (standard deviation of structural shocks)}\\\\\n ');
fprintf(fidTeX,'\\label{Table:Posterior:2}\\\\\n');
fprintf(fidTeX,'\\hline\\hline \\\\ \n');
fprintf(fidTeX,' & Prior distribution & Prior mean & Prior s.d. & Posterior mode & s.d. \\\\ \n');
fprintf(fidTeX,'\\hline \\endfirsthead \n');
fprintf(fidTeX,'\\caption{(continued)}\n ');
fprintf(fidTeX,'\\caption{(continued)}\\\\\n ');
fprintf(fidTeX,'\\label{Table:Posterior:2}\\\\\n');
fprintf(fidTeX,'\\hline\\hline \\\\ \n');
fprintf(fidTeX,' & Prior distribution & Prior mean & Prior s.d. & Posterior mode & s.d. \\\\ \n');
......@@ -260,12 +259,12 @@ if any(bayestopt_.pshape > 0) && options_.TeX %% Bayesian estimation (posterior
fprintf(fidTeX,' \n');
fprintf(fidTeX,'\\begin{center}\n');
fprintf(fidTeX,'\\begin{longtable}{l|lcccc} \n');
fprintf(fidTeX,'\\caption{Results from posterior maximization (standard deviation of measurement errors)}\n ');
fprintf(fidTeX,'\\caption{Results from posterior maximization (standard deviation of measurement errors)}\\\\\n ');
fprintf(fidTeX,'\\label{Table:Posterior:3}\\\\\n');
fprintf(fidTeX,'\\hline\\hline \\\\ \n');
fprintf(fidTeX,' & Prior distribution & Prior mean & Prior s.d. & Posterior mode & s.d. \\\\ \n');
fprintf(fidTeX,'\\hline \\endfirsthead \n');
fprintf(fidTeX,'\\caption{(continued)}\n ');
fprintf(fidTeX,'\\caption{(continued)}\\\\\n ');
fprintf(fidTeX,'\\label{Table:Posterior:3}\\\\\n');
fprintf(fidTeX,'\\hline\\hline \\\\ \n');
fprintf(fidTeX,' & Prior distribution & Prior mean & Prior s.d. & Posterior mode & s.d. \\\\ \n');
......@@ -299,12 +298,12 @@ if any(bayestopt_.pshape > 0) && options_.TeX %% Bayesian estimation (posterior
fprintf(fidTeX,' \n');
fprintf(fidTeX,'\\begin{center}\n');
fprintf(fidTeX,'\\begin{longtable}{l|lcccc} \n');
fprintf(fidTeX,'\\caption{Results from posterior parameters (correlation of structural shocks)}\n ');
fprintf(fidTeX,'\\caption{Results from posterior parameters (correlation of structural shocks)}\\\\\n ');
fprintf(fidTeX,'\\label{Table:Posterior:4}\\\\\n');
fprintf(fidTeX,'\\hline\\hline \\\\ \n');
fprintf(fidTeX,' & Prior distribution & Prior mean & Prior s.d. & Posterior mode & s.d. \\\\ \n');
fprintf(fidTeX,'\\hline \\endfirsthead \n');
fprintf(fidTeX,'\\caption{(continued)}\n ');
fprintf(fidTeX,'\\caption{(continued)}\\\\\n ');
fprintf(fidTeX,'\\label{Table:Posterior:4}\\\\\n');
fprintf(fidTeX,'\\hline\\hline \\\\ \n');
fprintf(fidTeX,' & Prior distribution & Prior mean & Prior s.d. & Posterior mode & s.d. \\\\ \n');
......@@ -339,12 +338,12 @@ if any(bayestopt_.pshape > 0) && options_.TeX %% Bayesian estimation (posterior
fprintf(fidTeX,' \n');
fprintf(fidTeX,'\\begin{center}\n');
fprintf(fidTeX,'\\begin{longtable}{l|lcccc} \n');
fprintf(fidTeX,'\\caption{Results from posterior parameters (correlation of measurement errors)}\n ');
fprintf(fidTeX,'\\caption{Results from posterior parameters (correlation of measurement errors)}\\\\\n ');
fprintf(fidTeX,'\\label{Table:Posterior:5}\\\\\n');
fprintf(fidTeX,'\\hline\\hline \\\\ \n');
fprintf(fidTeX,' & Prior distribution & Prior mean & Prior s.d. & Posterior mode & s.d. \\\\ \n');
fprintf(fidTeX,'\\hline \\endfirsthead \n');
fprintf(fidTeX,'\\caption{(continued)}\n ');
fprintf(fidTeX,'\\caption{(continued)}\\\\\n ');
fprintf(fidTeX,'\\label{Table:Posterior:5}\\\\\n');
fprintf(fidTeX,'\\hline\\hline \\\\ \n');
fprintf(fidTeX,' & Prior distribution & Prior mean & Prior s.d. & Posterior mode & s.d. \\\\ \n');
......@@ -380,12 +379,12 @@ elseif all(bayestopt_.pshape == 0) && options_.TeX %% MLE and GMM Latex output
fprintf(fidTeX,' \n');
fprintf(fidTeX,'\\begin{center}\n');
fprintf(fidTeX,'\\begin{longtable}{l|lcc} \n');
fprintf(fidTeX,['\\caption{Results from ' table_title ' maximization (parameters)}\n ']);
fprintf(fidTeX,['\\caption{Results from ' table_title ' maximization (parameters)}\\\\\n ']);
fprintf(fidTeX,['\\label{Table:' LaTeXtitle ':1}\\\\\n']);
fprintf(fidTeX,'\\hline\\hline \\\\ \n');
fprintf(fidTeX,' & Mode & s.d. & t-stat\\\\ \n');
fprintf(fidTeX,'\\hline \\endfirsthead \n');
fprintf(fidTeX,'\\caption{(continued)}\n ');
fprintf(fidTeX,'\\caption{(continued)}\\\\\n ');
fprintf(fidTeX,['\\label{Table:' LaTeXtitle ':1}\\\\\n']);
fprintf(fidTeX,'\\hline\\hline \\\\ \n');
fprintf(fidTeX,' & Mode & s.d. & t-stat\\\\ \n');
......@@ -416,12 +415,12 @@ elseif all(bayestopt_.pshape == 0) && options_.TeX %% MLE and GMM Latex output
fprintf(fidTeX,' \n');
fprintf(fidTeX,'\\begin{center}\n');
fprintf(fidTeX,'\\begin{longtable}{l|lcc} \n');
fprintf(fidTeX,['\\caption{Results from ' table_title ' maximization (standard deviation of structural shocks)}\n ']);
fprintf(fidTeX,['\\caption{Results from ' table_title ' maximization (standard deviation of structural shocks)}\\\\\n ']);
fprintf(fidTeX,['\\label{Table:' LaTeXtitle ':2}\\\\\n']);
fprintf(fidTeX,'\\hline\\hline \\\\ \n');
fprintf(fidTeX,' & Mode & s.d. & t-stat\\\\ \n');
fprintf(fidTeX,'\\hline \\endfirsthead \n');
fprintf(fidTeX,'\\caption{(continued)}\n ');
fprintf(fidTeX,'\\caption{(continued)}\\\\\n ');
fprintf(fidTeX,['\\label{Table:' LaTeXtitle ':2}\\\\\n']);
fprintf(fidTeX,'\\hline\\hline \\\\ \n');
fprintf(fidTeX,' & Mode & s.d. & t-stat\\\\ \n');
......@@ -453,12 +452,12 @@ elseif all(bayestopt_.pshape == 0) && options_.TeX %% MLE and GMM Latex output
fprintf(fidTeX,' \n');
fprintf(fidTeX,'\\begin{center}\n');
fprintf(fidTeX,'\\begin{longtable}{l|lcc} \n');
fprintf(fidTeX,['\\caption{Results from ' table_title ' maximization (standard deviation of measurement errors)}\n ']);
fprintf(fidTeX,['\\caption{Results from ' table_title ' maximization (standard deviation of measurement errors)}\\\\\n ']);
fprintf(fidTeX,['\\label{Table:' LaTeXtitle ':3}\\\\\n']);
fprintf(fidTeX,'\\hline\\hline \\\\ \n');
fprintf(fidTeX,' & Mode & s.d. & t-stat\\\\ \n');
fprintf(fidTeX,'\\hline \\endfirsthead \n');
fprintf(fidTeX,'\\caption{(continued)}\n ');
fprintf(fidTeX,'\\caption{(continued)}\\\\\n ');
fprintf(fidTeX,['\\label{Table:' LaTeXtitle ':3}\\\\\n']);
fprintf(fidTeX,'\\hline\\hline \\\\ \n');
fprintf(fidTeX,' & Mode & s.d. & t-stat\\\\ \n');
......@@ -490,12 +489,12 @@ elseif all(bayestopt_.pshape == 0) && options_.TeX %% MLE and GMM Latex output
fprintf(fidTeX,' \n');
fprintf(fidTeX,'\\begin{center}\n');
fprintf(fidTeX,'\\begin{longtable}{l|lcc} \n');
fprintf(fidTeX,['\\caption{Results from ' table_title ' maximization (correlation of structural shocks)}\n ']);
fprintf(fidTeX,['\\caption{Results from ' table_title ' maximization (correlation of structural shocks)}\\\\\n ']);
fprintf(fidTeX,['\\label{Table:' LaTeXtitle ':4}\\\\\n']);
fprintf(fidTeX,'\\hline\\hline \\\\ \n');
fprintf(fidTeX,' & Mode & s.d. & t-stat\\\\ \n');
fprintf(fidTeX,'\\hline \\endfirsthead \n');
fprintf(fidTeX,'\\caption{(continued)}\n ');
fprintf(fidTeX,'\\caption{(continued)}\\\\\n ');
fprintf(fidTeX,['\\label{Table:' LaTeXtitle ':4}\\\\\n']);
fprintf(fidTeX,'\\hline\\hline \\\\ \n');
fprintf(fidTeX,' & Mode & s.d. & t-stat\\\\ \n');
......@@ -528,12 +527,12 @@ elseif all(bayestopt_.pshape == 0) && options_.TeX %% MLE and GMM Latex output
fprintf(fidTeX,' \n');
fprintf(fidTeX,'\\begin{center}\n');
fprintf(fidTeX,'\\begin{longtable}{l|lcc} \n');
fprintf(fidTeX,['\\caption{Results from ' table_title ' maximization (correlation of measurement errors)}\n ']);
fprintf(fidTeX,['\\caption{Results from ' table_title ' maximization (correlation of measurement errors)}\\\\\n ']);
fprintf(fidTeX,['\\label{Table:' LaTeXtitle ':5}\\\\\n']);
fprintf(fidTeX,'\\hline\\hline \\\\ \n');
fprintf(fidTeX,' & Mode & s.d. & t-stat\\\\ \n');
fprintf(fidTeX,'\\hline \\endfirsthead \n');
fprintf(fidTeX,'\\caption{(continued)}\n ');
fprintf(fidTeX,'\\caption{(continued)}\\\\\n ');
fprintf(fidTeX,['\\label{Table:' LaTeXtitle ':5}\\\\\n']);
fprintf(fidTeX,'\\hline\\hline \\\\ \n');
fprintf(fidTeX,' & Mode & s.d. & t-stat\\\\ \n');
......
matlab/dsge_likelihood.m
View file @ 28eba237
......@@ -257,7 +257,7 @@ end
% Return, with endogenous penalty when possible, if dynare_resolve issues an error code (defined in resol).
if info(1) == 1 || info(1) == 2 || info(1) == 5 || info(1) == 7 || info(1) ...
== 8 || info(1) == 22 || info(1) == 24 || info(1) == 19
== 8 || info(1) == 22 || info(1) == 24 || info(1) == 19 || info(1) == 25
fval = objective_function_penalty_base+1;
info = info(1);
exit_flag = 0;
......@@ -627,7 +627,7 @@ if analytic_derivation,
analytic_deriv_info={analytic_derivation,DT,DYss,DOm,DH,DP,asy_Hess};
else
analytic_deriv_info={analytic_derivation,DT,DYss,DOm,DH,DP,D2T,D2Yss,D2Om,D2H,D2P};
clear DT DYss DOm DH DP D2T D2Yss D2Om D2H D2P,
clear DT DYss DOm DP D2T D2Yss D2Om D2H D2P,
end
else
analytic_deriv_info={0};
......
matlab/dsge_var_likelihood.m
View file @ 28eba237
......@@ -114,6 +114,8 @@ if dsge_prior_weight<(NumberOfParameters+NumberOfObservedVariables)/DynareDatase
fval = objective_function_penalty_base+abs(DynareDataset.info.ntobs*dsge_prior_weight-(NumberOfParameters+NumberOfObservedVariables));
exit_flag = 0;
info = 51;
info(2)=dsge_prior_weight;
info(3)=(NumberOfParameters+NumberOfObservedVariables)/DynareDataset.info.ntobs;
return
end
......@@ -127,7 +129,7 @@ end
% Return, with endogenous penalty when possible, if dynare_resolve issues an error code (defined in resol).
if info(1) == 1 || info(1) == 2 || info(1) == 5 || info(1) == 7 || info(1) ...
== 8 || info(1) == 22 || info(1) == 24
== 8 || info(1) == 22 || info(1) == 24 || info(1) == 25
fval = objective_function_penalty_base+1;
info = info(1);
exit_flag = 0;
......
matlab/dyn_first_order_solver.m
View file @ 28eba237
......@@ -261,7 +261,8 @@ else
Z11 = Z(indx_stable_root, indx_stable_root);
Z21 = Z(indx_explosive_root, indx_stable_root);
Z22 = Z(indx_explosive_root, indx_explosive_root);
[minus_gx,rc] = linsolve(Z22,Z21);
opts.TRANSA = false; % needed by Octave 4.0.0
[minus_gx,rc] = linsolve(Z22,Z21,opts);
if rc < 1e-9
% Z22 is near singular
info(1) = 5;
......@@ -273,7 +274,9 @@ else
opts.UT = true;
opts.TRANSA = true;
hx1 = linsolve(tt(indx_stable_root, indx_stable_root),Z11,opts)';
hx2 = linsolve(Z11,ss(indx_stable_root, indx_stable_root)')';
opts.UT = false; % needed by Octave 4.0.0
opts.TRANSA = false; % needed by Octave 4.0.0
hx2 = linsolve(Z11,ss(indx_stable_root, indx_stable_root)',opts)';
hx = hx1*hx2;
ghx = [hx(k1,:); gx(k2(nboth+1:end),:)];
end
......
matlab/dyn_ramsey_static.m
View file @ 28eba237
......@@ -36,9 +36,14 @@ steady_state = [];
params = M.params;
check = 0;
options_.steadystate.nocheck = 1;
% dyn_ramsey_static_1 is a subfunction
nl_func = @(x) dyn_ramsey_static_1(x,M,options_,oo);
if options_.steadystate_flag
% check_static_model is a subfunction
if check_static_model(oo.steady_state,M,options_,oo)
steady_state = oo.steady_state;
return
elseif options_.steadystate_flag
k_inst = [];
instruments = options_.instruments;
inst_nbr = size(options_.instruments,1);
......@@ -46,11 +51,10 @@ if options_.steadystate_flag
k_inst = [k_inst; strmatch(options_.instruments(i,:), ...
M.endo_names,'exact')];
end
ys = oo.steady_state;
if inst_nbr == 1
inst_val = csolve(nl_func,oo.steady_state(k_inst),'',options_.solve_tolf,100);
else
[inst_val,info1] = dynare_solve(nl_func,ys(k_inst),0);
[inst_val,info1] = dynare_solve(nl_func,oo.steady_state(k_inst),0);
end
ys(k_inst) = inst_val;
exo_ss = [oo.exo_steady_state oo.exo_det_steady_state];
......@@ -133,11 +137,11 @@ Uyy = reshape(Uyy,endo_nbr,endo_nbr);
% depends on multipliers to 0 to compute residuals
if (options_.bytecode)
[chck, res, junk] = bytecode('static',xx,[oo.exo_simul oo.exo_det_simul], ...
M.params, 'evaluate');
params, 'evaluate');
fJ = junk.g1;
else
[res,fJ] = feval([fname '_static'],xx,[oo.exo_simul oo.exo_det_simul], ...
M.params);
params);
end
% index of multipliers and corresponding equations
% the auxiliary variables before the Lagrange multipliers are treated
......@@ -145,9 +149,9 @@ end
aux_eq = [1:orig_endo_aux_nbr, orig_endo_aux_nbr+orig_eq_nbr+1:size(fJ,1)];
A = fJ(aux_eq,orig_endo_aux_nbr+1:end);
y = res(aux_eq);
aux_vars = -A\y;
mult = -A\y;
resids1 = y+A*aux_vars;
resids1 = y+A*mult;
if inst_nbr == 1
r1 = sqrt(resids1'*resids1);
else
......@@ -161,4 +165,22 @@ else
resids = [res(orig_endo_nbr+(1:orig_endo_nbr-inst_nbr)); r1];
end
rJ = [];
steady_state = [xx(1:orig_endo_aux_nbr); aux_vars];
\ No newline at end of file
if needs_set_auxiliary_variables
steady_state = s_a_v_func([xx(1:orig_endo_aux_nbr); mult]);
else
steady_state = [xx(1:orig_endo_aux_nbr); mult];
end
function result = check_static_model(ys,M,options_,oo)
result = false;
if (options_.bytecode)
[chck, res, junk] = bytecode('static',ys,[oo.exo_simul oo.exo_det_simul], ...
M.params, 'evaluate');
else
res = feval([M.fname '_static'],ys,[oo.exo_simul oo.exo_det_simul], ...
M.params);
end
if norm(res) < options_.solve_tolf
result = true;
end
matlab/dynare.m
View file @ 28eba237
......@@ -64,7 +64,11 @@ more off
% sets default format for save() command
if isoctave
if octave_ver_less_than('3.8')
default_save_options('-mat')
else
save_default_options('-mat')
end
end
if nargin < 1
......@@ -78,6 +82,7 @@ end
% Testing if file have extension
% If no extension default .mod is added
if isempty(strfind(fname,'.'))
fnamelength = length(fname);
fname1 = [fname '.dyn'];
d = dir(fname1);
if length(d) == 0
......@@ -90,8 +95,21 @@ else
&& ~strcmp(upper(fname(size(fname,2)-3:size(fname,2))),'.DYN')
error('DYNARE: argument must be a filename with .mod or .dyn extension')
end;
fnamelength = length(fname) - 4;
end;
if fnamelength + length('_set_auxiliary_variables') > namelengthmax()
error('The name of your MOD file is too long, please shorten it')
end
% Workaround for a strange bug with Octave: if there is any call to exist(fname)
% before the call to the preprocessor, then Octave will use the old copy of
% the .m instead of the newly generated one. Deleting the .m beforehand
% fixes the problem.
if isoctave && length(dir([fname(1:(end-4)) '.m'])) > 0
delete([fname(1:(end-4)) '.m'])
end
if ~isempty(strfind(fname,filesep))
fprintf('\nIt seems you are trying to call a mod-file not located in the "Current Folder". This is not possible (the %s symbol is not allowed in the name of the mod file).\n', filesep)
[pathtomodfile,basename,ext] = fileparts(fname);
......
matlab/dynare_config.m
View file @ 28eba237
......@@ -60,6 +60,7 @@ addpath([dynareroot '/ep/'])
addpath([dynareroot '/utilities/doc/'])
addpath([dynareroot '/utilities/tests/'])
addpath([dynareroot '/utilities/dates/'])
addpath([dynareroot '/utilities/dseries/'])
addpath([dynareroot '/utilities/dataset/'])
addpath([dynareroot '/utilities/general/'])
addpath([dynareroot '/reports/'])
......@@ -89,13 +90,15 @@ if ~isoctave && matlab_ver_less_than('7.4')
addpath([dynareroot '/missing/bsxfun'])
end
% ilu is missing in old versions of MATLAB and in Octave
if isoctave || matlab_ver_less_than('7.4')
% ilu is missing in old versions of MATLAB and in Octave < 4.0
if (isoctave && octave_ver_less_than('4.0')) || ...
(~isoctave && matlab_ver_less_than('7.4'))
addpath([dynareroot '/missing/ilu'])
end
% strjoin is missing in older versions of MATLAB and in Octave
if isoctave || matlab_ver_less_than('8.1')
% strjoin is missing in older versions of MATLAB and in Octave < 3.8
if (isoctave && octave_ver_less_than('3.8')) || ...
(~isoctave && matlab_ver_less_than('8.1'))
addpath([dynareroot '/missing/strjoin'])
end
......@@ -118,7 +121,7 @@ else
addpath(mexpath)
end
else
mexpath = [dynareroot '../mex/matlab/win32-7.5-8.2'];
mexpath = [dynareroot '../mex/matlab/win32-7.5-8.5'];
if exist(mexpath, 'dir')
addpath(mexpath)
end
......@@ -138,7 +141,7 @@ else
addpath(mexpath)
end
else
mexpath = [dynareroot '../mex/matlab/win64-7.8-8.2'];
mexpath = [dynareroot '../mex/matlab/win64-7.8-8.5'];
if exist(mexpath, 'dir')
addpath(mexpath)
end
......
matlab/dynare_estimation.m
View file @ 28eba237
......@@ -64,6 +64,11 @@ if (isnumeric(options_.mode_compute) && options_.mode_compute && options_.analyt
options_.analytic_derivation=1;
end
if options_.logged_steady_state
oo_.dr.ys=exp(oo_.dr.ys);
oo_.steady_state=exp(oo_.steady_state);
end
if nnobs > 1
for i=1:nnobs
......
matlab/dynare_estimation_1.m
View file @ 28eba237
......@@ -137,8 +137,7 @@ ub = bayestopt_.ub;
dr = oo_.dr;
%% load optimal_mh_scale parameter if previous run was with
%% mode_compute=6
% load optimal_mh_scale parameter if previous run was with mode_compute=6
mh_scale_fname = [M_.fname '_optimal_mh_scale_parameter.mat'];
if exist(mh_scale_fname)
if options_.mode_compute == 0
......@@ -147,7 +146,7 @@ if exist(mh_scale_fname)
clear tmp;
else
% remove the file if mode_compute ~= 0
delete('mh_scale_fname')
delete(mh_scale_fname)
end
end
......@@ -281,25 +280,22 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
epsilon = options_.gradient_epsilon;
% Change some options.
if isfield(options_,'optim_opt')
options_list = strsplit(options_.optim_opt,',');
number_of_options = length(options_list)/2;
o = 1;
while o<=number_of_options
switch strtrim(options_list{2*(o-1)+1})
case '''MaxIter'''
nit = str2num(options_list{2*(o-1)+2});
case '''InitialInverseHessian'''
H0 = eval(eval(options_list{2*(o-1)+2}));
case '''TolFun'''
crit = str2double(options_list{2*(o-1)+2});
case '''NumgradAlgorithm'''
numgrad = str2num(options_list{2*(o-1)+2});
case '''NumgradEpsilon'''
epsilon = str2double(options_list{2*(o-1)+2});
options_list = read_key_value_string(options_.optim_opt);
for i=1:rows(options_list)
switch options_list{i,1}
case 'MaxIter'
nit = options_list{i,2};
case 'InitialInverseHessian'
H0 = eval(options_list{i,2});
case 'TolFun'
crit = options_list{i,2};
case 'NumgradAlgorithm'
numgrad = options_list{i,2};
case 'NumgradEpsilon'
epsilon = options_list{i,2};
otherwise
warning(['csminwel: Unknown option (' options_list{2*(o-1)+1} ')!'])
warning(['csminwel: Unknown option (' options_list{i,1} ')!'])
end
o = o + 1;
end
end
% Set flag for analytical gradient.
......@@ -353,21 +349,19 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
gmhmaxlikOptions.varinit = 'prior';
end
if isfield(options_,'optim_opt')
options_list = strsplit(options_.optim_opt,',');
number_of_options = length(options_list)/2;
o = 1;
while o<=number_of_options
switch strtrim(options_list{2*(o-1)+1})
case '''NumberOfMh'''
gmhmaxlikOptions.iterations = str2num(options_list{2*(o-1)+2});
case '''ncov-mh'''
gmhmaxlikOptions.number = str2num(options_list{2*(o-1)+2});
case '''nscale'''
gmhmaxlikOptions.nscale = str2double(options_list{2*(o-1)+2});
case '''nclimb'''
gmhmaxlikOptions.nclimb = str2num(options_list{2*(o-1)+2});
case '''InitialCovarianceMatrix'''
switch eval(options_list{2*(o-1)+2})
options_list = read_key_value_string(options_.optim_opt);
for i=1:rows(options_list)
switch options_list{i,1}
case 'NumberOfMh'
gmhmaxlikOptions.iterations = options_list{i,2};
case 'ncov-mh'
gmhmaxlikOptions.number = options_list{i,2};
case 'nscale'
gmhmaxlikOptions.nscale = options_list{i,2};
case 'nclimb'
gmhmaxlikOptions.nclimb = options_list{i,2};
case 'InitialCovarianceMatrix'
switch options_list{i,2}
case 'previous'
if isempty(hh)
error('gmhmaxlik: No previous estimate of the Hessian matrix available!')
......@@ -375,19 +369,18 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
gmhmaxlikOptions.varinit = 'previous'
end
case {'prior', 'identity'}
gmhmaxlikOptions.varinit = eval(options_list{2*(o-1)+2});
gmhmaxlikOptions.varinit = options_list{i,2};
otherwise
error('gmhmaxlik: Unknown value for option ''InitialCovarianceMatrix''!')
end
case '''AcceptanceRateTarget'''
gmhmaxlikOptions.target = str2num(options_list{2*(o-1)+2});
case 'AcceptanceRateTarget'
gmhmaxlikOptions.target = options_list{i,2};
if gmhmaxlikOptions.target>1 || gmhmaxlikOptions.target<eps
error('gmhmaxlik: The value of option AcceptanceRateTarget should be a double between 0 and 1!')
end
otherwise
Warning(['gmhmaxlik: Unknown option (' options_list{2*(o-1)+1} ')!'])
warning(['gmhmaxlik: Unknown option (' options_list{i,1} ')!'])
end
o = o + 1;
end
end
% Evaluate the objective function.
......@@ -485,27 +478,24 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
% Dynare implementation of the simplex algorithm.
simplexOptions = options_.simplex;
if isfield(options_,'optim_opt')
options_list = strsplit(options_.optim_opt,',');
number_of_options = length(options_list)/2;
o = 1;
while o<=number_of_options
switch strtrim(options_list{2*(o-1)+1})
case '''MaxIter'''
simplexOptions.maxiter = str2num(options_list{2*(o-1)+2});
case '''TolFun'''
simplexOptions.tolerance.f = str2double(options_list{2*(o-1)+2});
case '''TolX'''
simplexOptions.tolerance.x = str2double(options_list{2*(o-1)+2});
case '''MaxFunEvals'''
simplexOptions.maxfcall = str2num(options_list{2*(o-1)+2});
case '''MaxFunEvalFactor'''
simplexOptions.maxfcallfactor = str2num(options_list{2*(o-1)+2});
case '''InitialSimplexSize'''
simplexOptions.delta_factor = str2double(options_list{2*(o-1)+2});
options_list = read_key_value_string(options_.optim_opt);
for i=1:rows(options_list)
switch options_list{i,1}
case 'MaxIter'
simplexOptions.maxiter = options_list{i,2};
case 'TolFun'
simplexOptions.tolerance.f = options_list{i,2};
case 'TolX'
simplexOptions.tolerance.x = options_list{i,2};
case 'MaxFunEvals'
simplexOptions.maxfcall = options_list{i,2};
case 'MaxFunEvalFactor'
simplexOptions.maxfcallfactor = options_list{i,2};
case 'InitialSimplexSize'
simplexOptions.delta_factor = options_list{i,2};
otherwise
warning(['simplex: Unknown option (' options_list{2*(o-1)+1} ')!'])
warning(['simplex: Unknown option (' options_list{i,1} ')!'])
end
o = o + 1;
end
end
[xparam1,fval,exitflag] = simplex_optimization_routine(objective_function,xparam1,simplexOptions,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
......@@ -515,23 +505,20 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
cmaesOptions = options_.cmaes;
% Modify defaults
if isfield(options_,'optim_opt')
options_list = strsplit(options_.optim_opt,',');
number_of_options = length(options_list)/2;
o = 1;
while o<=number_of_options
switch strtrim(options_list{2*(o-1)+1})
case '''MaxIter'''
cmaesOptions.MaxIter = str2num(options_list{2*(o-1)+2});
case '''TolFun'''
cmaesOptions.TolFun = str2double(options_list{2*(o-1)+2});
case '''TolX'''
cmaesOptions.TolX = str2double(options_list{2*(o-1)+2});
case '''MaxFunEvals'''
cmaesOptions.MaxFunEvals = str2num(options_list{2*(o-1)+2});
options_list = read_key_value_string(options_.optim_opt);
for i=1:rows(options_list)
switch options_list{i,1}
case 'MaxIter'
cmaesOptions.MaxIter = options_list{i,2};
case 'TolFun'
cmaesOptions.TolFun = options_list{i,2};
case 'TolX'
cmaesOptions.TolX = options_list{i,2};
case 'MaxFunEvals'
cmaesOptions.MaxFunEvals = options_list{i,2};
otherwise
warning(['cmaes: Unknown option (' options_list{2*(o-1)+1} ')!'])
warning(['cmaes: Unknown option (' options_list{i,1} ')!'])
end
o = o + 1;
end
end
warning('off','CMAES:NonfinitenessRange');
......@@ -542,30 +529,27 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
case 10
simpsaOptions = options_.simpsa;
if isfield(options_,'optim_opt')
options_list = strsplit(options_.optim_opt,',');
number_of_options = length(options_list)/2;
o = 1;
while o<=number_of_options
switch strtrim(options_list{2*(o-1)+1})
case '''MaxIter'''
simpsaOptions.MAX_ITER_TOTAL = str2num(options_list{2*(o-1)+2});
case '''TolFun'''
simpsaOptions.TOLFUN = str2double(options_list{2*(o-1)+2});
case '''TolX'''
tolx = str2double(options_list{2*(o-1)+2});
options_list = read_key_value_string(options_.optim_opt);
for i=1:rows(options_list)
switch options_list{i,1}
case 'MaxIter'
simpsaOptions.MAX_ITER_TOTAL = options_list{i,2};
case 'TolFun'
simpsaOptions.TOLFUN = options_list{i,2};
case 'TolX'
tolx = options_list{i,2};
if tolx<0
simpsaOptions = rmfield(simpsaOptions,'TOLX'); % Let cmaes choose the default.
simpsaOptions = rmfield(simpsaOptions,'TOLX'); % Let simpsa choose the default.
else
simpsaOptions.TOLX = tolx;
end
case '''EndTemparature'''
simpsaOptions.TEMP_END = str2double(options_list{2*(o-1)+2});
case '''MaxFunEvals'''
simpsaOptions.MAX_FUN_EVALS = str2num(options_list{2*(o-1)+2});
case 'EndTemparature'
simpsaOptions.TEMP_END = options_list{i,2};
case 'MaxFunEvals'
simpsaOptions.MAX_FUN_EVALS = options_list{i,2};
otherwise
warning(['simpsa: Unknown option (' options_list{2*(o-1)+1} ')!'])
warning(['simpsa: Unknown option (' options_list{i,1} ')!'])
end
o = o + 1;
end
end
simpsaOptionsList = options2cell(simpsaOptions);
......@@ -632,7 +616,6 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
[junk1, junk2, hh] = feval(objective_function,xparam1, ...
dataset_,options_,M_,estim_params_,bayestopt_,oo_);
options_.analytic_derivation = ana_deriv;
else
hh = reshape(hessian(objective_function,xparam1, ...
options_.gstep,dataset_,options_,M_,estim_params_,bayestopt_,oo_),nx,nx);
......@@ -806,7 +789,7 @@ if (any(bayestopt_.pshape >0 ) && options_.mh_replic) || ...
CutSample(M_, options_, estim_params_);
return
else
if ~options_.nodiagnostic && options_.mh_replic > 2000
if ~options_.nodiagnostic && options_.mh_replic>0
oo_= McMCDiagnostics(options_, estim_params_, M_,oo_);
end
%% Here i discard first half of the draws:
......@@ -837,8 +820,8 @@ if (any(bayestopt_.pshape >0 ) && options_.mh_replic) || ...
if options_.smoother || ~isempty(options_.filter_step_ahead) || options_.forecast
prior_posterior_statistics('posterior',dataset_);
end
xparam = get_posterior_parameters('mean');
M_ = set_all_parameters(xparam,estim_params_,M_);
xparam1 = get_posterior_parameters('mean');
M_ = set_all_parameters(xparam1,estim_params_,M_);
end
end
......
matlab/dynare_estimation_init.m
View file @ 28eba237
......@@ -85,7 +85,7 @@ if options_.lik_init == 1
if isempty(options_.qz_criterium)
options_.qz_criterium = 1-1e-6;
elseif options_.qz_criterium > 1-eps
error(['estimation: option qz_criterium is too large for estimating ' ...
error(['Estimation: option qz_criterium is too large for estimating ' ...
'a stationary model. If your model contains unit roots, use ' ...
'option diffuse_filter'])
end
......@@ -93,12 +93,6 @@ elseif isempty(options_.qz_criterium)
options_.qz_criterium = 1+1e-6;
end
% If the data are prefiltered then there must not be constants in the
% measurement equation of the DSGE model or in the DSGE-VAR model.
if options_.prefilter == 1
options_.noconstant = 1;
end
% Set options related to filtered variables.
if ~isequal(options_.filtered_vars,0) && isempty(options_.filter_step_ahead)
options_.filter_step_ahead = 1;
......@@ -130,13 +124,18 @@ if ~isempty(estim_params_) && ~isempty(options_.mode_file) && ~options_.mh_poste
number_of_estimated_parameters = length(xparam1);
mode_file = load(options_.mode_file);
if number_of_estimated_parameters>length(mode_file.xparam1)
% More estimated parameters than parameters in the mode file.
skipline()
disp(['The posterior mode file ' options_.mode_file ' has been generated using another specification of the model or another model!'])
disp(['Your mode file contains estimates for ' int2str(length(mode_file.xparam1)) ' parameters, while you are attempting to estimate ' int2str(number_of_estimated_parameters) ' parameters:'])
md = []; xd = [];
for i=1:number_of_estimated_parameters
id = strmatch(deblank(bayestopt_.name(i,:)),mode_file.parameter_names,'exact');
if isempty(id)
disp(['--> Estimated parameter ' bayestopt_.name{i} ' is not present in the loaded mod_file.'])
disp(['--> Estimated parameter ' bayestopt_.name{i} ' is not present in the loaded mode file (prior mean will be used, if possible).'])
else
xd = [xd; i];
md = [md; id];
end
end
for i=1:length(mode_file.xparam1)
......@@ -145,20 +144,30 @@ if ~isempty(estim_params_) && ~isempty(options_.mode_file) && ~options_.mh_poste
disp(['--> Parameter ' mode_file.parameter_names{i} ' is not estimated according to the current mod file.'])
end
end
error('Please change the mode_file option or the list of estimated parameters.')
if ~options_.mode_compute
% The posterior mode is not estimated.
error('Please change the mode_file option, the list of estimated parameters or set mode_compute>0.')
else
% The posterior mode is estimated, the Hessian evaluated at the mode is not needed so we set values for the parameters missing in the mode file using the prior mean.
if ~isempty(xd)
xparam1(xd) = mode_file.xparam1(md);
else
error('Please remove the mode_file option.')
end
end
elseif number_of_estimated_parameters<length(mode_file.xparam1)
% Less estimated parameters than parameters in the mode file.
skipline()
disp(['The posterior mode file ' options_.mode_file ' has been generated using another specification of the model or another model!'])
disp(['Your mode file contains estimates for ' int2str(length(mode_file.xparam1)) ' parameters, while you are attempting to estimate only ' int2str(number_of_estimated_parameters) ' parameters:'])
Id = [];
md = []; xd = [];
for i=1:number_of_estimated_parameters
id = strmatch(deblank(bayestopt_.name(i,:)),mode_file.parameter_names,'exact');
if isempty(id)
disp(['--> Estimated parameter ' deblank(bayestopt_.name(i,:)) ' is not present in the loaded mode file.'])
Id = [];
break
disp(['--> Estimated parameter ' deblank(bayestopt_.name(i,:)) ' is not present in the loaded mode file (prior mean will be used, if possible).'])
else
Id = [Id; id];
xd = [xd; i];
md = [md; id];
end
end
for i=1:length(mode_file.xparam1)
......@@ -167,59 +176,76 @@ if ~isempty(estim_params_) && ~isempty(options_.mode_file) && ~options_.mh_poste
disp(['--> Parameter ' mode_file.parameter_names{i} ' is not estimated according to the current mod file.'])
end
end
if isempty(Id)
% None of the estimated parameters are present in the mode_file.
error('Please change the mode_file option or the list of estimated parameters.')
else
% If possible, fix the mode_file.
if isequal(length(Id),number_of_estimated_parameters)
if ~options_.mode_compute
% The posterior mode is not estimated. If possible, fix the mode_file.
if isequal(length(xd),number_of_estimated_parameters)
disp('==> Fix mode file (remove unused parameters).')
mode_file.parameter_names = mode_file.parameter_names(Id,:);
mode_file.xparam1 = mode_file.xparam1(Id);
xparam1 = mode_file.xparam1(md);
if isfield(mode_file,'hh')
mode_file.hh = mode_file.hh(Id,Id);
hh = mode_file.hh(md,md);
end
else
error('Please change the mode_file option, the list of estimated parameters or set mode_compute>0.')
end
else
% The posterior mode is estimated, the Hessian evaluated at the mode is not needed so we set values for the parameters missing in the mode file using the prior mean.
if ~isempty(xd)
xparam1(xd) = mode_file.xparam1(md);
else
% None of the estimated parameters are present in the mode_file.
error('Please remove the mode_file option.')
end
end
else
% The number of declared estimated parameters match the number of parameters in the mode file.
% Check that the parameters in the mode file and according to the current mod file are identical.
if ~isfield(mode_file,'parameter_names')
disp(['The posterior mode file ' options_.mode_file ' has been generated using an older version of Dynare. It cannot be verified if it matches the present model. Proceed at your own risk.'])
mode_file.parameter_names=deblank(bayestopt_.name); %set names
end
if isequal(mode_file.parameter_names, bayestopt_.name)
% Ok! Nothing to do here.
xparam1 = mode_file.xparam1;
if isfield(mode_file,'hh')
hh = mode_file.hh;
end
else
skipline()
disp(['The posterior mode file ' options_.mode_file ' has been generated using another specification of the model or another model!'])
% Check if this only an ordering issue.
Id = [];
% Check if this only an ordering issue or if the missing parameters can be initialized with the prior mean.
md = []; xd = [];
for i=1:number_of_estimated_parameters
id = strmatch(deblank(bayestopt_.name(i,:)), mode_file.parameter_names,'exact');
if isempty(id)
disp(['--> Estimated parameter ' bayestopt_.name{i} ' is not present in the loaded mode file.'])
Id = [];
break
else
Id = [Id; id];
xd = [xd; i];
md = [md; id];
end
end
if isempty(Id)
% None of the estimated parameters are present in the mode_file.
error('Please change the mode_file option or the list of estimated parameters.')
else
% If possible, fix the mode_file.
if isequal(length(Id),number_of_estimated_parameters)
disp('==> Fix mode file (reorder the parameters).')
mode_file.parameter_names = mode_file.parameter_names(Id,:);
mode_file.xparam1 = mode_file.xparam1(Id);
if ~options_.mode_compute
% The posterior mode is not estimated
if isequal(length(xd), number_of_estimated_parameters)
% This is an ordering issue.
xparam1 = mode_file.xparam1(md);
if isfield(mode_file,'hh')
mode_file.hh = mode_file.hh(Id,Id);
hh = mode_file.hh(md,md);
end
else
error('Please change the mode_file option, the list of estimated parameters or set mode_compute>0.')
end
else
% The posterior mode is estimated, the Hessian evaluated at the mode is not needed so we set values for the parameters missing in the mode file using the prior mean.
if ~isempty(xd)
xparam1(xd) = mode_file.xparam1(md);
if isfield(mode_file,'hh')
hh(xd,xd) = mode_file.hh(md,md);
end
else
% None of the estimated parameters are present in the mode_file.
error('Please remove the mode_file option.')
end
end
end
xparam1 = mode_file.xparam1;
if isfield(mode_file,'hh')
hh = mode_file.hh;
end
skipline()
end
......@@ -252,7 +278,7 @@ end
if isempty(estim_params_)% If estim_params_ is empty (e.g. when running the smoother on a calibrated model)
if ~options_.smoother
error('ESTIMATION: the ''estimated_params'' block is mandatory (unless you are running a smoother)')
error('Estimation: the ''estimated_params'' block is mandatory (unless you are running a smoother)')
end
xparam1 = [];
bayestopt_.lb = [];
......@@ -336,11 +362,15 @@ oo_.dr.restrict_columns = [ic; length(k2)+(1:nspred-npred)'];
k3 = [];
k3p = [];
if options_.selected_variables_only
if options_.forecast > 0 && options_.mh_replic == 0 && ~options_.load_mh_file
fprintf('\nEstimation: The selected_variables_only option is incompatible with classical forecasts. It will be ignored.\n')
k3 = (1:M_.endo_nbr)';
k3p = (1:M_.endo_nbr)';
else
for i=1:size(var_list_,1)
k3 = [k3; strmatch(var_list_(i,:),M_.endo_names(dr.order_var,:), ...
'exact')];
k3p = [k3; strmatch(var_list_(i,:),M_.endo_names, ...
'exact')];
k3 = [k3; strmatch(var_list_(i,:),M_.endo_names(dr.order_var,:), 'exact')];
k3p = [k3; strmatch(var_list_(i,:),M_.endo_names, 'exact')];
end
end
else
k3 = (1:M_.endo_nbr)';
......@@ -391,7 +421,10 @@ else
bayestopt_.smoother_var_list);
end;
if options_.analytic_derivation,
if options_.analytic_derivation
if options_.lik_init == 3
error('analytic derivation is incompatible with diffuse filter')
end
options_.analytic_derivation = 1;
if ~(exist('sylvester3','file')==2),
dynareroot = strrep(which('dynare'),'dynare.m','');
......@@ -450,7 +483,7 @@ dataset_ = initialize_dataset(options_.datafile,options_.varobs,options_.first_o
options_.nobs = dataset_.info.ntobs;
% setting noconstant option
% setting steadystate_check_flag option
if options_.diffuse_filter
steadystate_check_flag = 0;
else
......@@ -470,4 +503,15 @@ if all(abs(oo_.steady_state(bayestopt_.mfys))<1e-9)
options_.noconstant = 1;
else
options_.noconstant = 0;
% If the data are prefiltered then there must not be constants in the
% measurement equation of the DSGE model or in the DSGE-VAR model.
if options_.prefilter
skipline()
disp('You have specified the option "prefilter" to demean your data but the')
disp('steady state of of the observed variables is non zero.')
disp('Either change the measurement equations, by centering the observed')
disp('variables in the model block, or drop the prefiltering.')
error('The option "prefilter" is inconsistent with the non-zero mean measurement equations in the model.')
end
end
matlab/dynare_identification.m
View file @ 28eba237
......@@ -207,15 +207,23 @@ else
name_tex = [cellstr(M_.exo_names_tex); cellstr(M_.param_names_tex)];
end
skipline()
disp(['==== Identification analysis ====' ]),
skipline()
if nparam<2,
options_ident.advanced=0;
advanced = options_ident.advanced;
disp('There is only one parameter to study for identitification.')
disp('The advanced option is re-set to 0.')
skipline()
end
options_ident = set_default_option(options_ident,'max_dim_cova_group',min([2,nparam-1]));
options_ident.max_dim_cova_group = min([options_ident.max_dim_cova_group,nparam-1]);
MaxNumberOfBytes=options_.MaxNumberOfBytes;
store_options_ident = options_ident;
skipline()
disp(['==== Identification analysis ====' ]),
skipline()
if iload <=0,
......@@ -326,8 +334,10 @@ if iload <=0,
disp('----------- ')
skipline()
return
end
else
parameters = 'Random_prior_params';
end
end
idehess_point.params=params;
% siH = idemodel_point.siH;
% siJ = idemoments_point.siJ;
......@@ -336,11 +346,11 @@ if iload <=0,
% normJ = max(abs(siJ)')';
% normLRE = max(abs(siLRE)')';
save([IdentifDirectoryName '/' M_.fname '_identif.mat'], 'idehess_point', 'idemoments_point','idemodel_point', 'idelre_point','store_options_ident')
save([IdentifDirectoryName '/' M_.fname '_' parameters '_identif.mat'], 'idehess_point', 'idemoments_point','idemodel_point', 'idelre_point','store_options_ident')
disp_identification(params, idemodel_point, idemoments_point, name, advanced);
if ~options_.nograph,
plot_identification(params,idemoments_point,idehess_point,idemodel_point,idelre_point,advanced,parameters,name,IdentifDirectoryName);
end
end
if SampleSize > 1,
skipline()
......@@ -522,7 +532,7 @@ if SampleSize > 1,
disp('Testing MC sample')
disp_identification(pdraws, idemodel, idemoments, name);
if ~options_.nograph,
plot_identification(pdraws,idemoments,idehess_point,idemodel,idelre,advanced,'MC sample - ',name, IdentifDirectoryName);
plot_identification(pdraws,idemoments,idehess_point,idemodel,idelre,advanced,'MC sample ',name, IdentifDirectoryName);
end
if advanced,
jcrit=find(idemoments.ino);
......
matlab/dynare_sensitivity.m
View file @ 28eba237
......@@ -23,6 +23,10 @@ function x0=dynare_sensitivity(options_gsa)
global M_ options_ oo_ bayestopt_ estim_params_
if options_.dsge_var
error('Identification does not support DSGE-VARs at the current stage')
end
fname_ = M_.fname;
if ~isfield(M_,'dname'),
M_.dname = M_.fname;
......@@ -156,6 +160,11 @@ options_gsa = set_default_option(options_gsa,'istart_rmse',options_.presample+1)
options_gsa = set_default_option(options_gsa,'alpha_rmse',0.001);
options_gsa = set_default_option(options_gsa,'alpha2_rmse',0);
if options_gsa.neighborhood_width,
options_gsa.pprior=0;
options_gsa.ppost=0;
end
if options_gsa.redform && options_gsa.neighborhood_width==0 && isempty(options_gsa.threshold_redform),
options_gsa.pprior=1;
options_gsa.ppost=0;
......@@ -281,10 +290,12 @@ if options_gsa.redform && ~isempty(options_gsa.namendo),% ...
% check existence of the SS_ANOVA toolbox
if isempty(options_gsa.threshold_redform) && ...
~(exist('gsa_sdp','file')==6 || exist('gsa_sdp','file')==2),
disp('Download Mapping routines at:')
disp('http://eemc.jrc.ec.europa.eu/softwareDYNARE-Dowload.htm')
disp(' ' )
error('Mapping routines missing!')
fprintf('\nThe "SS-ANOVA-R: MATLAB Toolbox for the estimation of Smoothing Spline ANOVA models with Recursive algorithms" is missing.\n')
fprintf('To obtain it, go to:\n\n')
fprintf('http://ipsc.jrc.ec.europa.eu/?id=790 \n\n')
fprintf('and follow the instructions there.\n')
fprintf('After obtaining the files, you need to unpack them and set a Matlab Path to those files.\n')
error('SS-ANOVA-R Toolbox missing!')
end
redform_map(OutputDirectoryName,options_gsa);
......
matlab/dynare_solve.m
View file @ 28eba237
......@@ -75,7 +75,7 @@ if ~isempty(i)
disp(' i) if all parameters occurring in these equations are defined')
disp(' ii) that no division by an endogenous variable initialized to 0 occurs')
info = 1;
x = NaN;
x = NaN(size(fvec));
return;
end
......
matlab/endogenous_prior_restrictions.m
View file @ 28eba237
......@@ -42,11 +42,11 @@ end
infos=[0 0];
varlist=Model.endo_names(DynareResults.dr.order_var,:);
varlist=varlist(DynareResults.dr.restrict_var_list,:);
T=1;
NT=1;
for j=1:size(endo_prior_restrictions.irf,1),
T=max(T,endo_prior_restrictions.irf{j,3});
NT=max(NT,endo_prior_restrictions.irf{j,3});
end
for t=1:T,
for t=1:NT,
RR = T^(t-1)*R;
for j=1:size(endo_prior_restrictions.irf,1),
if endo_prior_restrictions.irf{j,3}~=t,
......
matlab/ep/extended_path.m
View file @ 28eba237
......@@ -350,4 +350,7 @@ end% (while) loop over t
dyn_waitbar_close(hh);
oo_.endo_simul = oo_.steady_state;
if ~nargout
oo_.endo_simul = [initial_conditions, time_series];
dyn2vec;
end
\ No newline at end of file