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
  • ebenetce/dynare
  • chskcau/dynare-doc-fixes
28 results
Select Git revision
Loading items
Show changes
Show whitespace changes
Inline Side-by-side
Showing
with 904 additions and 47 deletions
matlab/+hank/compute_steady_state_residuals.m 0 → 100644
View file @ a29fb27b
% Copyright © 2025 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/>.
%
% COMPUTE_STEADY_STATE_RESIDUALS Computes residuals for a heterogeneous-agent model at the steady state
%
% This function evaluates the residuals of both aggregate and individual (heterogeneous-agent)
% equilibrium conditions given a steady-state object stored in `oo_.hank`. It uses previously
% computed policy functions and basis matrices to evaluate the model’s residuals pointwise.
%
% INPUTS
% M_ [struct] : Dynare model structure
% oo_ [struct] : Dynare output structure. Must contain the fields:
% - oo_.hank.ss : steady-state structure (see initialize_steady_state)
% - oo_.hank.sizes : tensor-product grid sizes
% - oo_.hank.mat : basis and interpolation matrices
%
% OUTPUTS
% F [matrix] : Residuals of the heterogeneous-agent equations (size: H_.endo_nbr × N_sp),
% where N_sp is the number of points in the full state space.
% G [vector] : Residuals of the aggregate equilibrium equations (size: M_.endo_nbr × 1)
%
% NOTES
% - For the aggregate block:
% - The individual state distribution `ss.d.hist` is reordered to match the policy grid.
% - Aggregate variables defined as expectations over the distribution are recomputed
% using `mat.pol.x_bar_dash` and the `Phi` basis matrix.
% - The residuals are evaluated using the sparse dynamic function `.sparse.dynamic_resid`.
%
% - For the heterogeneous block:
% - Constructs the full stacked endogenous (`yh`) and exogenous (`xh`) state vectors
% at times t-1, t, and t+1 using the stored policy functions and interpolation structure.
% - Calls `.sparse.dynamic_het1_resid` at each point in the tensor-product state space.
%
% dynamic_resid, dynamic_het1_resid
function [F,G] = compute_steady_state_residuals(M_, oo_)
ss = oo_.hank.ss;
sizes = oo_.hank.sizes;
mat = oo_.hank.mat;
% Rearrange ss.d.hist so that the state and shock orders is similar to the one
% used for policy functions. order_shocks(i) contains the position of
% ss.pol.shocks(i) in ss.d.shocks. In other words, ss.d.shocks(order_shocks)
% equals ss.pol.shocks. A similar logic applies to order_states.
order_shocks = cellfun(@(x) find(strcmp(x,ss.d.shocks)), ss.pol.shocks);
order_states = cellfun(@(x) find(strcmp(x,ss.d.states)), ss.pol.states);
dim_states = cellfun(@(s) sizes.d.states.(s), ss.d.states);
dim_shocks = cellfun(@(s) sizes.shocks.(s), ss.d.shocks);
hist = reshape(ss.d.hist, [dim_shocks dim_states]);
hist = permute(hist, [order_shocks sizes.n_e+order_states]);
% Compute aggregated heterogeneous variables
Ix = mat.pol.x_bar_dash*mat.Phi*reshape(hist,[],1);
% Computing aggregate residuals
% - Extended aggregate endogenous variable vector at the steady state - %
y = NaN(3*M_.endo_nbr,1);
for i=1:M_.orig_endo_nbr
var = M_.endo_names{i};
y(i) = ss.agg.(var);
y(M_.endo_nbr+i) = ss.agg.(var);
y(2*M_.endo_nbr+i) = ss.agg.(var);
end
% - Variables resulting from an aggregation operator - %
iIx = zeros(size(M_.heterogeneity_aggregates,1),1);
ix = zeros(size(M_.heterogeneity_aggregates,1),1);
for i=1:size(M_.heterogeneity_aggregates,1)
if (M_.heterogeneity_aggregates{i,1} == 'sum')
iIx(i) = M_.aux_vars(M_.heterogeneity_aggregates{i,2}).endo_index;
ix(i) = M_.heterogeneity_aggregates{i,3};
end
end
yagg = Ix(ix);
y(iIx) = yagg;
y(M_.endo_nbr+iIx) = yagg;
y(2*M_.endo_nbr+iIx) = yagg;
% - Exogenous variables - %
x = zeros(M_.exo_nbr,1);
% - Call to the residual function - %
G = feval([M_.fname '.sparse.dynamic_resid'], y, x, M_.params, [], yagg);
% Compute heterogeneous residuals
N_sp = sizes.N_e*sizes.pol.N_a;
H_ = M_.heterogeneity(1);
% Heterogeneous endogenous variable vector
yh = NaN(3*H_.endo_nbr, N_sp);
% Heterogeneous exogenous variable vector
xh = NaN(H_.exo_nbr, N_sp);
% Get the H_.endo_names indices for the steady-state ordering ss.pol.shocks and ss.pol.states
if isfield(ss.shocks, 'Pi')
order_shocks = cellfun(@(x) find(strcmp(x,H_.endo_names)), ss.pol.shocks);
else
order_shocks = cellfun(@(x) find(strcmp(x,H_.exo_names)), ss.pol.shocks);
end
% Get the H_.endo_names indices for the steady-state ordering ss.pol.shocks and ss.pol.states
order_states = cellfun(@(x) find(strcmp(x,H_.endo_names)), ss.pol.states);
% Get the indices of originally declared pol values in H_.endo_names
pol_names = H_.endo_names;
pol_names(ismember(H_.endo_names, ss.pol.shocks)) = [];
pol_ind = cellfun(@(x) find(strcmp(x,H_.endo_names)), pol_names);
% t-1
yh(order_states,:) = mat.pol.sm(1:sizes.n_a,:);
% t
yh(H_.endo_nbr+pol_ind,:) = mat.pol.x_bar;
if isfield(ss.shocks, 'Pi')
yh(H_.endo_nbr+order_shocks,:) = mat.pol.sm(sizes.n_a+1:end,:);
else
xh(order_shocks,:) = mat.pol.sm(sizes.n_a+1:end,:);
end
% t+1
yh(2*H_.endo_nbr+pol_ind,:) = mat.pol.x_bar_dash*mat.Phi_tilde_e;
% Call the residual function
F = NaN(H_.endo_nbr, N_sp);
for j=1:N_sp
F(:,j) = feval([M_.fname '.sparse.dynamic_het1_resid'], y, x, M_.params, [], yh(:,j), xh(:,j), []);
end
end
\ No newline at end of file
matlab/+hank/initialize_steady_state.m 0 → 100644
View file @ a29fb27b
% Copyright © 2025 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/>.
%
% INITIALIZE_STEADY_STATE Wrapper to validate, complete, and store steady-state input for HANK models.
%
% This function performs validation of the user-provided steady-state structure `ss`, constructs
% interpolation objects and basis matrices, computes missing policy functions for complementarity
% multipliers (if not supplied), and populates the global `oo_.hank` structure with the results.
%
% INPUTS
% M_ [struct] : Dynare model structure
% options_ [struct] : Dynare options structure
% oo_ [struct] : Dynare output structure to which results will be written
% ss [struct] : User-supplied steady-state structure (partial or complete)
% OUTPUTS
% oo_ [struct] : Updated Dynare output structure containing:
% - oo_.hank.ss : validated and completed steady-state structure
% - oo_.hank.sizes : grid size structure
% - oo_.hank.mat : interpolation and policy function matrices
function oo_ = initialize_steady_state(M_, options_, oo_, ss)
assert(isscalar(M_.heterogeneity), 'Heterogeneous-agent models with more than one heterogeneity dimension are not allowed yet!');
% Check steady-state input
[out_ss, sizes] = hank.check_steady_state_input(M_, options_, ss);
% Build useful basis and interpolation matrices
H_ = M_.heterogeneity(1);
mat = build_grids_and_interpolation_matrices(out_ss, sizes);
% Compute policy matrices without the complementarity conditions multipliers
[mat.pol.x_bar, mat.pol.x_bar_dash] = compute_pol_matrices(H_, out_ss, sizes, mat);
% Compute the policy values of complementarity conditions multipliers
mult_values = compute_pol_mcp_mul(M_, out_ss, sizes, mat);
% Update ss.pol.values accordingly
f = fieldnames(mult_values);
for i=1:numel(f)
out_ss.pol.values.(f{i}) = mult_values.(f{i});
end
% Update the policy matrices
[mat.pol.x_bar, mat.pol.x_bar_dash] = compute_pol_matrices(H_, out_ss, sizes, mat);
% Store the results in oo_
oo_.hank.ss = out_ss;
oo_.hank.sizes = sizes;
oo_.hank.mat = mat;
end
\ No newline at end of file
matlab/+hank/private/build_grids_and_interpolation_matrices.m 0 → 100644
View file @ a29fb27b
% Copyright © 2025 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/>.
%
% BUILD_GRIDS_AND_INTERPOLATION_MATRICES Constructs interpolation and basis matrices
% for evaluating policy functions and distributions in heterogeneous-agent models.
%
% This function builds all numerical structures needed to evaluate policy residuals,
% compute expectations, and perform interpolation over the state space. It constructs:
% - tensor-product basis matrices (`Phi`)
% - expectation operators (`Phi_tilde_e`)
% - grid point matrices (`pol.sm`, `d.sm`)
% - LU decomposition of the interpolation basis for efficient linear solves.
%
% INPUTS
% ss [struct] : Steady-state structure validated by `check_steady_state_input`
% sizes [struct] : Structure containing grid dimensions and counts
%
% OUTPUT
% mat [struct] : Structure containing matrices used in solution and simulation:
% - Phi : basis matrix for evaluating interpolated functions over distribution grid
% - Phi_tilde : expanded identity matrix for state variables (used in solving)
% - Phi_tilde_e : expectation operator incorporating both interpolation and transition shocks
% - Mu : Kronecker product of exogenous transition or quadrature weights
% - pol.sm : full tensor grid of state and shock combinations for policy evaluation
% - d.sm : full tensor grid of state and shock combinations for distribution evaluation
% - L_Phi_tilde,
% U_Phi_tilde,
% P_Phi_tilde : LU decomposition of `Phi_tilde` (used to project policy values)
%
% Notes:
% - The structure `ss` must contain valid grids and values in `ss.pol.values`, `ss.pol.states`,
% `ss.pol.shocks`, `ss.d.grids`, and `ss.shocks`.
% - Assumes linear interpolation and separable basis construction across state variables.
% - `Phi_tilde_e` applies the expected operator `E_t[f(x')]` for the simulation step.
function mat = build_grids_and_interpolation_matrices(ss, sizes)
mat = struct;
% Compute Φ and ̃Φ with useful basis matrices
Phi = speye(sizes.N_e);
Phi_tilde = speye(sizes.N_e);
B = struct;
B_tilde_bar = struct;
for i=sizes.n_a:-1:1
s = ss.pol.states{i};
[ind, w] = find_bracket_linear_weight(ss.pol.grids.(s), ss.d.grids.(s));
B.(s) = lin_spli(ind, w, sizes.pol.states.(s));
[ind, w] = find_bracket_linear_weight(ss.pol.grids.(s), reshape(ss.pol.values.(s), 1, []));
B_tilde_bar.(s) = lin_spli(ind, w, sizes.pol.states.(s));
Phi = kron(B.(s), Phi);
Phi_tilde = kron(speye(sizes.pol.states.(s)), Phi_tilde);
end
mat.Phi = Phi;
mat.Phi_tilde = Phi_tilde;
% Compute ̃Φᵉ and Mu
N_sp = sizes.pol.N_a*sizes.N_e;
hadamard_prod = ones(sizes.pol.N_a,N_sp);
prod_a_1_to_jm1 = 1;
prod_a_jp1_to_n_a = sizes.pol.N_a;
for j=1:sizes.n_a
s = ss.pol.states{j};
prod_a_jp1_to_n_a = prod_a_jp1_to_n_a/sizes.pol.states.(s);
hadamard_prod = hadamard_prod .* kron(kron(ones(prod_a_1_to_jm1,1), B_tilde_bar.(s)), ones(prod_a_jp1_to_n_a,1));
prod_a_1_to_jm1 = prod_a_1_to_jm1*sizes.pol.states.(s);
end
Mu = eye(1);
if isfield(ss.shocks, 'Pi')
for i=sizes.n_e:-1:1
e = ss.pol.shocks{i};
Mu = kron(ss.shocks.Pi.(e), Mu);
end
Phi_tilde_e = kron(hadamard_prod, ones(sizes.N_e,1)).*kron(ones(sizes.pol.N_a), Mu');
else
for i=sizes.n_e:-1:1
e = ss.pol.shocks{i};
Mu = kron(ss.shocks.w.(e), Mu);
end
Phi_tilde_e = kron(hadamard_prod, Mu);
end
mat.Mu = Mu;
mat.Phi_tilde_e = Phi_tilde_e;
% Compute state and shock grids for the policy and distribution state grids
mat.pol.sm = set_state_matrix(ss, sizes, "pol");
mat.d.sm = set_state_matrix(ss, sizes, "d");
% LU decomposition of Phi_tilde
[L_Phi_tilde,U_Phi_tilde,P_Phi_tilde] = lu(Phi_tilde);
mat.L_Phi_tilde = L_Phi_tilde;
mat.U_Phi_tilde = U_Phi_tilde;
mat.P_Phi_tilde = P_Phi_tilde;
end
\ No newline at end of file
matlab/perfect-foresight-models/linear_approximation_accuracy.m → matlab/+hank/private/check_consistency.m
View file @ a29fb27b
function err = linear_approximation_accuracy(options_, M_, oo_) % Copyright © 2025 Dynare Team
% Evaluates the accuracy of the linear approximation when solving perfect foresight models, by
% reporting the max absolute value of the dynamic residuals.
%
% INPUTS
% - options_ [struct] contains various options.
% - M_ [struct] contains a description of the model.
% - oo_ [struct] contains results.
%
% OUTPUTS
% - err [double] n*1 vector, evaluation of the accuracy (n is the number of equations).
% Copyright © 2015-2017 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
...@@ -26,30 +14,24 @@ function err = linear_approximation_accuracy(options_, M_, oo_) ...@@ -26,30 +14,24 @@ function err = linear_approximation_accuracy(options_, M_, oo_)
% %
% You should have received a copy of the GNU General Public License % You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>. % along with Dynare. If not, see <https://www.gnu.org/licenses/>.
%
lead_lag_incidence = M_.lead_lag_incidence; % Check consistency between two sets of variable names.
%
ny = M_.endo_nbr; % INPUTS
% - f1 [cell array of strings]: list of variable names to be checked for consistency
maximum_lag = M_.maximum_lag; % - f2 [cell array of strings]: reference list of valid variable names
% - f1_name [char]: name of the first set (for error messages)
periods = options_.periods; % - f2_name [char]: name of the second set (for error messages)
steady_state = oo_.steady_state; %
params = M_.params; % OUTPUTS
endo_simul = oo_.endo_simul; % - (none) This function throws an error if there are elements in `f1` that are not present in `f2`.
exo_simul = oo_.exo_simul; %
% DESCRIPTION
model_dynamic = str2func([M_.fname,'.dynamic']); % Checks that all elements of `f1` are included in `f2`. If not, throws a descriptive error
% mentioning the missing variables and the corresponding structure names for easier debugging.
residuals = zeros(ny,periods); function check_consistency(f1, f2, f1_name, f2_name)
f1_in_f2 = ismember(f1,f2);
Y = endo_simul(:); if ~all(f1_in_f2)
error('Misspecified steady-state input `ss`: the following variables are mentioned in `%s`, but are missing in `%s`: %s', f1_name, f2_name, strjoin(f1(~f1_in_f2)));
i_cols = find(lead_lag_incidence')+(maximum_lag-1)*ny; end
for it = (maximum_lag+1):(maximum_lag+periods)
residuals(:,it-1) = model_dynamic(Y(i_cols), exo_simul, params, steady_state,it);
i_cols = i_cols + ny;
end end
\ No newline at end of file
err = transpose(max(abs(transpose(residuals))));
\ No newline at end of file
matlab/+hank/private/check_fun.m 0 → 100644
View file @ a29fb27b
% Copyright © 2025 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/>.
%
% Apply a validation function to a set of fields within a structure.
%
% INPUTS
% - f [structure]: structure containing fields to be validated
% - f_name [char]: name of the structure (for error messages)
% - symbs [cell array of strings]: list of field names to validate
% - fun [function handle]: validation function applied to each field
% - text [char]: descriptive text for the error message in case of validation failure
% - eq [optional, scalar]: if provided, the output of `fun` is compared to `eq`
%
% OUTPUTS
% - (none) This function throws an error if any field fails the validation.
%
% DESCRIPTION
% For each symbol in `symbs`, applies the validation function `fun` to the
% corresponding field in `f`. If an optional `eq` argument is provided, the
% output of `fun` is compared to `eq`. Throws a descriptive error listing all
% fields that do not satisfy the validation criteria.
function check_fun(f, f_name, symbs, fun, text, eq)
elt_check = cellfun(@(s) fun(f.(s)), symbs);
if nargin == 6
elt_check = elt_check == eq;
end
if ~all(elt_check)
error('Misspecified steady-state input `ss`: the following fields in `%s` %s: %s.', f_name, text, strjoin(symbs(~elt_check)));
end
end
\ No newline at end of file
matlab/+hank/private/check_isfield.m 0 → 100644
View file @ a29fb27b
% Copyright © 2025 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/>.
%
% Check the existence of a field in a structure.
%
% INPUTS
% - f [char]: name of the field to check
% - s [structure]: structure in which the field should be present
% - f_name [optional, char]: name to display in the error message (defaults to `f`)
% - details [optional, char]: additional details to append to the error message
%
% OUTPUTS
% - (none) This function throws an error if the specified field is missing.
%
% DESCRIPTION
% Checks whether the field `f` exists in the structure `s`.
% If not, throws an error indicating the missing field.
% If provided, `details` is appended to the error message for additional context.
function check_isfield(f, s, f_name, details)
if nargin < 4
details = '';
end
if nargin < 3
f_name = f;
end
if ~isfield(s, f)
error('Misspecified steady-state input `ss`: the `%s` field is missing.%s', f_name, details);
end
end
\ No newline at end of file
matlab/+hank/private/check_isstruct.m 0 → 100644
View file @ a29fb27b
% Copyright © 2025 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/>.
%
% Check that an input is a structure.
%
% INPUTS
% - f [any type]: variable to check
% - f_name [char]: name to display in the error message
%
% OUTPUTS
% - (none) This function throws an error if the input is not a structure.
%
% DESCRIPTION
% Checks whether `f` is a structure.
% If not, throws an error indicating that the field `f_name` should be a structure.
function check_isstruct(f, f_name)
if ~isstruct(f)
error('Misspecified steady-state input `ss`: the `%s` field should be a structure.', f_name);
end
end
\ No newline at end of file
matlab/+hank/private/check_missingredundant.m 0 → 100644
View file @ a29fb27b
% Copyright © 2025 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/>.
%
% Check for missing and redundant fields in a structure.
%
% INPUTS
% - f [structure]: structure to check
% - f_name [char]: name to display in error/warning messages
% - symbs [cell array of char]: list of expected field names
% - nowarningredundant [logical]: if true, suppresses warnings about redundant fields
%
% OUTPUTS
% - (none) This function throws an error if any expected field is missing,
% and issues a warning if redundant fields are found (unless nowarningredundant is true).
%
% DESCRIPTION
% Verifies that all fields listed in `symbs` are present in the structure `f`.
% Throws an error if any expected field is missing.
% If `nowarningredundant` is false, checks if `f` contains fields not listed
% in `symbs`, and issues a warning if redundant fields are found.
function check_missingredundant(f, f_name, symbs, nowarningredundant)
fields = fieldnames(f);
symbs_in_fields = ismember(symbs, fields);
if ~all(symbs_in_fields)
error('Misspecified steady-state input `ss`. The following fields are missing in `%s`: %s.', f_name, strjoin(symbs(~symbs_in_fields)));
end
if ~nowarningredundant
fields_in_symbs = ismember(fields, symbs);
if ~all(fields_in_symbs)
warning('Steady-state input `ss`. The following fields are redundant in `%s`: %s.', f_name, strjoin(fields(~fields_in_symbs)));
end
end
end
\ No newline at end of file
matlab/+hank/private/check_permutation.m 0 → 100644
View file @ a29fb27b
% Copyright © 2025 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/>.
%
% Check and returns a permutation order for variables.
%
% INPUTS
% - s [structure]: structure containing the permutation information
% - f_name [char]: name of the field within `s` that should contain the permutation list
% - s_name [char]: full name of the structure `s` for error messages
% - symbs [cell array of char]: list of expected variable names
%
% OUTPUTS
% - out_order [cell array or string array]: permutation order of variables
%
% DESCRIPTION
% Checks that the field `f_name` exists in structure `s` and contains a valid
% permutation (i.e., a cell array of character vectors or a string array).
% If the field is missing, returns the default order given by `symbs`.
% Throws an error if the specified variables are not consistent with `symbs`.
function out_order = check_permutation(s, f_name, s_name, symbs)
if ~isfield(s, f_name)
out_order = symbs;
else
f = s.(f_name);
if ~isstring(f) && ~iscellstr(f)
error('Misspecified steady-state input `ss`: the `%s.%s` field should be a cell array of character vectors or a string array.', s_name, f_name);
end
err_var = setdiff(f, symbs);
if ~isempty(err_var)
error('Misspecified steady-state input `ss`: the set of variables of the `%s.%s` field is not consistent with the information in M_ and the other fields in the steady-state structure `ss`. Problematic variables: %s.', s_name, f_name, strjoin(err_var));
end
out_order = f;
end
end
\ No newline at end of file
matlab/+hank/private/compute_pol_matrices.m 0 → 100644
View file @ a29fb27b
% Copyright © 2025 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/>.
%
% COMPUTE_POL_MATRICES Construct interpolated policy function matrices for state transitions.
%
% Given a steady-state structure `ss` with discretized policy functions, this function builds:
% - `x_bar`: matrix of policy function values (flattened over the joint state space)
% - `x_bar_dash`: transformed version of `x_bar` projected onto the basis used for expectations
%
% INPUTS
% H_ [struct] : Heterogeneous agent substructure from Dynare’s `M_`
% ss [struct] : Validated steady-state structure
% sizes [struct] : Structure with grid sizes
% mat [struct] : Structure containing transformation matrices (e.g., `U_Phi_tilde`)
%
% OUTPUTS
% x_bar [matrix] : Matrix of policy function values (num_vars x N_e * N_a)
% x_bar_dash [matrix] : Projected matrix in basis representation for expectations
function [x_bar, x_bar_dash] = compute_pol_matrices(H_, ss, sizes, mat)
% Computing x_bar
N_sp = sizes.N_e*sizes.pol.N_a;
x_names = H_.endo_names;
x_names(ismember(x_names, ss.pol.shocks)) = [];
x_bar = NaN(numel(x_names), N_sp);
for i=1:numel(x_names)
x = x_names{i};
if isfield(ss.pol.values, x)
x_bar(i,:) = reshape(ss.pol.values.(x),1,[]);
end
end
% Computing x_bar^#
x_bar_dash = x_bar/mat.U_Phi_tilde;
x_bar_dash = x_bar_dash/mat.L_Phi_tilde;
x_bar_dash = x_bar_dash*mat.P_Phi_tilde;
end
\ No newline at end of file
matlab/+hank/private/compute_pol_mcp_mul.m 0 → 100644
View file @ a29fb27b
% Copyright © 2025 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/>.
%
% COMPUTE_POL_MCP_MUL Compute policy function values for multipliers associated with binding
% complementarity constraints in HANK models.
%
% This routine evaluates residuals at grid points where inequality constraints are binding, and
% infers the corresponding values of Lagrange multipliers from the model’s dynamic equations. It
% distinguishes lower-bound and upper-bound constraints, and only evaluates residuals where needed.
%
% INPUTS
% M_ [struct] : Dynare model structure
% ss [struct] : Steady-state input validated by `check_steady_state_input`
% sizes [struct] : Structure containing grid sizes, as returned by `check_steady_state_input`
% mat [struct] : Structure with interpolation matrices and policy matrices
%
% OUTPUT
% mult_values [struct] : Structure extending `ss.pol.values` with computed policy function
% values for multipliers (e.g. `MULT_L_a`, `MULT_U_c`)
%
% Internally, this function uses `M_.fname.sparse.dynamic_het1_resid` to evaluate residuals,
% and relies on complementarity mappings from `M_.fname.dynamic_het1_complementarity_conditions`.
function mult_values = compute_pol_mcp_mul(M_, ss, sizes, mat)
% Get complementarity condition bounds
[lb, ub] = feval(sprintf('%s.dynamic_het1_complementarity_conditions', M_.fname), M_.params);
% Get constrained variables indices and names
H_ = M_.heterogeneity(1);
lb_constrained_var_names = H_.endo_names(lb ~= -Inf);
ub_constrained_var_names = H_.endo_names(ub ~= Inf);
lb = lb(lb ~= -Inf);
ub = ub(ub ~= Inf);
% Get the associated multipliers
lb_mult = arrayfun(@(x) "MULT_L_"+x,lb_constrained_var_names);
ub_mult = arrayfun(@(x) "MULT_U_"+x,ub_constrained_var_names);
mult_ind = arrayfun(@(x) x.endo_index, H_.aux_vars);
mult_names = H_.endo_names(mult_ind);
% Get the associated equations
eq_nbr = arrayfun(@(x) x.eq_nbr, H_.aux_vars);
eq_nbr_lb = arrayfun(@(x) eq_nbr(x == mult_names), lb_mult);
eq_nbr_ub = arrayfun(@(x) eq_nbr(x == mult_names), ub_mult);
% Initialize the multipliers policy functions
mult_values = struct;
for i=1:numel(lb_mult)
mult_name = lb_mult(i);
if ~isfield(ss.pol.values, mult_name)
mult_values.(mult_name) = zeros(sizes.N_e, sizes.pol.N_a);
end
end
for i=1:numel(ub_mult)
mult_name = ub_mult(i);
if ~isfield(ss.pol.values, mult_name)
mult_values.(mult_name) = zeros(sizes.N_e, sizes.pol.N_a);
end
end
% Get the discretized policy function values of multipliers
% Extended aggregate endogenous variable vector at the steady state
y = zeros(3*M_.endo_nbr,1);
for i=1:M_.orig_endo_nbr
var = M_.endo_names{i};
y(i) = ss.agg.(var);
y(M_.endo_nbr+i) = ss.agg.(var);
y(2*M_.endo_nbr+i) = ss.agg.(var);
end
% Aggregate exogenous variable vector at the steady state
x = zeros(M_.exo_nbr,1);
% Heterogeneous endogenous variable vector
yh = zeros(3*H_.endo_nbr,1);
% Heterogeneous exogenous variable vector
xh = zeros(H_.exo_nbr,1);
% Get the H_.endo_names indices for the steady-state ordering ss.pol.shocks and ss.pol.states
if isfield(ss.shocks, 'Pi')
order_shocks = cellfun(@(x) find(strcmp(x,H_.endo_names)), ss.pol.shocks);
else
order_shocks = cellfun(@(x) find(strcmp(x,H_.exo_names)), ss.pol.shocks);
end
% Get the H_.endo_names indices for the steady-state ordering ss.pol.shocks and ss.pol.states
order_states = cellfun(@(x) find(strcmp(x,H_.endo_names)), ss.pol.states);
% Get the indices of declared ss.pol.values element in H_.endo_names
all_pol_names = H_.endo_names;
all_pol_names(ismember(all_pol_names, ss.pol.shocks)) = [];
all_pol_ind = cellfun(@(x) find(strcmp(x,H_.endo_names)), all_pol_names);
% Get the indices of required declared pol values in H_.endo_names
min_pol_names = H_.endo_names(1:H_.orig_endo_nbr);
min_pol_names(ismember(min_pol_names, ss.pol.shocks)) = [];
min_pol_ind = cellfun(@(x) find(strcmp(x,H_.endo_names)), min_pol_names);
% For each lower-bound constrained variables var, find the indices of values for
% which ss.pol.values <= lb. For each of these indices, store minus the residual
% of the equation associated with var. Proceed in a similar way for
% upper-bound constrained variables, but store plus the residual.
for i=1:numel(lb_constrained_var_names)
var = lb_constrained_var_names{i};
mult_name = lb_mult(i);
if ~isfield(ss.pol.values, mult_name)
ind = find(ss.pol.values.(var) <= lb(i));
eq = eq_nbr_lb(i);
for j=1:numel(ind)
ind_sm = ind(j);
% t-1
yh(order_states) = mat.pol.sm(1:sizes.n_a, ind_sm);
% t
yh(H_.endo_nbr+min_pol_ind) = cellfun(@(x) ss.pol.values.(x)(ind_sm), min_pol_names);
if isfield(ss.shocks, 'Pi')
yh(H_.endo_nbr+order_shocks) = mat.pol.sm(sizes.n_a+1:end,ind_sm);
else
xh(order_shocks) = mat.pol.sm(sizes.n_a+1:end,ind_sm);
end
% t+1
yh(2*H_.endo_nbr+all_pol_ind) = mat.pol.x_bar_dash*mat.Phi_tilde_e(:,ind_sm);
% Residual function call
r = feval([M_.fname '.sparse.dynamic_het1_resid'], y, x, M_.params, [], yh, xh, []);
mult_values.(mult_name)(ind_sm) = -r(eq);
end
end
end
for i=1:numel(ub_constrained_var_names)
var = ub_constrained_var_names{i};
mult_name = ub_mult(i);
if ~isfield(ss.pol.values, mult_name)
ind = find(ss.pol.values.(var) >= ub(i));
eq = eq_nbr_ub(i);
for j=1:numel(ind)
ind_sm = ind(j);
% t-1
yh(order_states) = mat.pol.sm(1:sizes.n_a, ind_sm);
% t
yh(H_.endo_nbr+min_pol_ind) = cellfun(@(x) ss.pol.values.(x)(ind_sm), min_pol_names);
if isfield(ss.shocks, 'Pi')
yh(H_.endo_nbr+order_shocks) = mat.pol.sm(sizes.n_a+1:end,ind_sm);
else
xh(order_shocks) = mat.pol.sm(sizes.n_a+1:end,ind_sm);
end
% t+1
yh(2*H_.endo_nbr+all_pol_ind) = mat.pol.x_bar_dash*mat.Phi_tilde_e(:,ind_sm);
% Residual function call
r = feval([M_.fname '.sparse.dynamic_het1_resid'], y, x, M_.params, [], yh, xh, []);
mult_values.(mult_name)(ind_sm) = r(eq);
end
end
end
end
\ No newline at end of file
matlab/+hank/private/lin_spli.m 0 → 100644
View file @ a29fb27b
% Copyright © 2025 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/>.
%
% LIN_SPLI Constructs a sparse linear interpolation matrix for 1D bracketed data.
%
% This function builds a sparse matrix `B` of size m×n that maps `n` query points,
% located between known grid points, to the surrounding `m`-dimensional grid
% using linear interpolation weights.
%
% It is typically used in Dynare HANK routines after calling `bracket_linear_weight`
% (via MEX or native MATLAB) to form interpolation matrices over endogenous grids.
%
% INPUTS
% ind [n×1 int32] : Lower bracket indices for each query (i.e., index `ilow` s.t.
% x(ind(i)) ≤ xq(i) ≤ x(ind(i)+1))
%
% w [n×1 double] : Linear interpolation weights (relative position in interval):
% w(i) = (x(ind(i)+1) - xq(i)) / (x(ind(i)+1) - x(ind(i)))
%
% m [int scalar] : Size of the full grid (number of basis functions, i.e., length of `x`)
%
% OUTPUT
% B [m×n sparse] : Sparse interpolation matrix such that:
% f_interp = B * f_grid
% where `f_grid` is a column vector of function values on the grid.
%
% NOTES
% - Each column of `B` contains exactly two nonzero entries: `w(i)` and `1 - w(i)`
% - Assumes linear interpolation over uniform or non-uniform 1D grids
function B = lin_spli(ind, w, m)
i = [ind;ind+1];
n = length(ind);
j = repmat(int32(1:n)',2,1);
v = [w;1-w];
B = sparse(i,j,v,m,n);
end
\ No newline at end of file
matlab/+hank/private/set_state_matrix.m 0 → 100644
View file @ a29fb27b
% Copyright © 2025 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/>.
%
% SET_STATE_MATRIX Builds the full tensor-product grid of state and shock variables.
%
% This function constructs a matrix containing all combinations of state and shock values
% used for evaluating policy functions or distributions. The result is a matrix of size
% (n_a + n_e) × (N_a * N_e), where each column corresponds to one grid point in the joint
% state space.
%
% INPUTS
% ss [struct] : Steady-state structure containing grids and state/shock names
% sizes [struct] : Structure containing sizes of endogenous states and exogenous shocks
% field [char] : Field name of `ss` ('pol' or 'd'), indicating whether to construct
% the grid for policy functions or for the distribution
%
% OUTPUT
% sm [matrix] : Matrix of size (n_a + n_e) × (N_a * N_e), where each column contains
% one combination of state and shock values, ordered lexicographically.
%
% Notes:
% - The first `n_a` rows correspond to endogenous states; the next `n_e` to exogenous shocks.
% - Ordering is compatible with Kronecker-product-based interpolation logic.
% - Used in residual evaluation, interpolation, and expectation computation routines.
function sm = set_state_matrix(ss, sizes, field)
% Useful size vectors and variables
f = ss.(field);
N = sizes.(field).N_a*sizes.N_e;
% Setting the states matrix
sm = zeros(N, sizes.n_a+sizes.n_e);
% Filling the state matrix
prod_a_1_to_jm1 = 1;
prod_a_jp1_to_n_a = N;
for j=1:sizes.n_a
state = f.states{j};
prod_a_jp1_to_n_a = prod_a_jp1_to_n_a/sizes.(field).states.(state);
sm(:,j) = kron(ones(prod_a_1_to_jm1,1), kron(f.grids.(state), ones(prod_a_jp1_to_n_a,1)));
prod_a_1_to_jm1 = prod_a_1_to_jm1*sizes.(field).states.(state);
end
prod_theta_1_to_jm1 = sizes.(field).N_a;
prod_theta_jp1_to_n_theta = sizes.N_e;
for j=1:sizes.n_e
shock = f.shocks{j};
prod_theta_jp1_to_n_theta = prod_theta_jp1_to_n_theta/sizes.shocks.(shock);
sm(:,sizes.n_a+j) = kron(ones(prod_theta_1_to_jm1,1), kron(ss.shocks.grids.(shock), ones(prod_theta_jp1_to_n_theta,1)));
prod_theta_1_to_jm1 = prod_theta_1_to_jm1*sizes.shocks.(shock);
end
sm = sm';
end
\ No newline at end of file
matlab/+identification/analysis.m
View file @ a29fb27b
...@@ -24,7 +24,7 @@ function [ide_moments, ide_spectrum, ide_minimal, ide_hess, ide_reducedform, ide ...@@ -24,7 +24,7 @@ function [ide_moments, ide_spectrum, ide_minimal, ide_hess, ide_reducedform, ide
% * indpstderr [stderrparam_nbr by 1] % * indpstderr [stderrparam_nbr by 1]
% index of stderr parameters for which identification is checked % index of stderr parameters for which identification is checked
% * indpcorr [corrparam_nbr by 2] % * indpcorr [corrparam_nbr by 2]
% matrix of corr parmeters for which identification is checked % matrix of corr parameters for which identification is checked
% * options_ident [structure] % * options_ident [structure]
% identification options % identification options
% * dataset_info [structure] % * dataset_info [structure]
...@@ -70,7 +70,6 @@ function [ide_moments, ide_spectrum, ide_minimal, ide_hess, ide_reducedform, ide ...@@ -70,7 +70,6 @@ function [ide_moments, ide_spectrum, ide_minimal, ide_hess, ide_reducedform, ide
% * identification.checks_via_subsets % * identification.checks_via_subsets
% * isoctave % * isoctave
% * identification.get_jacobians (previously getJJ) % * identification.get_jacobians (previously getJJ)
% * matlab_ver_less_than
% * prior_bounds % * prior_bounds
% * resol % * resol
% * set_all_parameters % * set_all_parameters
...@@ -297,18 +296,27 @@ if info(1) == 0 %no errors in solution ...@@ -297,18 +296,27 @@ if info(1) == 0 %no errors in solution
options_.analytic_derivation = -2; %this sets asy_Hess=1 in dsge_likelihood.m options_.analytic_derivation = -2; %this sets asy_Hess=1 in dsge_likelihood.m
[info, oo_, options_, M_] = stoch_simul(M_, options_, oo_, options_.varobs); [info, oo_, options_, M_] = stoch_simul(M_, options_, oo_, options_.varobs);
dataset_ = dseries(oo_.endo_simul(options_.varobs_id,100+1:end)',dates('1Q1'), options_.varobs); %get information on moments dataset_ = dseries(oo_.endo_simul(options_.varobs_id,100+1:end)',dates('1Q1'), options_.varobs); %get information on moments
% set info on missing data
if dataset_info.missing.state
[dataset_info.missing.aindex, dataset_info.missing.number_of_observations, dataset_info.missing.no_more_missing_observations, dataset_info.missing.vindex] = ...
describe_missing_data(dataset_.data);
else
dataset_info.missing.aindex = num2cell(transpose(repmat(1:dataset_.vobs,dataset_.nobs,1)),1);
dataset_info.missing.no_more_missing_observations = 1;
end
derivatives_info.no_DLIK = 1; derivatives_info.no_DLIK = 1;
bounds = prior_bounds(bayestopt_, options_.prior_trunc); %reset bounds as lb and ub must only be operational during mode-finding bounds = prior_bounds(bayestopt_, options_.prior_trunc); %reset bounds as lb and ub must only be operational during mode-finding
%note that for order>1 we do not provide any information on DT,DYss,DOM in derivatives_info, such that dsge_likelihood creates an error. Therefore the computation will be based on simulated_moment_uncertainty for order>1. %note that for order>1 we do not provide any information on DT,DYss,DOM in derivatives_info, such that dsge_likelihood creates an error. Therefore the computation will be based on simulated_moment_uncertainty for order>1.
[~, info, ~, ~, AHess, ~, ~, M_, options_, ~, oo_.dr] = dsge_likelihood(params', dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, bounds, oo_.dr, oo_.steady_state,oo_.exo_steady_state, oo_.exo_det_steady_state, derivatives_info); %non-used output variables need to be set for octave for some reason [~, info, ~, ~, AHess] = dsge_likelihood(params', dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, bounds, oo_.dr, oo_.steady_state,oo_.exo_steady_state, oo_.exo_det_steady_state, derivatives_info);
%note that for the order of parameters in AHess we have: stderr parameters come first, corr parameters second, model parameters third. the order within these blocks corresponds to the order specified in the estimated_params block %note that for the order of parameters in AHess we have: stderr parameters come first, corr parameters second, model parameters third. the order within these blocks corresponds to the order specified in the estimated_params block
options_.analytic_derivation = analytic_derivation; %reset option options_.analytic_derivation = analytic_derivation; %reset option
AHess = -AHess; %take negative of hessian AHess = -AHess; %take negative of Hessian
if min(eig(AHess))<-tol_rank if min(eig(AHess))<-tol_rank
error('identification.analysis: Analytic Hessian is not positive semi-definite!') error('identification.analysis: Analytic Hessian is not positive semi-definite!')
end end
ide_hess.AHess = AHess; %store asymptotic Hessian ide_hess.AHess = AHess; %store asymptotic Hessian
%normalize asymptotic hessian %normalize asymptotic Hessian
deltaM = sqrt(diag(AHess)); deltaM = sqrt(diag(AHess));
iflag = any((deltaM.*deltaM)==0); %check if all second-order derivatives wrt to a single parameter are nonzero iflag = any((deltaM.*deltaM)==0); %check if all second-order derivatives wrt to a single parameter are nonzero
tildaM = AHess./((deltaM)*(deltaM')); %this normalization is for numerical purposes tildaM = AHess./((deltaM)*(deltaM')); %this normalization is for numerical purposes
...@@ -344,8 +352,8 @@ if info(1) == 0 %no errors in solution ...@@ -344,8 +352,8 @@ if info(1) == 0 %no errors in solution
siTMP = si_dMOMENTS./repmat(sd,[1 totparam_nbr]); siTMP = si_dMOMENTS./repmat(sd,[1 totparam_nbr]);
MIM = (siTMP'*VV(:,id))*(DD(id,id)\(WW(:,id)'*siTMP)); MIM = (siTMP'*VV(:,id))*(DD(id,id)\(WW(:,id)'*siTMP));
clear siTMP; clear siTMP;
ide_hess.AHess = MIM; %store asymptotic hessian ide_hess.AHess = MIM; %store asymptotic Hessian
%normalize asymptotic hessian %normalize asymptotic Hessian
deltaM = sqrt(diag(MIM)); deltaM = sqrt(diag(MIM));
iflag = any((deltaM.*deltaM)==0); iflag = any((deltaM.*deltaM)==0);
tildaM = MIM./((deltaM)*(deltaM')); tildaM = MIM./((deltaM)*(deltaM'));
......
matlab/+identification/bruteforce.m
View file @ a29fb27b
...@@ -59,7 +59,7 @@ tittxt1=strrep(tittxt1, '.', ''); ...@@ -59,7 +59,7 @@ tittxt1=strrep(tittxt1, '.', '');
cosnJ = zeros(totparam_nbr,max_dim_cova_group); %initialize cosnJ = zeros(totparam_nbr,max_dim_cova_group); %initialize
pars{totparam_nbr,max_dim_cova_group}=[]; %initialize pars{totparam_nbr,max_dim_cova_group}=[]; %initialize
for ll = 1:max_dim_cova_group for ll = 1:max_dim_cova_group
h = dyn_waitbar(0,['Brute force collinearity for ' int2str(ll) ' parameters.']); h = waitbar.run(0,['Brute force collinearity for ' int2str(ll) ' parameters.']);
for ii = 1:totparam_nbr for ii = 1:totparam_nbr
tmp = find([1:totparam_nbr]~=ii); tmp = find([1:totparam_nbr]~=ii);
tmp2 = nchoosek(tmp,ll); %find all possible combinations, ind16 could speed this up tmp2 = nchoosek(tmp,ll); %find all possible combinations, ind16 could speed this up
...@@ -81,9 +81,9 @@ for ll = 1:max_dim_cova_group ...@@ -81,9 +81,9 @@ for ll = 1:max_dim_cova_group
else else
pars{ii,ll} = NaN(1,ll); pars{ii,ll} = NaN(1,ll);
end end
dyn_waitbar(ii/totparam_nbr,h) waitbar.run(ii/totparam_nbr,h)
end end
dyn_waitbar_close(h); waitbar.close(h);
if TeX if TeX
filename = [OutputDirectoryName '/' fname '_collin_patterns_',tittxt1,'_' int2str(ll) '.tex']; filename = [OutputDirectoryName '/' fname '_collin_patterns_',tittxt1,'_' int2str(ll) '.tex'];
fidTeX = fopen(filename,'w'); fidTeX = fopen(filename,'w');
......
matlab/+identification/checks.m
View file @ a29fb27b
...@@ -10,7 +10,7 @@ function [condX, rankX, ind0, indno, ixno, Mco, Pco, jweak, jweak_pair] = checks ...@@ -10,7 +10,7 @@ function [condX, rankX, ind0, indno, ixno, Mco, Pco, jweak, jweak_pair] = checks
% test_flag = 0: Sample information matrix (Ahess) % test_flag = 0: Sample information matrix (Ahess)
% test_flag = 1: Jacobian of Moments (J), reduced-form (dTAU) or dynamic model (dLRE) % test_flag = 1: Jacobian of Moments (J), reduced-form (dTAU) or dynamic model (dLRE)
% test_flag = 2: Jacobian of minimal system (D) % test_flag = 2: Jacobian of minimal system (D)
% test_flag = 3: Gram matrix (hessian or correlation type matrix) of spectrum (G) % test_flag = 3: Gram matrix (Hessian or correlation type matrix) of spectrum (G)
% ------------------------------------------------------------------------- % -------------------------------------------------------------------------
% OUTPUTS % OUTPUTS
% * cond [double] condition number of X % * cond [double] condition number of X
......
matlab/+identification/checks_via_subsets.m
View file @ a29fb27b
function [ide_dynamic, ide_reducedform, ide_moments, ide_spectrum, ide_minimal] = checks_via_subsets(ide_dynamic, ide_reducedform, ide_moments, ide_spectrum, ide_minimal, totparam_nbr, modparam_nbr, options_ident,error_indicator) function [ide_dynamic, ide_reducedform, ide_moments, ide_spectrum, ide_minimal] = checks_via_subsets(ide_dynamic, ide_reducedform, ide_moments, ide_spectrum, ide_minimal, totparam_nbr, modparam_nbr, options_ident,error_indicator)
%[ide_dynamic, ide_reducedform, ide_moments, ide_spectrum, ide_minimal] = checks_via_subsets(ide_dynamic, ide_reducedform, ide_moments, ide_spectrum, ide_minimal, totparam_nbr, modparam_nbr, options_ident,error_indicator) %[ide_dynamic, ide_reducedform, ide_moments, ide_spectrum, ide_minimal] = checks_via_subsets(ide_dynamic, ide_reducedform, ide_moments, ide_spectrum, ide_minimal, totparam_nbr, modparam_nbr, options_ident,error_indicator)
% ------------------------------------------------------------------------- % -------------------------------------------------------------------------
% Finds problematic sets of paramters via checking the necessary rank condition % Finds problematic sets of parameters via checking the necessary rank condition
% of the Jacobians for all possible combinations of parameters. The rank is % of the Jacobians for all possible combinations of parameters. The rank is
% computed via an inbuild function based on the SVD, similar to matlab's % computed via an inbuilt function based on the SVD, similar to MATLAB's
% rank. The idea is that once we have the Jacobian for all parameters, we % rank. The idea is that once we have the Jacobian for all parameters, we
% can easily set up Jacobians containing all combinations of parameters by % can easily set up Jacobians containing all combinations of parameters by
% picking the relevant columns/elements of the full Jacobian. Then the rank % picking the relevant columns/elements of the full Jacobian. Then the rank
% of these smaller Jacobians indicates whether this paramter combination is % of these smaller Jacobians indicates whether this parameter combination is
% identified or not. To speed up computations: % identified or not. To speed up computations:
% (1) single parameters are removed from possible higher-order sets, % (1) single parameters are removed from possible higher-order sets,
% (2) for parameters that are collinear, i.e. rank failure for 2 element sets, % (2) for parameters that are collinear, i.e. rank failure for 2 element sets,
% we replace the second parameter by the first one, and then compute % we replace the second parameter by the first one, and then compute
% higher-order combinations [uncommented] % higher-order combinations [uncommented]
% (3) all lower-order problematic sets are removed from higher-order sets % (3) all lower-order problematic sets are removed from higher-order sets
% by an inbuild function % by an inbuilt function
% (4) we could replace nchoosek by a mex version, e.g. VChooseK % (4) we could replace nchoosek by a mex version, e.g. VChooseK
% (https://de.mathworks.com/matlabcentral/fileexchange/26190-vchoosek) as % (https://de.mathworks.com/matlabcentral/fileexchange/26190-vchoosek) as
% nchoosek could be the bottleneck in terms of speed (and memory for models % nchoosek could be the bottleneck in terms of speed (and memory for models
...@@ -322,7 +322,7 @@ end ...@@ -322,7 +322,7 @@ end
indtotparam = unique([indparam_dDYNAMIC indparam_dREDUCEDFORM indparam_dMOMENTS indparam_dSPECTRUM indparam_dMINIMAL]); indtotparam = unique([indparam_dDYNAMIC indparam_dREDUCEDFORM indparam_dMOMENTS indparam_dSPECTRUM indparam_dMINIMAL]);
for j=2:min(length(indtotparam),max_dim_subsets_groups) % Check j-element subsets for j=2:min(length(indtotparam),max_dim_subsets_groups) % Check j-element subsets
h = dyn_waitbar(0,['Brute force collinearity for ' int2str(j) ' parameters.']); h = waitbar.run(0,['Brute force collinearity for ' int2str(j) ' parameters.']);
%Step1: get all possible unique subsets of j elements %Step1: get all possible unique subsets of j elements
if ~no_identification_dynamic ... if ~no_identification_dynamic ...
|| (~no_identification_reducedform && ~error_indicator.identification_reducedform)... || (~no_identification_reducedform && ~error_indicator.identification_reducedform)...
...@@ -424,7 +424,7 @@ for j=2:min(length(indtotparam),max_dim_subsets_groups) % Check j-element subset ...@@ -424,7 +424,7 @@ for j=2:min(length(indtotparam),max_dim_subsets_groups) % Check j-element subset
end end
end end
end end
dyn_waitbar(k/maxk,h) waitbar.run(k/maxk,h)
end end
%Step 4: Compare rank conditions for all possible subsets. If rank condition is violated, then the corresponding numbers of the parameters are stored %Step 4: Compare rank conditions for all possible subsets. If rank condition is violated, then the corresponding numbers of the parameters are stored
...@@ -473,7 +473,7 @@ for j=2:min(length(indtotparam),max_dim_subsets_groups) % Check j-element subset ...@@ -473,7 +473,7 @@ for j=2:min(length(indtotparam),max_dim_subsets_groups) % Check j-element subset
% indtotparam(ismember(indtotparam,idx2(:,2))) = []; % indtotparam(ismember(indtotparam,idx2(:,2))) = [];
% end % end
% end % end
dyn_waitbar_close(h); waitbar.close(h);
end end
%% Save output variables %% Save output variables
......
matlab/+identification/get_perturbation_params_derivs.m
View file @ a29fb27b
function DERIVS = get_perturbation_params_derivs(M_, options_, estim_params_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state, indpmodel, indpstderr, indpcorr, d2flag) function DERIVS = get_perturbation_params_derivs(M_, options_, estim_params_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state, indpmodel, indpstderr, indpcorr, d2flag)
% DERIVS = get_perturbation_params_derivs(M_, options_, estim_params_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state, indpmodel, indpstderr, indpcorr, d2flag) % DERIVS = get_perturbation_params_derivs(M_, options_, estim_params_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state, indpmodel, indpstderr, indpcorr, d2flag)
% previously getH.m in dynare 4.5 % previously getH.m in Dynare 4.5
% ------------------------------------------------------------------------- % -------------------------------------------------------------------------
% Computes derivatives (with respect to parameters) of % Computes derivatives (with respect to parameters) of
% (1) steady-state (ys) and covariance matrix of shocks (Sigma_e) % (1) steady-state (ys) and covariance matrix of shocks (Sigma_e)
% (2) dynamic model jacobians (g1, g2, g3) % (2) dynamic model Jacobians (g1, g2, g3)
% (3) perturbation solution matrices: % (3) perturbation solution matrices:
% * order==1: ghx,ghu % * order==1: ghx,ghu
% * order==2: ghx,ghu,ghxx,ghxu,ghuu,ghs2 % * order==2: ghx,ghu,ghxx,ghxu,ghuu,ghs2
...@@ -40,7 +40,7 @@ function DERIVS = get_perturbation_params_derivs(M_, options_, estim_params_, dr ...@@ -40,7 +40,7 @@ function DERIVS = get_perturbation_params_derivs(M_, options_, estim_params_, dr
% dYss: [endo_nbr by modparam_nbr] in DR order % dYss: [endo_nbr by modparam_nbr] in DR order
% Jacobian (wrt model parameters only) of steady state, i.e. ys(order_var,:) % Jacobian (wrt model parameters only) of steady state, i.e. ys(order_var,:)
% dSigma_e: [exo_nbr by exo_nbr by totparam_nbr] in declaration order % dSigma_e: [exo_nbr by exo_nbr by totparam_nbr] in declaration order
% Jacobian (wrt to all paramters) of covariance matrix of shocks, i.e. Sigma_e % Jacobian (wrt to all parameters) of covariance matrix of shocks, i.e. Sigma_e
% dg1: [endo_nbr by yy0ex0_nbr by modparam_nbr] in DR order % dg1: [endo_nbr by yy0ex0_nbr by modparam_nbr] in DR order
% Parameter Jacobian of first derivative (wrt dynamic model variables) of dynamic model (wrt to model parameters only) % Parameter Jacobian of first derivative (wrt dynamic model variables) of dynamic model (wrt to model parameters only)
% dg2: [endo_nbr by yy0ex0_nbr^2*modparam_nbr] in DR order % dg2: [endo_nbr by yy0ex0_nbr^2*modparam_nbr] in DR order
...@@ -56,7 +56,7 @@ function DERIVS = get_perturbation_params_derivs(M_, options_, estim_params_, dr ...@@ -56,7 +56,7 @@ function DERIVS = get_perturbation_params_derivs(M_, options_, estim_params_, dr
% dghu: [endo_nbr by exo_nbr by totparam_nbr] in DR order % dghu: [endo_nbr by exo_nbr by totparam_nbr] in DR order
% Jacobian (wrt to all parameters) of first-order perturbation solution matrix ghu % Jacobian (wrt to all parameters) of first-order perturbation solution matrix ghu
% dOm: [endo_nbr by endo_nbr by totparam_nbr] in DR order % dOm: [endo_nbr by endo_nbr by totparam_nbr] in DR order
% Jacobian (wrt to all paramters) of Om = ghu*Sigma_e*transpose(ghu) % Jacobian (wrt to all parameters) of Om = ghu*Sigma_e*transpose(ghu)
% dghxx [endo_nbr by nspred*nspred by totparam_nbr] in DR order % dghxx [endo_nbr by nspred*nspred by totparam_nbr] in DR order
% Jacobian (wrt to all parameters) of second-order perturbation solution matrix ghxx % Jacobian (wrt to all parameters) of second-order perturbation solution matrix ghxx
% dghxu [endo_nbr by nspred*exo_nbr by totparam_nbr] in DR order % dghxu [endo_nbr by nspred*exo_nbr by totparam_nbr] in DR order
...@@ -295,7 +295,7 @@ if analytic_derivation_mode == -1 ...@@ -295,7 +295,7 @@ if analytic_derivation_mode == -1
% - covariance matrix of shocks: dSigma_e % - covariance matrix of shocks: dSigma_e
% - perturbation solution matrices: dghx, dghu, dghxx, dghxu, dghuu, dghs2, dghxxx, dghxxu, dghxuu, dghuuu, dghxss, dghuss % - perturbation solution matrices: dghx, dghu, dghxx, dghxu, dghuu, dghs2, dghxxx, dghxxu, dghxuu, dghuuu, dghxss, dghuss
%Parameter Jacobian of covariance matrix and solution matrices (wrt selected stderr, corr and model paramters) %Parameter Jacobian of covariance matrix and solution matrices (wrt selected stderr, corr and model parameters)
dSig_gh = identification.fjaco(numerical_objective_fname, xparam1, 'perturbation_solution', estim_params_, M_, options_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state); dSig_gh = identification.fjaco(numerical_objective_fname, xparam1, 'perturbation_solution', estim_params_, M_, options_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state);
ind_Sigma_e = (1:exo_nbr^2); ind_Sigma_e = (1:exo_nbr^2);
ind_ghx = ind_Sigma_e(end) + (1:endo_nbr*nspred); ind_ghx = ind_Sigma_e(end) + (1:endo_nbr*nspred);
...@@ -327,7 +327,7 @@ if analytic_derivation_mode == -1 ...@@ -327,7 +327,7 @@ if analytic_derivation_mode == -1
DERIVS.dghxss = reshape(dSig_gh(ind_ghxss,:), [endo_nbr nspred totparam_nbr]); %in tensor notation, wrt selected parameters DERIVS.dghxss = reshape(dSig_gh(ind_ghxss,:), [endo_nbr nspred totparam_nbr]); %in tensor notation, wrt selected parameters
DERIVS.dghuss = reshape(dSig_gh(ind_ghuss,:), [endo_nbr exo_nbr totparam_nbr]); %in tensor notation, wrt selected parameters DERIVS.dghuss = reshape(dSig_gh(ind_ghuss,:), [endo_nbr exo_nbr totparam_nbr]); %in tensor notation, wrt selected parameters
end end
% Parameter Jacobian of Om=ghu*Sigma_e*ghu' and Correlation_matrix (wrt selected stderr, corr and model paramters) % Parameter Jacobian of Om=ghu*Sigma_e*ghu' and Correlation_matrix (wrt selected stderr, corr and model parameters)
DERIVS.dOm = zeros(endo_nbr,endo_nbr,totparam_nbr); %initialize in tensor notation DERIVS.dOm = zeros(endo_nbr,endo_nbr,totparam_nbr); %initialize in tensor notation
DERIVS.dCorrelation_matrix = zeros(exo_nbr,exo_nbr,totparam_nbr); %initialize in tensor notation DERIVS.dCorrelation_matrix = zeros(exo_nbr,exo_nbr,totparam_nbr); %initialize in tensor notation
if ~isempty(indpstderr) %derivatives of ghu wrt stderr parameters are zero by construction if ~isempty(indpstderr) %derivatives of ghu wrt stderr parameters are zero by construction
...@@ -372,7 +372,7 @@ if analytic_derivation_mode == -1 ...@@ -372,7 +372,7 @@ if analytic_derivation_mode == -1
end end
if d2flag if d2flag
% Hessian (wrt paramters) of steady state and first-order solution matrices ghx and Om % Hessian (wrt parameters) of steady state and first-order solution matrices ghx and Om
% note that hessian_sparse.m (contrary to hessian.m) does not take symmetry into account, but focuses already on unique values % note that hessian_sparse.m (contrary to hessian.m) does not take symmetry into account, but focuses already on unique values
options_.order = 1; %make sure only first order options_.order = 1; %make sure only first order
d2Yss_KalmanA_Om = identification.hessian_sparse(numerical_objective_fname, xparam1, gstep, 'Kalman_Transition', estim_params_, M_, options_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state); d2Yss_KalmanA_Om = identification.hessian_sparse(numerical_objective_fname, xparam1, gstep, 'Kalman_Transition', estim_params_, M_, options_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state);
...@@ -1530,7 +1530,7 @@ function g22 = get_2nd_deriv_mat(gpp,i,j,npar) ...@@ -1530,7 +1530,7 @@ function g22 = get_2nd_deriv_mat(gpp,i,j,npar)
% output: % output:
% g22: [npar by npar] Hessian matrix (wrt parameters) of Jacobian of dynamic model for equation i % g22: [npar by npar] Hessian matrix (wrt parameters) of Jacobian of dynamic model for equation i
% rows: first parameter in Hessian % rows: first parameter in Hessian
% columns: second paramater in Hessian % columns: second parameter in Hessian
g22=zeros(npar,npar); g22=zeros(npar,npar);
is=find(gpp(:,1)==i); is=find(gpp(:,1)==i);
...@@ -1543,7 +1543,7 @@ return ...@@ -1543,7 +1543,7 @@ return
function r22 = get_all_resid_2nd_derivs(rpp,m,npar) function r22 = get_all_resid_2nd_derivs(rpp,m,npar)
% inputs: % inputs:
% - rpp: [#second_order_residual_terms by 4] double Hessian matrix (wrt paramters) of model equations % - rpp: [#second_order_residual_terms by 4] double Hessian matrix (wrt parameters) of model equations
% rows: respective derivative term % rows: respective derivative term
% 1st column: equation number of the term appearing % 1st column: equation number of the term appearing
% 2nd column: number of the first parameter in derivative % 2nd column: number of the first parameter in derivative
...@@ -1570,7 +1570,7 @@ return ...@@ -1570,7 +1570,7 @@ return
function h2 = get_hess_deriv(hp,i,j,m,npar) function h2 = get_hess_deriv(hp,i,j,m,npar)
% inputs: % inputs:
% - hp: [#first_order_Hessian_terms by 5] double Jacobian matrix (wrt paramters) of dynamic Hessian % - hp: [#first_order_Hessian_terms by 5] double Jacobian matrix (wrt parameters) of dynamic Hessian
% rows: respective derivative term % rows: respective derivative term
% 1st column: equation number of the term appearing % 1st column: equation number of the term appearing
% 2nd column: column number of first variable in Hessian of the dynamic model % 2nd column: column number of first variable in Hessian of the dynamic model
......
matlab/+identification/plot.m
View file @ a29fb27b
...@@ -25,7 +25,7 @@ function plot(M_, params, idemoments, idehess, idemodel, idelre, advanced, tittx ...@@ -25,7 +25,7 @@ function plot(M_, params, idemoments, idehess, idemodel, idelre, advanced, tittx
% SPECIAL REQUIREMENTS % SPECIAL REQUIREMENTS
% None % None
% Copyright © 2008-2023 Dynare Team % Copyright © 2008-2025 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
...@@ -62,8 +62,7 @@ if SampleSize == 1 ...@@ -62,8 +62,7 @@ if SampleSize == 1
subplot(211) subplot(211)
mmm = (idehess.ide_strength_dMOMENTS); mmm = (idehess.ide_strength_dMOMENTS);
[~, is] = sort(mmm); [~, is] = sort(mmm);
if ~all(isnan(idehess.ide_strength_dMOMENTS_prior)) ... if ~all(isnan(idehess.ide_strength_dMOMENTS_prior))
&& ~(nparam == 1 && ~isoctave && matlab_ver_less_than('9.7')) % MATLAB < R2019b does not accept bar(1, [2 3])
bar(1:nparam,log([idehess.ide_strength_dMOMENTS(:,is)' idehess.ide_strength_dMOMENTS_prior(:,is)'])) bar(1:nparam,log([idehess.ide_strength_dMOMENTS(:,is)' idehess.ide_strength_dMOMENTS_prior(:,is)']))
else else
bar(1:nparam,log(idehess.ide_strength_dMOMENTS(:,is)')) bar(1:nparam,log(idehess.ide_strength_dMOMENTS(:,is)'))
...@@ -113,8 +112,7 @@ if SampleSize == 1 ...@@ -113,8 +112,7 @@ if SampleSize == 1
end end
subplot(212) subplot(212)
if ~all(isnan(idehess.deltaM_prior)) ... if ~all(isnan(idehess.deltaM_prior))
&& ~(nparam == 1 && ~isoctave && matlab_ver_less_than('9.7')) % MATLAB < R2019b does not accept bar(1, [2 3])
bar(1:nparam, log([idehess.deltaM(is) idehess.deltaM_prior(is)])) bar(1:nparam, log([idehess.deltaM(is) idehess.deltaM_prior(is)]))
else else
bar(1:nparam, log([idehess.deltaM(is)])) bar(1:nparam, log([idehess.deltaM(is)]))
......
matlab/+identification/run.m
View file @ a29fb27b
...@@ -11,9 +11,9 @@ function [pdraws, STO_REDUCEDFORM, STO_MOMENTS, STO_DYNAMIC, STO_si_dDYNAMIC, ST ...@@ -11,9 +11,9 @@ function [pdraws, STO_REDUCEDFORM, STO_MOMENTS, STO_DYNAMIC, STO_si_dDYNAMIC, ST
% to put identification in your mod file, otherwise the preprocessor won't provide all necessary objects % to put identification in your mod file, otherwise the preprocessor won't provide all necessary objects
% ========================================================================= % =========================================================================
% INPUTS % INPUTS
% * M_ [structure] Matlab's structure describing the model % * M_ [structure] MATLAB's structure describing the model
% * oo_ [structure] Matlab's structure describing the results % * oo_ [structure] MATLAB's structure describing the results
% * options_ [structure] Matlab's structure describing the current options % * options_ [structure] MATLAB's structure describing the current options
% * bayestopt_ [structure] describing the priors % * bayestopt_ [structure] describing the priors
% * estim_params_ [structure] characterizing parameters to be estimated % * estim_params_ [structure] characterizing parameters to be estimated
% * options_ident [structure] identification options % * options_ident [structure] identification options
...@@ -37,8 +37,8 @@ function [pdraws, STO_REDUCEDFORM, STO_MOMENTS, STO_DYNAMIC, STO_si_dDYNAMIC, ST ...@@ -37,8 +37,8 @@ function [pdraws, STO_REDUCEDFORM, STO_MOMENTS, STO_DYNAMIC, STO_si_dDYNAMIC, ST
% This function calls % This function calls
% * checkpath % * checkpath
% * identification.display % * identification.display
% * dyn_waitbar % * waitbar.run
% * dyn_waitbar_close % * waitbar.close
% * get_all_parameters % * get_all_parameters
% * get_posterior_parameters % * get_posterior_parameters
% * get_the_name % * get_the_name
...@@ -210,7 +210,7 @@ end ...@@ -210,7 +210,7 @@ end
% 5: i) option 2 for non-stationary elements by setting their initial variance in the forecast error matrix to 10 on the diagonal and all co-variances to 0 and % 5: i) option 2 for non-stationary elements by setting their initial variance in the forecast error matrix to 10 on the diagonal and all co-variances to 0 and
% ii) option 1 for the stationary elements % ii) option 1 for the stationary elements
options_ident = set_default_option(options_ident,'analytic_derivation',1); options_ident = set_default_option(options_ident,'analytic_derivation',1);
% 1: analytic derivation of gradient and hessian of likelihood in dsge_likelihood.m, only works for stationary models, i.e. kalman_algo<3 % 1: analytic derivation of gradient and Hessian of likelihood in dsge_likelihood.m, only works for stationary models, i.e. kalman_algo<3
options_ident = set_default_option(options_ident,'order',1); options_ident = set_default_option(options_ident,'order',1);
% 1: first-order perturbation approximation, identification is based on linear state space system % 1: first-order perturbation approximation, identification is based on linear state space system
% 2: second-order perturbation approximation, identification is based on second-order pruned state space system % 2: second-order perturbation approximation, identification is based on second-order pruned state space system
...@@ -240,7 +240,7 @@ end ...@@ -240,7 +240,7 @@ end
% check for external draws, i.e. set pdraws0 for a gsa analysis % check for external draws, i.e. set pdraws0 for a gsa analysis
if options_ident.gsa_sample_file if options_ident.gsa_sample_file
GSAFolder = checkpath('gsa',dname); GSAFolder = CheckPath('gsa',dname);
if options_ident.gsa_sample_file==1 if options_ident.gsa_sample_file==1
load([GSAFolder,filesep,fname,'_prior'],'lpmat','lpmat0','istable'); load([GSAFolder,filesep,fname,'_prior'],'lpmat','lpmat0','istable');
elseif options_ident.gsa_sample_file==2 elseif options_ident.gsa_sample_file==2
...@@ -297,9 +297,14 @@ options_.prior_mc = options_ident.prior_mc; ...@@ -297,9 +297,14 @@ options_.prior_mc = options_ident.prior_mc;
options_.schur_vec_tol = options_ident.schur_vec_tol; options_.schur_vec_tol = options_ident.schur_vec_tol;
options_.nomoments = 0; options_.nomoments = 0;
options_.analytic_derivation=options_ident.analytic_derivation; options_.analytic_derivation=options_ident.analytic_derivation;
% 1: analytic derivation of gradient and hessian of likelihood in dsge_likelihood.m, only works for stationary models, i.e. kalman_algo<3 % 1: analytic derivation of gradient and Hessian of likelihood in dsge_likelihood.m, only works for stationary models, i.e. kalman_algo<3
options_ = set_default_option(options_,'datafile',''); options_ = set_default_option(options_,'datafile','');
options_.mode_compute = 0; options_.mode_compute = 0;
if strcmp('slice',options_.posterior_sampler_options.posterior_sampling_method)
if strfind(options_.posterior_sampler_options.sampling_opt,'slice_initialize_with_mode'',1')
options_.posterior_sampler_options.sampling_opt=strrep(options_.posterior_sampler_options.sampling_opt,'slice_initialize_with_mode'',1','slice_initialize_with_mode'',0'); % reset option to prevent crash in check_posterior_sampler_options.m if mode_compute is set to 0, see https://git.dynare.org/Dynare/dynare/-/issues/1957
end
end
options_.plot_priors = 0; options_.plot_priors = 0;
options_.smoother = 1; options_.smoother = 1;
options_.options_ident = []; options_.options_ident = [];
...@@ -349,7 +354,7 @@ if prior_exist % use estimated_params block ...@@ -349,7 +354,7 @@ if prior_exist % use estimated_params block
if ~isempty(estim_params_.var_exo) if ~isempty(estim_params_.var_exo)
indpstderr = estim_params_.var_exo(:,1); %values correspond to varexo declaration order, row number corresponds to order in estimated_params indpstderr = estim_params_.var_exo(:,1); %values correspond to varexo declaration order, row number corresponds to order in estimated_params
end end
indpcorr=[]; %initialize matrix for corr paramters indpcorr=[]; %initialize matrix for corr parameters
if ~isempty(estim_params_.corrx) if ~isempty(estim_params_.corrx)
indpcorr = estim_params_.corrx(:,1:2); %values correspond to varexo declaration order, row number corresponds to order in estimated_params indpcorr = estim_params_.corrx(:,1:2); %values correspond to varexo declaration order, row number corresponds to order in estimated_params
end end
...@@ -485,7 +490,7 @@ if iload <=0 ...@@ -485,7 +490,7 @@ if iload <=0
kk=0; kk=0;
while kk<50 && info(1) while kk<50 && info(1)
kk=kk+1; kk=kk+1;
params = Prior.draw(); params = Prior.draw()'; %column vector is expected
options_ident.tittxt = 'Random_prior_params'; %title text for graphs and figures options_ident.tittxt = 'Random_prior_params'; %title text for graphs and figures
% perform identification analysis % perform identification analysis
[ide_moments_point, ide_spectrum_point, ide_minimal_point, ide_hess_point, ide_reducedform_point, ide_dynamic_point, ~, info, error_indicator_point] = ... [ide_moments_point, ide_spectrum_point, ide_minimal_point, ide_hess_point, ide_reducedform_point, ide_dynamic_point, ~, info, error_indicator_point] = ...
...@@ -504,7 +509,10 @@ if iload <=0 ...@@ -504,7 +509,10 @@ if iload <=0
else else
% found a (random) point that solves the model % found a (random) point that solves the model
fprintf('Found a random draw from the priors that solves the model:\n'); fprintf('Found a random draw from the priors that solves the model:\n');
disp(params); labels = name;
headers = {'Name', 'Value'};
lh = cellofchararraymaxlength(labels)+2;
dyntable(options_, 'Feasible draw', headers, labels, params', lh, 10, 6);
fprintf('Identification now continues for this draw.'); fprintf('Identification now continues for this draw.');
parameters = 'Random_prior_params'; parameters = 'Random_prior_params';
parameters_TeX = 'Random prior parameter draw'; parameters_TeX = 'Random prior parameter draw';
...@@ -525,7 +533,7 @@ if iload <=0 ...@@ -525,7 +533,7 @@ if iload <=0
if SampleSize > 1 if SampleSize > 1
% initializations for Monte Carlo Analysis % initializations for Monte Carlo Analysis
fprintf('\nMonte Carlo Testing\n'); fprintf('\nMonte Carlo Testing\n');
h = dyn_waitbar(0,'Monte Carlo identification checks ...'); h = waitbar.run(0,'Monte Carlo identification checks ...');
iteration = 0; % initialize counter for admissable draws iteration = 0; % initialize counter for admissable draws
run_index = 0; % initialize counter for admissable draws after saving previous draws to file(s) run_index = 0; % initialize counter for admissable draws after saving previous draws to file(s)
file_index = 0; % initialize counter for files (if options_.MaxNumberOfBytes is reached, we store results in files) file_index = 0; % initialize counter for files (if options_.MaxNumberOfBytes is reached, we store results in files)
...@@ -739,13 +747,13 @@ if iload <=0 ...@@ -739,13 +747,13 @@ if iload <=0
run_index = 0; % reset index run_index = 0; % reset index
end end
if SampleSize > 1 && mod(iteration,3) if SampleSize > 1 && mod(iteration,3)
dyn_waitbar(iteration/SampleSize, h, ['MC identification checks ', int2str(iteration), '/', int2str(SampleSize)]); waitbar.run(iteration/SampleSize, h, ['MC identification checks ', int2str(iteration), '/', int2str(SampleSize)]);
end end
end end
end end
if SampleSize > 1 if SampleSize > 1
dyn_waitbar_close(h); waitbar.close(h);
normalize_STO_DYNAMIC = std(STO_DYNAMIC,0,2); normalize_STO_DYNAMIC = std(STO_DYNAMIC,0,2);
if ~options_MC.no_identification_reducedform if ~options_MC.no_identification_reducedform
normalize_STO_REDUCEDFORM = std(STO_REDUCEDFORM,0,2); normalize_STO_REDUCEDFORM = std(STO_REDUCEDFORM,0,2);
......