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 890 additions and 358 deletions
matlab/+hank/private/check_fun.m 0 → 100644
View file @ 4a32acc2
% 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 @ 4a32acc2
% 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 @ 4a32acc2
% 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/convertAimCodeToInfo.m → matlab/+hank/private/check_missingredundant.m
View file @ 4a32acc2
function [info] = convertAimCodeToInfo(aimCode) % Copyright © 2025 Dynare Team
% function [info] = convertAimCodeToInfo(aimCode)
% Returns an appropriate code for print_info
%
% INPUTS
% aimCode [integer] code returned by AIM
% (aimCode==1) e='Aim: unique solution.';
% (aimCode==2) e='Aim: roots not correctly computed by real_schur.';
% (aimCode==3) e='Aim: too many big roots.';
% (aimCode==35) e='Aim: too many big roots, and q(:,right) is singular.';
% (aimCode==4) e='Aim: too few big roots.';
% (aimCode==45) e='Aim: too few big roots, and q(:,right) is singular.';
% (aimCode==5) e='Aim: q(:,right) is singular.';
% (aimCode==61) e='Aim: too many exact shiftrights.';
% (aimCode==62) e='Aim: too many numeric shiftrights.';
% (aimCode==63) e='Aim: A is NAN or INF.';
% (aimCode==64) e='Aim: Problem in SPEIG.';
%
% OUTPUTS
% info [integer] Code to be used to print error in print_info.m
% Copyright © 2011-2017 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
...@@ -35,30 +14,34 @@ function [info] = convertAimCodeToInfo(aimCode) ...@@ -35,30 +14,34 @@ function [info] = convertAimCodeToInfo(aimCode)
% %
% 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/>.
%
switch aimCode % Check for missing and redundant fields in a structure.
case 1 %
info = 0; % no problem encountered % INPUTS
case 2 % - f [structure]: structure to check
info = 102; % - f_name [char]: name to display in error/warning messages
case 3 % - symbs [cell array of char]: list of expected field names
info = 103; % - nowarningredundant [logical]: if true, suppresses warnings about redundant fields
case 35 %
info = 135; % OUTPUTS
case 4 % - (none) This function throws an error if any expected field is missing,
info = 104; % and issues a warning if redundant fields are found (unless nowarningredundant is true).
case 45 %
info = 145; % DESCRIPTION
case 5 % Verifies that all fields listed in `symbs` are present in the structure `f`.
info = 105; % Throws an error if any expected field is missing.
case 61 % If `nowarningredundant` is false, checks if `f` contains fields not listed
info = 161; % in `symbs`, and issues a warning if redundant fields are found.
case 62 function check_missingredundant(f, f_name, symbs, nowarningredundant)
info = 162; fields = fieldnames(f);
case 63 symbs_in_fields = ismember(symbs, fields);
info = 163; if ~all(symbs_in_fields)
case 64 error('Misspecified steady-state input `ss`. The following fields are missing in `%s`: %s.', f_name, strjoin(symbs(~symbs_in_fields)));
info = 164; end
otherwise if ~nowarningredundant
info = 1; 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 end
\ No newline at end of file
matlab/+hank/private/check_permutation.m 0 → 100644
View file @ 4a32acc2
% 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 @ 4a32acc2
% 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 @ 4a32acc2
% 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 @ 4a32acc2
% 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 @ 4a32acc2
% 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 → matlab/+identification/analysis.m
View file @ 4a32acc2
function [ide_moments, ide_spectrum, ide_minimal, ide_hess, ide_reducedform, ide_dynamic, derivatives_info, info, error_indicator] = identification_analysis(params, indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, init) function [ide_moments, ide_spectrum, ide_minimal, ide_hess, ide_reducedform, ide_dynamic, derivatives_info, info, error_indicator] = analysis(M_,options_,oo_,bayestopt_,estim_params_,params, indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, init)
% [ide_moments, ide_spectrum, ide_minimal, ide_hess, ide_reducedform, ide_dynamic, derivatives_info, info, error_indicator] = identification_analysis(params, indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, init) % [ide_moments, ide_spectrum, ide_minimal, ide_hess, ide_reducedform, ide_dynamic, derivatives_info, info, error_indicator] = analysis(M_,options_,oo_,bayestopt_,estim_params_,params, indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, init)
% ------------------------------------------------------------------------- % -------------------------------------------------------------------------
% This function wraps all identification analysis, i.e. it % This function wraps all identification analysis, i.e. it
% (1) wraps functions for the theoretical identification analysis based on moments (Iskrev, 2010), % (1) wraps functions for the theoretical identification analysis based on moments (Iskrev, 2010),
...@@ -12,6 +12,11 @@ function [ide_moments, ide_spectrum, ide_minimal, ide_hess, ide_reducedform, ide ...@@ -12,6 +12,11 @@ function [ide_moments, ide_spectrum, ide_minimal, ide_hess, ide_reducedform, ide
% moments, spectrum, reduced-form solution and dynamic model derivatives % moments, spectrum, reduced-form solution and dynamic model derivatives
% ========================================================================= % =========================================================================
% INPUTS % INPUTS
% * M_ [structure] describing the model
% * options_ [structure] describing the options
% * oo_ [structure] storing the results
% * bayestopt_ [structure] describing the priors
% * estim_params_ [structure] characterizing parameters to be estimated
% * params [mc_sample_nbr by totparam_nbr] % * params [mc_sample_nbr by totparam_nbr]
% parameter values for identification checks % parameter values for identification checks
% * indpmodel [modparam_nbr by 1] % * indpmodel [modparam_nbr by 1]
...@@ -19,7 +24,7 @@ function [ide_moments, ide_spectrum, ide_minimal, ide_hess, ide_reducedform, ide ...@@ -19,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]
...@@ -53,19 +58,18 @@ function [ide_moments, ide_spectrum, ide_minimal, ide_hess, ide_reducedform, ide ...@@ -53,19 +58,18 @@ function [ide_moments, ide_spectrum, ide_minimal, ide_hess, ide_reducedform, ide
% indicator on problems % indicator on problems
% ------------------------------------------------------------------------- % -------------------------------------------------------------------------
% This function is called by % This function is called by
% * dynare_identification.m % * identification.run
% ------------------------------------------------------------------------- % -------------------------------------------------------------------------
% This function calls % This function calls
% * [M_.fname,'.dynamic'] % * [M_.fname,'.dynamic']
% * dseries % * dseries
% * dsge_likelihood.m % * dsge_likelihood.m
% * dyn_vech % * dyn_vech
% * ident_bruteforce % * identification.bruteforce
% * identification_checks % * identification.checks
% * identification_checks_via_subsets % * identification.checks_via_subsets
% * isoctave % * isoctave
% * get_identification_jacobians (previously getJJ) % * identification.get_jacobians (previously getJJ)
% * matlab_ver_less_than
% * prior_bounds % * prior_bounds
% * resol % * resol
% * set_all_parameters % * set_all_parameters
...@@ -73,7 +77,7 @@ function [ide_moments, ide_spectrum, ide_minimal, ide_hess, ide_reducedform, ide ...@@ -73,7 +77,7 @@ function [ide_moments, ide_spectrum, ide_minimal, ide_hess, ide_reducedform, ide
% * stoch_simul % * stoch_simul
% * vec % * vec
% ========================================================================= % =========================================================================
% Copyright © 2008-2021 Dynare Team % Copyright © 2008-2023 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
...@@ -91,7 +95,6 @@ function [ide_moments, ide_spectrum, ide_minimal, ide_hess, ide_reducedform, ide ...@@ -91,7 +95,6 @@ function [ide_moments, ide_spectrum, ide_minimal, ide_hess, ide_reducedform, ide
% along with Dynare. If not, see <https://www.gnu.org/licenses/>. % along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% ========================================================================= % =========================================================================
global oo_ M_ options_ bayestopt_ estim_params_
persistent ind_dMOMENTS ind_dREDUCEDFORM ind_dDYNAMIC persistent ind_dMOMENTS ind_dREDUCEDFORM ind_dDYNAMIC
% persistent indices are necessary, because in a MC loop the numerical threshold % persistent indices are necessary, because in a MC loop the numerical threshold
% used may provide vectors of different length, leading to crashes in MC loops % used may provide vectors of different length, leading to crashes in MC loops
...@@ -116,7 +119,7 @@ if ~isempty(estim_params_) ...@@ -116,7 +119,7 @@ if ~isempty(estim_params_)
M_ = set_all_parameters(params,estim_params_,M_); M_ = set_all_parameters(params,estim_params_,M_);
end end
%get options (see dynare_identification.m for description of options) %get options (see identification.run.m for description of options)
nlags = options_ident.ar; nlags = options_ident.ar;
advanced = options_ident.advanced; advanced = options_ident.advanced;
replic = options_ident.replic; replic = options_ident.replic;
...@@ -138,7 +141,7 @@ error_indicator.identification_spectrum=0; ...@@ -138,7 +141,7 @@ error_indicator.identification_spectrum=0;
if info(1) == 0 %no errors in solution if info(1) == 0 %no errors in solution
% Compute parameter Jacobians for identification analysis % Compute parameter Jacobians for identification analysis
[MEAN, dMEAN, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dMOMENTS, dSPECTRUM, dSPECTRUM_NO_MEAN, dMINIMAL, derivatives_info] = get_identification_jacobians(estim_params_, M_, oo_, options_, options_ident, indpmodel, indpstderr, indpcorr, indvobs); [~, ~, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dMOMENTS, dSPECTRUM, dSPECTRUM_NO_MEAN, dMINIMAL, derivatives_info] = identification.get_jacobians(estim_params_, M_, options_, options_ident, indpmodel, indpstderr, indpcorr, indvobs, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
if isempty(dMINIMAL) if isempty(dMINIMAL)
% Komunjer and Ng is not computed if (1) minimality conditions are not fullfilled or (2) there are more shocks and measurement errors than observables, so we need to reset options % Komunjer and Ng is not computed if (1) minimality conditions are not fullfilled or (2) there are more shocks and measurement errors than observables, so we need to reset options
error_indicator.identification_minimal = 1; error_indicator.identification_minimal = 1;
...@@ -202,7 +205,7 @@ if info(1) == 0 %no errors in solution ...@@ -202,7 +205,7 @@ if info(1) == 0 %no errors in solution
options_ident_local.no_identification_spectrum = 1; %do not recompute dSPECTRUM options_ident_local.no_identification_spectrum = 1; %do not recompute dSPECTRUM
options_ident_local.ar = nlags; %store new lag number options_ident_local.ar = nlags; %store new lag number
options_.ar = nlags; %store new lag number options_.ar = nlags; %store new lag number
[~, ~, ~, ~, ~, ~, MOMENTS, dMOMENTS, ~, ~, ~, ~] = get_identification_jacobians(estim_params_, M_, oo_, options_, options_ident_local, indpmodel, indpstderr, indpcorr, indvobs); [~, ~, ~, ~, ~, ~, MOMENTS, dMOMENTS] = identification.get_jacobians(estim_params_, M_, options_, options_ident_local, indpmodel, indpstderr, indpcorr, indvobs, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
ind_dMOMENTS = (find(max(abs(dMOMENTS'),[],1) > tol_deriv)); %new index with non-zero rows ind_dMOMENTS = (find(max(abs(dMOMENTS'),[],1) > tol_deriv)); %new index with non-zero rows
end end
...@@ -293,25 +296,34 @@ if info(1) == 0 %no errors in solution ...@@ -293,25 +296,34 @@ 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.
[fval, info, cost_flag, DLIK, AHess, ys, trend_coeff, M_, options_, bayestopt_, oo_] = dsge_likelihood(params', dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, bounds, oo_, derivatives_info); %non-used output variables need to be set for octave for some reason [~, 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
%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
if iflag || rank(AHess)>rank(tildaM) if iflag || rank(AHess)>rank(tildaM)
[ide_hess.cond, ide_hess.rank, ide_hess.ind0, ide_hess.indno, ide_hess.ino, ide_hess.Mco, ide_hess.Pco] = identification_checks(AHess, 0, tol_rank, tol_sv, totparam_nbr); [ide_hess.cond, ide_hess.rank, ide_hess.ind0, ide_hess.indno, ide_hess.ino, ide_hess.Mco, ide_hess.Pco] = identification.checks(AHess, 0, tol_rank, tol_sv, totparam_nbr);
else %use normalized version if possible else %use normalized version if possible
[ide_hess.cond, ide_hess.rank, ide_hess.ind0, ide_hess.indno, ide_hess.ino, ide_hess.Mco, ide_hess.Pco] = identification_checks(tildaM, 0, tol_rank, tol_sv, totparam_nbr); [ide_hess.cond, ide_hess.rank, ide_hess.ind0, ide_hess.indno, ide_hess.ino, ide_hess.Mco, ide_hess.Pco] = identification.checks(tildaM, 0, tol_rank, tol_sv, totparam_nbr);
end end
indok = find(max(ide_hess.indno,[],1)==0); indok = find(max(ide_hess.indno,[],1)==0);
ide_uncert_unnormaliz(indok) = sqrt(diag(inv(AHess(indok,indok))))'; ide_uncert_unnormaliz(indok) = sqrt(diag(inv(AHess(indok,indok))))';
...@@ -321,7 +333,7 @@ if info(1) == 0 %no errors in solution ...@@ -321,7 +333,7 @@ if info(1) == 0 %no errors in solution
diag_chh = sum(si_dREDUCEDFORM(:,ind1)'.*temp1)'; diag_chh = sum(si_dREDUCEDFORM(:,ind1)'.*temp1)';
ind1 = ind1(ind1>stderrparam_nbr+corrparam_nbr); ind1 = ind1(ind1>stderrparam_nbr+corrparam_nbr);
cdynamic = si_dDYNAMIC(:,ind1-stderrparam_nbr-corrparam_nbr)*((AHess(ind1,ind1))\si_dDYNAMIC(:,ind1-stderrparam_nbr-corrparam_nbr)'); cdynamic = si_dDYNAMIC(:,ind1-stderrparam_nbr-corrparam_nbr)*((AHess(ind1,ind1))\si_dDYNAMIC(:,ind1-stderrparam_nbr-corrparam_nbr)');
flag_score = 1; %this is used for the title in plot_identification.m flag_score = 1; %this is used for the title in identification.plot.m
catch catch
%Asymptotic Hessian via simulation %Asymptotic Hessian via simulation
if options_.order > 1 if options_.order > 1
...@@ -332,7 +344,7 @@ if info(1) == 0 %no errors in solution ...@@ -332,7 +344,7 @@ if info(1) == 0 %no errors in solution
options_.periods = periods+100; options_.periods = periods+100;
end end
replic = max([replic, length(ind_dMOMENTS)*3]); replic = max([replic, length(ind_dMOMENTS)*3]);
cmm = simulated_moment_uncertainty(ind_dMOMENTS, periods, replic,options_,M_,oo_); %covariance matrix of moments cmm = identification.simulated_moment_uncertainty(ind_dMOMENTS, periods, replic,options_,M_,oo_); %covariance matrix of moments
sd = sqrt(diag(cmm)); sd = sqrt(diag(cmm));
cc = cmm./(sd*sd'); cc = cmm./(sd*sd');
[VV,DD,WW] = eig(cc); [VV,DD,WW] = eig(cc);
...@@ -340,15 +352,15 @@ if info(1) == 0 %no errors in solution ...@@ -340,15 +352,15 @@ 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'));
if iflag || rank(MIM)>rank(tildaM) if iflag || rank(MIM)>rank(tildaM)
[ide_hess.cond, ide_hess.rank, ide_hess.ind0, ide_hess.indno, ide_hess.ino, ide_hess.Mco, ide_hess.Pco] = identification_checks(MIM, 0, tol_rank, tol_sv, totparam_nbr); [ide_hess.cond, ide_hess.rank, ide_hess.ind0, ide_hess.indno, ide_hess.ino, ide_hess.Mco, ide_hess.Pco] = identification.checks(MIM, 0, tol_rank, tol_sv, totparam_nbr);
else %use normalized version if possible else %use normalized version if possible
[ide_hess.cond, ide_hess.rank, ide_hess.ind0, ide_hess.indno, ide_hess.ino, ide_hess.Mco, ide_hess.Pco] = identification_checks(tildaM, 0, tol_rank, tol_sv, totparam_nbr); [ide_hess.cond, ide_hess.rank, ide_hess.ind0, ide_hess.indno, ide_hess.ino, ide_hess.Mco, ide_hess.Pco] = identification.checks(tildaM, 0, tol_rank, tol_sv, totparam_nbr);
end end
indok = find(max(ide_hess.indno,[],1)==0); indok = find(max(ide_hess.indno,[],1)==0);
ind1 = find(ide_hess.ind0); ind1 = find(ide_hess.ind0);
...@@ -359,7 +371,7 @@ if info(1) == 0 %no errors in solution ...@@ -359,7 +371,7 @@ if info(1) == 0 %no errors in solution
if ~isempty(indok) if ~isempty(indok)
ide_uncert_unnormaliz(indok) = (sqrt(diag(inv(tildaM(indok,indok))))./deltaM(indok))'; %sqrt(diag(inv(MIM(indok,indok))))'; ide_uncert_unnormaliz(indok) = (sqrt(diag(inv(tildaM(indok,indok))))./deltaM(indok))'; %sqrt(diag(inv(MIM(indok,indok))))';
end end
flag_score = 0; %this is used for the title in plot_identification.m flag_score = 0; %this is used for the title in identification.plot.m
end % end of computing sample information matrix for identification strength measure end % end of computing sample information matrix for identification strength measure
ide_strength_dMOMENTS(indok) = (1./(ide_uncert_unnormaliz(indok)'./abs(params(indok)'))); %this is s_i in Ratto and Iskrev (2011, p.13) ide_strength_dMOMENTS(indok) = (1./(ide_uncert_unnormaliz(indok)'./abs(params(indok)'))); %this is s_i in Ratto and Iskrev (2011, p.13)
...@@ -371,7 +383,7 @@ if info(1) == 0 %no errors in solution ...@@ -371,7 +383,7 @@ if info(1) == 0 %no errors in solution
if size(quant,1)==1 if size(quant,1)==1
si_dMOMENTSnorm = abs(quant).*normaliz_prior_std; si_dMOMENTSnorm = abs(quant).*normaliz_prior_std;
else else
si_dMOMENTSnorm = vnorm(quant).*normaliz_prior_std; si_dMOMENTSnorm = identification.vnorm(quant).*normaliz_prior_std;
end end
iy = find(diag_chh); iy = find(diag_chh);
ind_dREDUCEDFORM = ind_dREDUCEDFORM(iy); ind_dREDUCEDFORM = ind_dREDUCEDFORM(iy);
...@@ -381,7 +393,7 @@ if info(1) == 0 %no errors in solution ...@@ -381,7 +393,7 @@ if info(1) == 0 %no errors in solution
if size(quant,1)==1 if size(quant,1)==1
si_dREDUCEDFORMnorm = abs(quant).*normaliz_prior_std; si_dREDUCEDFORMnorm = abs(quant).*normaliz_prior_std;
else else
si_dREDUCEDFORMnorm = vnorm(quant).*normaliz_prior_std; si_dREDUCEDFORMnorm = identification.vnorm(quant).*normaliz_prior_std;
end end
else else
si_dREDUCEDFORMnorm = []; si_dREDUCEDFORMnorm = [];
...@@ -395,7 +407,7 @@ if info(1) == 0 %no errors in solution ...@@ -395,7 +407,7 @@ if info(1) == 0 %no errors in solution
if size(quant,1)==1 if size(quant,1)==1
si_dDYNAMICnorm = abs(quant).*normaliz_prior_std(stderrparam_nbr+corrparam_nbr+1:end); si_dDYNAMICnorm = abs(quant).*normaliz_prior_std(stderrparam_nbr+corrparam_nbr+1:end);
else else
si_dDYNAMICnorm = vnorm(quant).*normaliz_prior_std(stderrparam_nbr+corrparam_nbr+1:end); si_dDYNAMICnorm = identification.vnorm(quant).*normaliz_prior_std(stderrparam_nbr+corrparam_nbr+1:end);
end end
else else
si_dDYNAMICnorm=[]; si_dDYNAMICnorm=[];
...@@ -461,11 +473,11 @@ if info(1) == 0 %no errors in solution ...@@ -461,11 +473,11 @@ if info(1) == 0 %no errors in solution
ide_moments.MOMENTS = MOMENTS; ide_moments.MOMENTS = MOMENTS;
if advanced if advanced
% here we do not normalize (i.e. we set norm_dMOMENTS=1) as the OLS in ident_bruteforce is very sensitive to norm_dMOMENTS % here we do not normalize (i.e. we set norm_dMOMENTS=1) as the OLS in identification.bruteforce is very sensitive to norm_dMOMENTS
[ide_moments.pars, ide_moments.cosndMOMENTS] = ident_bruteforce(dMOMENTS(ind_dMOMENTS,:), max_dim_cova_group, options_.TeX, options_ident.name_tex, options_ident.tittxt, tol_deriv); [ide_moments.pars, ide_moments.cosndMOMENTS] = identification.bruteforce(M_.dname,M_.fname,dMOMENTS(ind_dMOMENTS,:), max_dim_cova_group, options_.TeX, options_ident.name_tex, options_ident.tittxt, tol_deriv);
end end
%here we focus on the unnormalized S and V, which is then used in plot_identification.m and for prior_mc > 1 %here we focus on the unnormalized S and V, which is then used in identification.plot.m and for prior_mc > 1
[~, S, V] = svd(dMOMENTS(ind_dMOMENTS,:),0); [~, S, V] = svd(dMOMENTS(ind_dMOMENTS,:),0);
if size(S,1) == 1 if size(S,1) == 1
S = S(1); % edge case that S is not a matrix but a row vector S = S(1); % edge case that S is not a matrix but a row vector
...@@ -518,9 +530,9 @@ if info(1) == 0 %no errors in solution ...@@ -518,9 +530,9 @@ if info(1) == 0 %no errors in solution
%% Perform identification checks, i.e. find out which parameters are involved %% Perform identification checks, i.e. find out which parameters are involved
if checks_via_subsets if checks_via_subsets
% identification_checks_via_subsets is only for debugging % identification.checks_via_subsets is only for debugging
[ide_dynamic, ide_reducedform, ide_moments, ide_spectrum, ide_minimal] = ... [ide_dynamic, ide_reducedform, ide_moments, ide_spectrum, ide_minimal] = ...
identification_checks_via_subsets(ide_dynamic, ide_reducedform, ide_moments, ide_spectrum, ide_minimal, totparam_nbr, modparam_nbr, options_ident, error_indicator); identification.checks_via_subsets(ide_dynamic, ide_reducedform, ide_moments, ide_spectrum, ide_minimal, totparam_nbr, modparam_nbr, options_ident, error_indicator);
if ~error_indicator.identification_minimal if ~error_indicator.identification_minimal
ide_minimal.minimal_state_space=1; ide_minimal.minimal_state_space=1;
else else
...@@ -528,19 +540,19 @@ if info(1) == 0 %no errors in solution ...@@ -528,19 +540,19 @@ if info(1) == 0 %no errors in solution
end end
else else
[ide_dynamic.cond, ide_dynamic.rank, ide_dynamic.ind0, ide_dynamic.indno, ide_dynamic.ino, ide_dynamic.Mco, ide_dynamic.Pco, ide_dynamic.jweak, ide_dynamic.jweak_pair] = ... [ide_dynamic.cond, ide_dynamic.rank, ide_dynamic.ind0, ide_dynamic.indno, ide_dynamic.ino, ide_dynamic.Mco, ide_dynamic.Pco, ide_dynamic.jweak, ide_dynamic.jweak_pair] = ...
identification_checks(dDYNAMIC(ind_dDYNAMIC,:)./norm_dDYNAMIC, 1, tol_rank, tol_sv, modparam_nbr); identification.checks(dDYNAMIC(ind_dDYNAMIC,:)./norm_dDYNAMIC, 1, tol_rank, tol_sv, modparam_nbr);
if ~options_ident.no_identification_reducedform && ~error_indicator.identification_reducedform if ~options_ident.no_identification_reducedform && ~error_indicator.identification_reducedform
[ide_reducedform.cond, ide_reducedform.rank, ide_reducedform.ind0, ide_reducedform.indno, ide_reducedform.ino, ide_reducedform.Mco, ide_reducedform.Pco, ide_reducedform.jweak, ide_reducedform.jweak_pair] = ... [ide_reducedform.cond, ide_reducedform.rank, ide_reducedform.ind0, ide_reducedform.indno, ide_reducedform.ino, ide_reducedform.Mco, ide_reducedform.Pco, ide_reducedform.jweak, ide_reducedform.jweak_pair] = ...
identification_checks(dREDUCEDFORM(ind_dREDUCEDFORM,:)./norm_dREDUCEDFORM, 1, tol_rank, tol_sv, totparam_nbr); identification.checks(dREDUCEDFORM(ind_dREDUCEDFORM,:)./norm_dREDUCEDFORM, 1, tol_rank, tol_sv, totparam_nbr);
end end
if ~options_ident.no_identification_moments && ~error_indicator.identification_moments if ~options_ident.no_identification_moments && ~error_indicator.identification_moments
[ide_moments.cond, ide_moments.rank, ide_moments.ind0, ide_moments.indno, ide_moments.ino, ide_moments.Mco, ide_moments.Pco, ide_moments.jweak, ide_moments.jweak_pair] = ... [ide_moments.cond, ide_moments.rank, ide_moments.ind0, ide_moments.indno, ide_moments.ino, ide_moments.Mco, ide_moments.Pco, ide_moments.jweak, ide_moments.jweak_pair] = ...
identification_checks(dMOMENTS(ind_dMOMENTS,:)./norm_dMOMENTS, 1, tol_rank, tol_sv, totparam_nbr); identification.checks(dMOMENTS(ind_dMOMENTS,:)./norm_dMOMENTS, 1, tol_rank, tol_sv, totparam_nbr);
end end
if ~options_ident.no_identification_minimal if ~options_ident.no_identification_minimal
if ~error_indicator.identification_minimal if ~error_indicator.identification_minimal
[ide_minimal.cond, ide_minimal.rank, ide_minimal.ind0, ide_minimal.indno, ide_minimal.ino, ide_minimal.Mco, ide_minimal.Pco, ide_minimal.jweak, ide_minimal.jweak_pair] = ... [ide_minimal.cond, ide_minimal.rank, ide_minimal.ind0, ide_minimal.indno, ide_minimal.ino, ide_minimal.Mco, ide_minimal.Pco, ide_minimal.jweak, ide_minimal.jweak_pair] = ...
identification_checks(dMINIMAL(ind_dMINIMAL,:)./norm_dMINIMAL, 2, tol_rank, tol_sv, totparam_nbr); identification.checks(dMINIMAL(ind_dMINIMAL,:)./norm_dMINIMAL, 2, tol_rank, tol_sv, totparam_nbr);
ide_minimal.minimal_state_space=1; ide_minimal.minimal_state_space=1;
else else
ide_minimal.minimal_state_space=0; ide_minimal.minimal_state_space=0;
...@@ -548,7 +560,7 @@ if info(1) == 0 %no errors in solution ...@@ -548,7 +560,7 @@ if info(1) == 0 %no errors in solution
end end
if ~options_ident.no_identification_spectrum && ~error_indicator.identification_spectrum if ~options_ident.no_identification_spectrum && ~error_indicator.identification_spectrum
[ide_spectrum.cond, ide_spectrum.rank, ide_spectrum.ind0, ide_spectrum.indno, ide_spectrum.ino, ide_spectrum.Mco, ide_spectrum.Pco, ide_spectrum.jweak, ide_spectrum.jweak_pair] = ... [ide_spectrum.cond, ide_spectrum.rank, ide_spectrum.ind0, ide_spectrum.indno, ide_spectrum.ino, ide_spectrum.Mco, ide_spectrum.Pco, ide_spectrum.jweak, ide_spectrum.jweak_pair] = ...
identification_checks(tilda_dSPECTRUM, 3, tol_rank, tol_sv, totparam_nbr); identification.checks(tilda_dSPECTRUM, 3, tol_rank, tol_sv, totparam_nbr);
end end
end end
end end
matlab/ident_bruteforce.m → matlab/+identification/bruteforce.m
View file @ 4a32acc2
function [pars, cosnJ] = ident_bruteforce(J, max_dim_cova_group, TeX, name_tex, tittxt, tol_deriv) function [pars, cosnJ] = bruteforce(dname,fname,J, max_dim_cova_group, TeX, name_tex, tittxt, tol_deriv)
% function [pars, cosnJ] = ident_bruteforce(J,n,TeX, pnames_TeX,tittxt) % [pars, cosnJ] = bruteforce(dname,fname,J, max_dim_cova_group, TeX, name_tex, tittxt, tol_deriv)
% ------------------------------------------------------------------------- % -------------------------------------------------------------------------
% given the Jacobian matrix J of moment derivatives w.r.t. parameters % given the Jacobian matrix J of moment derivatives w.r.t. parameters
% computes, for each column of J, the groups of columns from 1 to n that % computes, for each column of J, the groups of columns from 1 to n that
...@@ -18,9 +18,9 @@ function [pars, cosnJ] = ident_bruteforce(J, max_dim_cova_group, TeX, name_tex, ...@@ -18,9 +18,9 @@ function [pars, cosnJ] = ident_bruteforce(J, max_dim_cova_group, TeX, name_tex,
% cosnJ : cosn of each column with the selected group of columns % cosnJ : cosn of each column with the selected group of columns
% ------------------------------------------------------------------------- % -------------------------------------------------------------------------
% This function is called by % This function is called by
% * identification_analysis.m % * identification.analysis.m
% ========================================================================= % =========================================================================
% Copyright © 2009-2019 Dynare Team % Copyright © 2009-2023 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
...@@ -38,20 +38,18 @@ function [pars, cosnJ] = ident_bruteforce(J, max_dim_cova_group, TeX, name_tex, ...@@ -38,20 +38,18 @@ function [pars, cosnJ] = ident_bruteforce(J, max_dim_cova_group, TeX, name_tex,
% along with Dynare. If not, see <https://www.gnu.org/licenses/>. % along with Dynare. If not, see <https://www.gnu.org/licenses/>.
% ========================================================================= % =========================================================================
global M_ options_ OutputDirectoryName = CheckPath('identification',dname);
OutputDirectoryName = CheckPath('identification',M_.dname);
totparam_nbr = size(J,2); % number of parameters totparam_nbr = size(J,2); % number of parameters
if nargin<2 || isempty(max_dim_cova_group) if nargin<4 || isempty(max_dim_cova_group)
max_dim_cova_group = 4; % max n-tuple max_dim_cova_group = 4; % max n-tuple
end end
if nargin<3 || isempty(TeX) if nargin<5 || isempty(TeX)
TeX = 0; % no Tex output TeX = 0; % no Tex output
tittxt=''; tittxt='';
end end
if nargin < 6 if nargin < 8
tol_deriv = 1.e-8; tol_deriv = 1.e-8;
end end
...@@ -69,7 +67,7 @@ for ll = 1:max_dim_cova_group ...@@ -69,7 +67,7 @@ for ll = 1:max_dim_cova_group
cosnJ2=zeros(size(tmp2,1),1); cosnJ2=zeros(size(tmp2,1),1);
b=[]; b=[];
for jj = 1:size(tmp2,1) for jj = 1:size(tmp2,1)
[cosnJ2(jj,1), b(:,jj)] = cosn([J(:,ii),J(:,tmp2(jj,:))]); [cosnJ2(jj,1), b(:,jj)] = identification.cosn([J(:,ii),J(:,tmp2(jj,:))]);
end end
cosnJ(ii,ll) = max(cosnJ2(:,1)); cosnJ(ii,ll) = max(cosnJ2(:,1));
if cosnJ(ii,ll)>tol_deriv if cosnJ(ii,ll)>tol_deriv
...@@ -87,7 +85,7 @@ for ll = 1:max_dim_cova_group ...@@ -87,7 +85,7 @@ for ll = 1:max_dim_cova_group
end end
dyn_waitbar_close(h); dyn_waitbar_close(h);
if TeX if TeX
filename = [OutputDirectoryName '/' M_.fname '_collin_patterns_',tittxt1,'_' int2str(ll) '.tex']; filename = [OutputDirectoryName '/' fname '_collin_patterns_',tittxt1,'_' int2str(ll) '.tex'];
fidTeX = fopen(filename,'w'); fidTeX = fopen(filename,'w');
fprintf(fidTeX,'%% TeX-table generated by ident_bruteforce (Dynare).\n'); fprintf(fidTeX,'%% TeX-table generated by ident_bruteforce (Dynare).\n');
fprintf(fidTeX,['%% Collinearity patterns with ',int2str(ll),' parameter(s): ',tittxt,'\n']); fprintf(fidTeX,['%% Collinearity patterns with ',int2str(ll),' parameter(s): ',tittxt,'\n']);
...@@ -114,12 +112,12 @@ for ll = 1:max_dim_cova_group ...@@ -114,12 +112,12 @@ for ll = 1:max_dim_cova_group
plist=''; plist='';
for ii=1:ll for ii=1:ll
if ~isnan(pars{i,ll}(ii)) if ~isnan(pars{i,ll}(ii))
plist = [plist ' $' name_tex{pars{i,ll}(ii)} '\;\; $ ']; plist = [plist ' ' name_tex{pars{i,ll}(ii)} '\;\; '];
else else
plist = [plist ' ---- ']; plist = [plist ' ---- '];
end end
end end
fprintf(fidTeX,'$%s$ & [%s] & %7.3f \\\\ \n',... fprintf(fidTeX,'%s & [%s] & %7.3f \\\\ \n',...
name_tex{i},... name_tex{i},...
plist,... plist,...
cosnJ(i,ll)); cosnJ(i,ll));
......
matlab/identification_checks.m → matlab/+identification/checks.m
View file @ 4a32acc2
function [condX, rankX, ind0, indno, ixno, Mco, Pco, jweak, jweak_pair] = identification_checks(X, test_flag, tol_rank, tol_sv, param_nbr) function [condX, rankX, ind0, indno, ixno, Mco, Pco, jweak, jweak_pair] = checks(X, test_flag, tol_rank, tol_sv, param_nbr)
% function [condX, rankX, ind0, indno, ixno, Mco, Pco, jweak, jweak_pair] = identification_checks(X, test_flag, tol_rank, tol_sv, param_nbr) % function [condX, rankX, ind0, indno, ixno, Mco, Pco, jweak, jweak_pair] = checks(X, test_flag, tol_rank, tol_sv, param_nbr)
% ------------------------------------------------------------------------- % -------------------------------------------------------------------------
% Checks rank criteria of identification tests and finds out parameter sets % Checks rank criteria of identification tests and finds out parameter sets
% that are not identifiable via the nullspace, pairwise correlation % that are not identifiable via the nullspace, pairwise correlation
...@@ -10,7 +10,7 @@ function [condX, rankX, ind0, indno, ixno, Mco, Pco, jweak, jweak_pair] = identi ...@@ -10,7 +10,7 @@ function [condX, rankX, ind0, indno, ixno, Mco, Pco, jweak, jweak_pair] = identi
% 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
...@@ -24,10 +24,10 @@ function [condX, rankX, ind0, indno, ixno, Mco, Pco, jweak, jweak_pair] = identi ...@@ -24,10 +24,10 @@ function [condX, rankX, ind0, indno, ixno, Mco, Pco, jweak, jweak_pair] = identi
% * jweak_pair [(vech) matrix] gives 1 if a couple parameters has Pco=1 (with tolerance tol_rank) % * jweak_pair [(vech) matrix] gives 1 if a couple parameters has Pco=1 (with tolerance tol_rank)
% ------------------------------------------------------------------------- % -------------------------------------------------------------------------
% This function is called by % This function is called by
% * identification_analysis.m % * identification.analysis.m
% ------------------------------------------------------------------------- % -------------------------------------------------------------------------
% This function calls % This function calls
% * cosn % * identification.cosn
% * dyn_vech % * dyn_vech
% * vnorm % * vnorm
% ========================================================================= % =========================================================================
...@@ -75,7 +75,7 @@ end ...@@ -75,7 +75,7 @@ end
% find non-zero columns at machine precision % find non-zero columns at machine precision
if size(Xpar,1) > 1 if size(Xpar,1) > 1
ind1 = find(vnorm(Xpar) >= eps); ind1 = find(identification.vnorm(Xpar) >= eps);
else else
ind1 = find(abs(Xpar) >= eps); % if only one parameter ind1 = find(abs(Xpar) >= eps); % if only one parameter
end end
...@@ -86,7 +86,7 @@ else ...@@ -86,7 +86,7 @@ else
Xparnonzero = Xpar(:,ind1); % focus on non-zero columns Xparnonzero = Xpar(:,ind1); % focus on non-zero columns
end end
[eu, ee2, ee1] = svd( [Xparnonzero Xrest], 0 ); [~, ~, ee1] = svd( [Xparnonzero Xrest], 0 );
condX = cond([Xparnonzero Xrest]); condX = cond([Xparnonzero Xrest]);
rankX = rank(X, tol_rank); rankX = rank(X, tol_rank);
icheck = 0; %initialize flag which is equal to 0 if we already found all single parameters that are not identified icheck = 0; %initialize flag which is equal to 0 if we already found all single parameters that are not identified
...@@ -94,7 +94,7 @@ if param_nbr > 0 && (rankX<rankrequired) ...@@ -94,7 +94,7 @@ if param_nbr > 0 && (rankX<rankrequired)
% search for singular values associated to ONE individual parameter % search for singular values associated to ONE individual parameter
% Compute an orthonormal basis for the null space using the columns of ee1 that correspond % Compute an orthonormal basis for the null space using the columns of ee1 that correspond
% to singular values equal to zero and associated to an individual parameter % to singular values equal to zero and associated to an individual parameter
ee0 = [rankX+1:size([Xparnonzero Xrest],2)]; %look into last columns with singular values of problematic parameter sets (except single parameters) ee0 = rankX+1:size([Xparnonzero Xrest],2); %look into last columns with singular values of problematic parameter sets (except single parameters)
ind11 = ones(length(ind1),1); %initialize ind11 = ones(length(ind1),1); %initialize
for j=1:length(ee0) for j=1:length(ee0)
% check if nullspace is spanned by only one parameter % check if nullspace is spanned by only one parameter
...@@ -118,7 +118,7 @@ if icheck ...@@ -118,7 +118,7 @@ if icheck
else else
Xparnonzero = Xpar(:,ind1); % focus on non-zero columns Xparnonzero = Xpar(:,ind1); % focus on non-zero columns
end end
[eu, ee2, ee1] = svd( [Xparnonzero Xrest], 0 ); [~, ~, ee1] = svd( [Xparnonzero Xrest], 0 );
condX = cond([Xparnonzero Xrest]); condX = cond([Xparnonzero Xrest]);
rankX = rank(X,tol_rank); rankX = rank(X,tol_rank);
end end
...@@ -141,7 +141,7 @@ if test_flag == 0 || test_flag == 3 % G is a Gram matrix and hence should be a c ...@@ -141,7 +141,7 @@ if test_flag == 0 || test_flag == 3 % G is a Gram matrix and hence should be a c
else else
Mco = NaN(param_nbr,1); Mco = NaN(param_nbr,1);
for ii = 1:size(Xparnonzero,2) for ii = 1:size(Xparnonzero,2)
Mco(ind1(ii),:) = cosn([Xparnonzero(:,ii) , Xparnonzero(:,find([1:1:size(Xparnonzero,2)]~=ii)), Xrest]); Mco(ind1(ii),:) = identification.cosn([Xparnonzero(:,ii) , Xparnonzero(:,find([1:1:size(Xparnonzero,2)]~=ii)), Xrest]);
end end
end end
...@@ -151,9 +151,9 @@ if param_nbr>0 && (rankX<rankrequired || min(1-Mco)<tol_rank) ...@@ -151,9 +151,9 @@ if param_nbr>0 && (rankX<rankrequired || min(1-Mco)<tol_rank)
if length(ind1)<param_nbr if length(ind1)<param_nbr
% single parameters with zero columns % single parameters with zero columns
ixno = ixno + 1; ixno = ixno + 1;
indno(ixno,:) = (~ismember([1:param_nbr],ind1)); indno(ixno,:) = (~ismember(1:param_nbr,ind1));
end end
ee0 = [rankX+1:size([Xparnonzero Xrest],2)]; %look into last columns with singular values of problematic parameter sets (except single parameters) ee0 = rankX+1:size([Xparnonzero Xrest],2); %look into last columns with singular values of problematic parameter sets (except single parameters)
for j=1:length(ee0) for j=1:length(ee0)
% linearly dependent parameters % linearly dependent parameters
ixno = ixno + 1; ixno = ixno + 1;
...@@ -170,13 +170,13 @@ end ...@@ -170,13 +170,13 @@ end
jweak = zeros(1,param_nbr); jweak = zeros(1,param_nbr);
jweak_pair = zeros(param_nbr,param_nbr); jweak_pair = zeros(param_nbr,param_nbr);
if test_flag ~= 0 || test_flag ~= 0 if test_flag ~= 0
% these tests only apply to Jacobians, not to Gram matrices, i.e. Hessian-type or 'covariance' matrices % these tests only apply to Jacobians, not to Gram matrices, i.e. Hessian-type or 'covariance' matrices
Pco = NaN(param_nbr,param_nbr); Pco = NaN(param_nbr,param_nbr);
for ii = 1:size(Xparnonzero,2) for ii = 1:size(Xparnonzero,2)
Pco(ind1(ii),ind1(ii)) = 1; Pco(ind1(ii),ind1(ii)) = 1;
for jj = ii+1:size(Xparnonzero,2) for jj = ii+1:size(Xparnonzero,2)
Pco(ind1(ii),ind1(jj)) = cosn([Xparnonzero(:,ii),Xparnonzero(:,jj),Xrest]); Pco(ind1(ii),ind1(jj)) = identification.cosn([Xparnonzero(:,ii),Xparnonzero(:,jj),Xrest]);
Pco(ind1(jj),ind1(ii)) = Pco(ind1(ii),ind1(jj)); Pco(ind1(jj),ind1(ii)) = Pco(ind1(ii),ind1(jj));
end end
end end
......
matlab/identification_checks_via_subsets.m → matlab/+identification/checks_via_subsets.m
View file @ 4a32acc2
function [ide_dynamic, ide_reducedform, ide_moments, ide_spectrum, ide_minimal] = identification_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] = identification_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
...@@ -50,7 +50,7 @@ function [ide_dynamic, ide_reducedform, ide_moments, ide_spectrum, ide_minimal] ...@@ -50,7 +50,7 @@ function [ide_dynamic, ide_reducedform, ide_moments, ide_spectrum, ide_minimal]
% * rank: [integer] rank of Jacobian % * rank: [integer] rank of Jacobian
% ------------------------------------------------------------------------- % -------------------------------------------------------------------------
% This function is called by % This function is called by
% * identification_analysis.m % * identification.analysis.m
% ========================================================================= % =========================================================================
% Copyright © 2019-2021 Dynare Team % Copyright © 2019-2021 Dynare Team
% %
...@@ -161,7 +161,7 @@ end ...@@ -161,7 +161,7 @@ end
% initialize for spectrum criteria % initialize for spectrum criteria
if ~no_identification_spectrum && ~error_indicator.identification_spectrum if ~no_identification_spectrum && ~error_indicator.identification_spectrum
dSPECTRUM = ide_spectrum.tilda_dSPECTRUM; %tilda dSPECTRUM is normalized dSPECTRUM matrix in identification_analysis.m dSPECTRUM = ide_spectrum.tilda_dSPECTRUM; %tilda dSPECTRUM is normalized dSPECTRUM matrix in identification.analysis.m
%alternative normalization %alternative normalization
%dSPECTRUM = ide_spectrum.dSPECTRUM; %dSPECTRUM = ide_spectrum.dSPECTRUM;
%dSPECTRUM(ide_spectrum.ind_dSPECTRUM,:) = dSPECTRUM(ide_spectrum.ind_dSPECTRUM,:)./ide_spectrum.norm_dSPECTRUM; %normalize %dSPECTRUM(ide_spectrum.ind_dSPECTRUM,:) = dSPECTRUM(ide_spectrum.ind_dSPECTRUM,:)./ide_spectrum.norm_dSPECTRUM; %normalize
......
matlab/cosn.m → matlab/+identification/cosn.m
View file @ 4a32acc2
...@@ -17,7 +17,7 @@ function [co, b, yhat] = cosn(H) ...@@ -17,7 +17,7 @@ function [co, b, yhat] = cosn(H)
% * y [n by 1] predicted endogenous values given ols estimation % * y [n by 1] predicted endogenous values given ols estimation
% ------------------------------------------------------------------------- % -------------------------------------------------------------------------
% This function is called by % This function is called by
% * identification_checks.m % * identification.checks.m
% * ident_bruteforce.m % * ident_bruteforce.m
% ========================================================================= % =========================================================================
% Copyright © 2008-2019 Dynare Team % Copyright © 2008-2019 Dynare Team
......
matlab/disp_identification.m → matlab/+identification/display.m
View file @ 4a32acc2
function disp_identification(pdraws, ide_reducedform, ide_moments, ide_spectrum, ide_minimal, name, options_ident) function display(pdraws, ide_reducedform, ide_moments, ide_spectrum, ide_minimal, name, options_ident)
% disp_identification(pdraws, ide_reducedform, ide_moments, ide_spectrum, ide_minimal, name, options_ident) % display(pdraws, ide_reducedform, ide_moments, ide_spectrum, ide_minimal, name, options_ident)
% ------------------------------------------------------------------------- % -------------------------------------------------------------------------
% This function displays all identification analysis to the command line % This function displays all identification analysis to the command line
% ========================================================================= % =========================================================================
...@@ -26,7 +26,7 @@ function disp_identification(pdraws, ide_reducedform, ide_moments, ide_spectrum, ...@@ -26,7 +26,7 @@ function disp_identification(pdraws, ide_reducedform, ide_moments, ide_spectrum,
% * all output is printed on the command line % * all output is printed on the command line
% ------------------------------------------------------------------------- % -------------------------------------------------------------------------
% This function is called by % This function is called by
% * dynare_identification.m % * identification.run
% ========================================================================= % =========================================================================
% Copyright © 2010-2021 Dynare Team % Copyright © 2010-2021 Dynare Team
% %
...@@ -54,7 +54,7 @@ tol_rank = options_ident.tol_rank; ...@@ -54,7 +54,7 @@ tol_rank = options_ident.tol_rank;
checks_via_subsets = options_ident.checks_via_subsets; checks_via_subsets = options_ident.checks_via_subsets;
%% Display settings %% Display settings
disp([' ']), disp(' '),
fprintf('Note that differences in the criteria could be due to numerical settings,\n') fprintf('Note that differences in the criteria could be due to numerical settings,\n')
fprintf('numerical errors or the method used to find problematic parameter sets.\n') fprintf('numerical errors or the method used to find problematic parameter sets.\n')
fprintf('Settings:\n') fprintf('Settings:\n')
...@@ -157,7 +157,7 @@ for jide = 1:4 ...@@ -157,7 +157,7 @@ for jide = 1:4
disp(' !!!WARNING!!!'); disp(' !!!WARNING!!!');
if SampleSize>1 if SampleSize>1
if non_minimal_state_space_error if non_minimal_state_space_error
fprintf(['\n The minimal state space could not be computed for %u out of %u cases.\n'],SampleSize-EffectiveSampleSize,SampleSize); fprintf('\n The minimal state space could not be computed for %u out of %u cases.\n',SampleSize-EffectiveSampleSize,SampleSize);
end end
if jide==2 if jide==2
if sum(ide.ino & ide.minimal_state_space)>0 if sum(ide.ino & ide.minimal_state_space)>0
...@@ -207,7 +207,7 @@ for jide = 1:4 ...@@ -207,7 +207,7 @@ for jide = 1:4
end end
end end
%% display problematic parameters computed by identification_checks_via_subsets %% display problematic parameters computed by identification.checks_via_subsets
elseif checks_via_subsets elseif checks_via_subsets
if ide.rank < size(Jacob,2) if ide.rank < size(Jacob,2)
no_warning_message_display = 0; no_warning_message_display = 0;
......
matlab/fjaco.m → matlab/+identification/fjaco.m
View file @ 4a32acc2
function fjac = fjaco(f,x,varargin) function fjac = fjaco(f,x,varargin)
% fjac = fjaco(f,x,varargin)
% FJACO Computes two-sided finite difference Jacobian % FJACO Computes two-sided finite difference Jacobian
% USAGE % USAGE
% fjac = fjaco(f,x,P1,P2,...) % fjac = fjaco(f,x,P1,P2,...)
...@@ -10,7 +10,7 @@ function fjac = fjaco(f,x,varargin) ...@@ -10,7 +10,7 @@ function fjac = fjaco(f,x,varargin)
% OUTPUT % OUTPUT
% fjac : finite difference Jacobian % fjac : finite difference Jacobian
% %
% Copyright © 2010-2017,2019-2020 Dynare Team % Copyright © 2010-2017,2019-2023 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
...@@ -30,8 +30,8 @@ function fjac = fjaco(f,x,varargin) ...@@ -30,8 +30,8 @@ function fjac = fjaco(f,x,varargin)
ff=feval(f,x,varargin{:}); ff=feval(f,x,varargin{:});
tol = eps.^(1/3); %some default value tol = eps.^(1/3); %some default value
if strcmp(func2str(f),'get_perturbation_params_derivs_numerical_objective') || strcmp(func2str(f),'identification_numerical_objective') if strcmp(func2str(f),'identification.get_perturbation_params_derivs_numerical_objective') || strcmp(func2str(f),'identification.numerical_objective')
tol= varargin{5}.dynatol.x; tol= varargin{4}.dynatol.x;
end end
h = tol.*max(abs(x),1); h = tol.*max(abs(x),1);
xh1=x+h; xh0=x-h; xh1=x+h; xh0=x-h;
...@@ -40,12 +40,12 @@ fjac = NaN(length(ff),length(x)); ...@@ -40,12 +40,12 @@ fjac = NaN(length(ff),length(x));
for j=1:length(x) for j=1:length(x)
xx = x; xx = x;
xx(j) = xh1(j); f1=feval(f,xx,varargin{:}); xx(j) = xh1(j); f1=feval(f,xx,varargin{:});
if isempty(f1) && (strcmp(func2str(f),'get_perturbation_params_derivs_numerical_objective') || strcmp(func2str(f),'identification_numerical_objective') ) if isempty(f1) && (strcmp(func2str(f),'identification.get_perturbation_params_derivs_numerical_objective') || strcmp(func2str(f),'identification.numerical_objective') )
[~,info]=feval(f,xx,varargin{:}); [~,info]=feval(f,xx,varargin{:});
disp_info_error_identification_perturbation(info,j); disp_info_error_identification_perturbation(info,j);
end end
xx(j) = xh0(j); f0=feval(f,xx,varargin{:}); xx(j) = xh0(j); f0=feval(f,xx,varargin{:});
if isempty(f0) && (strcmp(func2str(f),'get_perturbation_params_derivs_numerical_objective') || strcmp(func2str(f),'identification_numerical_objective') ) if isempty(f0) && (strcmp(func2str(f),'identification.get_perturbation_params_derivs_numerical_objective') || strcmp(func2str(f),'identification.numerical_objective') )
[~,info]=feval(f,xx,varargin{:}); [~,info]=feval(f,xx,varargin{:});
disp_info_error_identification_perturbation(info,j) disp_info_error_identification_perturbation(info,j)
end end
...@@ -57,9 +57,9 @@ feval(f,x,varargin{:}); ...@@ -57,9 +57,9 @@ feval(f,x,varargin{:});
%Auxiliary functions %Auxiliary functions
function disp_info_error_identification_perturbation(info,j) function disp_info_error_identification_perturbation(info,j)
% there are errors in the solution algorithm % there are errors in the solution algorithm
probl_par = get_the_name(j,varargin{5}.TeX,varargin{3},varargin{2},varargin{5}); probl_par = get_the_name(j,varargin{4}.TeX,varargin{3},varargin{2},varargin{4}.varobs);
skipline() skipline()
message = get_error_message(info,varargin{5}); message = get_error_message(info,varargin{4});
fprintf('Parameter error in numerical two-sided difference method:\n') fprintf('Parameter error in numerical two-sided difference method:\n')
fprintf('Cannot solve the model for %s (info = %d, %s)\n', probl_par, info(1), message); fprintf('Cannot solve the model for %s (info = %d, %s)\n', probl_par, info(1), message);
fprintf('Possible solutions:\n') fprintf('Possible solutions:\n')
......
matlab/get_identification_jacobians.m → matlab/+identification/get_jacobians.m
View file @ 4a32acc2
function [MEAN, dMEAN, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dMOMENTS, dSPECTRUM, dSPECTRUM_NO_MEAN, dMINIMAL, derivatives_info] = get_identification_jacobians(estim_params, M, oo, options, options_ident, indpmodel, indpstderr, indpcorr, indvobs) function [MEAN, dMEAN, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dMOMENTS, dSPECTRUM, dSPECTRUM_NO_MEAN, dMINIMAL, derivatives_info] = get_jacobians(estim_params, M_, options_, options_ident, indpmodel, indpstderr, indpcorr, indvobs, dr, endo_steady_state, exo_steady_state, exo_det_steady_state)
% function [MEAN, dMEAN, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dMOMENTS, dSPECTRUM, dMINIMAL, derivatives_info] = get_identification_jacobians(estim_params, M, oo, options, options_ident, indpmodel, indpstderr, indpcorr, indvobs) % [MEAN, dMEAN, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dMOMENTS, dSPECTRUM, dSPECTRUM_NO_MEAN, dMINIMAL, derivatives_info] = get_jacobians(estim_params, M_, options_, options_ident, indpmodel, indpstderr, indpcorr, indvobs, dr, endo_steady_state, exo_steady_state, exo_det_steady_state)
% previously getJJ.m in Dynare 4.5 % previously getJJ.m in Dynare 4.5
% Sets up the Jacobians needed for identification analysis % Sets up the Jacobians needed for identification analysis
% ========================================================================= % =========================================================================
% INPUTS % INPUTS
% estim_params: [structure] storing the estimation information % estim_params: [structure] storing the estimation information
% M: [structure] storing the model information % M_: [structure] storing the model information
% oo: [structure] storing the reduced-form solution results % options_: [structure] storing the options
% options: [structure] storing the options
% options_ident: [structure] storing the options for identification % options_ident: [structure] storing the options for identification
% indpmodel: [modparam_nbr by 1] index of estimated parameters in M_.params; % indpmodel: [modparam_nbr by 1] index of estimated parameters in M_.params;
% corresponds to model parameters (no stderr and no corr) % corresponds to model parameters (no stderr and no corr)
...@@ -22,6 +21,10 @@ function [MEAN, dMEAN, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dM ...@@ -22,6 +21,10 @@ function [MEAN, dMEAN, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dM
% in the estimated_params block; if estimated_params block % in the estimated_params block; if estimated_params block
% is not available, then no corr parameters are selected % is not available, then no corr parameters are selected
% indvobs: [obs_nbr by 1] index of observed (VAROBS) variables % indvobs: [obs_nbr by 1] index of observed (VAROBS) variables
% dr [structure] Reduced form model.
% endo_steady_state [vector] steady state value for endogenous variables
% exo_steady_state [vector] steady state value for exogenous variables
% exo_det_steady_state [vector] steady state value for exogenous deterministic variables
% ------------------------------------------------------------------------- % -------------------------------------------------------------------------
% OUTPUTS % OUTPUTS
% %
...@@ -81,7 +84,7 @@ function [MEAN, dMEAN, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dM ...@@ -81,7 +84,7 @@ function [MEAN, dMEAN, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dM
% %
% ------------------------------------------------------------------------- % -------------------------------------------------------------------------
% This function is called by % This function is called by
% * identification_analysis.m % * identification.analysis.m
% ------------------------------------------------------------------------- % -------------------------------------------------------------------------
% This function calls % This function calls
% * commutation % * commutation
...@@ -91,7 +94,7 @@ function [MEAN, dMEAN, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dM ...@@ -91,7 +94,7 @@ function [MEAN, dMEAN, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dM
% * fjaco % * fjaco
% * get_perturbation_params_derivs (previously getH) % * get_perturbation_params_derivs (previously getH)
% * get_all_parameters % * get_all_parameters
% * identification_numerical_objective (previously thet2tau) % * identification.numerical_objective (previously thet2tau)
% * pruned_state_space_system % * pruned_state_space_system
% * vec % * vec
% ========================================================================= % =========================================================================
...@@ -123,19 +126,19 @@ useautocorr = options_ident.useautocorr; ...@@ -123,19 +126,19 @@ useautocorr = options_ident.useautocorr;
grid_nbr = options_ident.grid_nbr; grid_nbr = options_ident.grid_nbr;
kronflag = options_ident.analytic_derivation_mode; kronflag = options_ident.analytic_derivation_mode;
% get fields from M % get fields from M_
endo_nbr = M.endo_nbr; endo_nbr = M_.endo_nbr;
exo_nbr = M.exo_nbr; exo_nbr = M_.exo_nbr;
fname = M.fname; fname = M_.fname;
lead_lag_incidence = M.lead_lag_incidence; lead_lag_incidence = M_.lead_lag_incidence;
nspred = M.nspred; nspred = M_.nspred;
nstatic = M.nstatic; nstatic = M_.nstatic;
params = M.params; params = M_.params;
Sigma_e = M.Sigma_e; Sigma_e = M_.Sigma_e;
stderr_e = sqrt(diag(Sigma_e)); stderr_e = sqrt(diag(Sigma_e));
% set all selected values: stderr and corr come first, then model parameters % set all selected values: stderr and corr come first, then model parameters
xparam1 = get_all_parameters(estim_params, M); %try using estimated_params block xparam1 = get_all_parameters(estim_params, M_); %try using estimated_params block
if isempty(xparam1) if isempty(xparam1)
%if there is no estimated_params block, consider all stderr and all model parameters, but no corr parameters %if there is no estimated_params block, consider all stderr and all model parameters, but no corr parameters
xparam1 = [stderr_e', params']; xparam1 = [stderr_e', params'];
...@@ -150,94 +153,94 @@ obs_nbr = length(indvobs); ...@@ -150,94 +153,94 @@ obs_nbr = length(indvobs);
d2flag = 0; % do not compute second parameter derivatives d2flag = 0; % do not compute second parameter derivatives
% Get Jacobians (wrt selected params) of steady state, dynamic model derivatives and perturbation solution matrices for all endogenous variables % Get Jacobians (wrt selected params) of steady state, dynamic model derivatives and perturbation solution matrices for all endogenous variables
oo.dr.derivs = get_perturbation_params_derivs(M, options, estim_params, oo, indpmodel, indpstderr, indpcorr, d2flag); dr.derivs = identification.get_perturbation_params_derivs(M_, options_, estim_params, dr, endo_steady_state, exo_steady_state, exo_det_steady_state, indpmodel, indpstderr, indpcorr, d2flag);
[I,~] = find(lead_lag_incidence'); %I is used to select nonzero columns of the Jacobian of endogenous variables in dynamic model files [I,~] = find(lead_lag_incidence'); %I is used to select nonzero columns of the Jacobian of endogenous variables in dynamic model files
yy0 = oo.dr.ys(I); %steady state of dynamic (endogenous and auxiliary variables) in lead_lag_incidence order yy0 = dr.ys(I); %steady state of dynamic (endogenous and auxiliary variables) in lead_lag_incidence order
Yss = oo.dr.ys(oo.dr.order_var); % steady state in DR order Yss = dr.ys(dr.order_var); % steady state in DR order
if order == 1 if order == 1
[~, g1 ] = feval([fname,'.dynamic'], yy0, oo.exo_steady_state', params, oo.dr.ys, 1); [~, g1 ] = feval([fname,'.dynamic'], yy0, exo_steady_state', params, dr.ys, 1);
%g1 is [endo_nbr by yy0ex0_nbr first derivative (wrt all dynamic variables) of dynamic model equations, i.e. df/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order %g1 is [endo_nbr by yy0ex0_nbr first derivative (wrt all dynamic variables) of dynamic model equations, i.e. df/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order
DYNAMIC = [Yss; DYNAMIC = [Yss;
vec(g1)]; %add steady state and put rows of g1 in DR order vec(g1)]; %add steady state and put rows of g1 in DR order
dDYNAMIC = [oo.dr.derivs.dYss; dDYNAMIC = [dr.derivs.dYss;
reshape(oo.dr.derivs.dg1,size(oo.dr.derivs.dg1,1)*size(oo.dr.derivs.dg1,2),size(oo.dr.derivs.dg1,3)) ]; %reshape dg1 in DR order and add steady state reshape(dr.derivs.dg1,size(dr.derivs.dg1,1)*size(dr.derivs.dg1,2),size(dr.derivs.dg1,3)) ]; %reshape dg1 in DR order and add steady state
REDUCEDFORM = [Yss; REDUCEDFORM = [Yss;
vec(oo.dr.ghx); vec(dr.ghx);
dyn_vech(oo.dr.ghu*Sigma_e*transpose(oo.dr.ghu))]; %in DR order dyn_vech(dr.ghu*Sigma_e*transpose(dr.ghu))]; %in DR order
dREDUCEDFORM = zeros(endo_nbr*nspred+endo_nbr*(endo_nbr+1)/2, totparam_nbr); dREDUCEDFORM = zeros(endo_nbr*nspred+endo_nbr*(endo_nbr+1)/2, totparam_nbr);
for j=1:totparam_nbr for j=1:totparam_nbr
dREDUCEDFORM(:,j) = [vec(oo.dr.derivs.dghx(:,:,j)); dREDUCEDFORM(:,j) = [vec(dr.derivs.dghx(:,:,j));
dyn_vech(oo.dr.derivs.dOm(:,:,j))]; dyn_vech(dr.derivs.dOm(:,:,j))];
end end
dREDUCEDFORM = [ [zeros(endo_nbr, stderrparam_nbr+corrparam_nbr) oo.dr.derivs.dYss]; dREDUCEDFORM ]; % add steady state dREDUCEDFORM = [ [zeros(endo_nbr, stderrparam_nbr+corrparam_nbr) dr.derivs.dYss]; dREDUCEDFORM ]; % add steady state
elseif order == 2 elseif order == 2
[~, g1, g2 ] = feval([fname,'.dynamic'], yy0, oo.exo_steady_state', params, oo.dr.ys, 1); [~, g1, g2 ] = feval([fname,'.dynamic'], yy0, exo_steady_state', params, dr.ys, 1);
%g1 is [endo_nbr by yy0ex0_nbr first derivative (wrt all dynamic variables) of dynamic model equations, i.e. df/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order %g1 is [endo_nbr by yy0ex0_nbr first derivative (wrt all dynamic variables) of dynamic model equations, i.e. df/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order
%g2 is [endo_nbr by yy0ex0_nbr^2] second derivative (wrt all dynamic variables) of dynamic model equations, i.e. d(df/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order %g2 is [endo_nbr by yy0ex0_nbr^2] second derivative (wrt all dynamic variables) of dynamic model equations, i.e. d(df/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order
DYNAMIC = [Yss; DYNAMIC = [Yss;
vec(g1); vec(g1);
vec(g2)]; %add steady state and put rows of g1 and g2 in DR order vec(g2)]; %add steady state and put rows of g1 and g2 in DR order
dDYNAMIC = [oo.dr.derivs.dYss; dDYNAMIC = [dr.derivs.dYss;
reshape(oo.dr.derivs.dg1,size(oo.dr.derivs.dg1,1)*size(oo.dr.derivs.dg1,2),size(oo.dr.derivs.dg1,3)); %reshape dg1 in DR order reshape(dr.derivs.dg1,size(dr.derivs.dg1,1)*size(dr.derivs.dg1,2),size(dr.derivs.dg1,3)); %reshape dg1 in DR order
reshape(oo.dr.derivs.dg2,size(oo.dr.derivs.dg1,1)*size(oo.dr.derivs.dg1,2)^2,size(oo.dr.derivs.dg1,3))]; %reshape dg2 in DR order reshape(dr.derivs.dg2,size(dr.derivs.dg1,1)*size(dr.derivs.dg1,2)^2,size(dr.derivs.dg1,3))]; %reshape dg2 in DR order
REDUCEDFORM = [Yss; REDUCEDFORM = [Yss;
vec(oo.dr.ghx); vec(dr.ghx);
dyn_vech(oo.dr.ghu*Sigma_e*transpose(oo.dr.ghu)); dyn_vech(dr.ghu*Sigma_e*transpose(dr.ghu));
vec(oo.dr.ghxx); vec(dr.ghxx);
vec(oo.dr.ghxu); vec(dr.ghxu);
vec(oo.dr.ghuu); vec(dr.ghuu);
vec(oo.dr.ghs2)]; %in DR order vec(dr.ghs2)]; %in DR order
dREDUCEDFORM = zeros(endo_nbr*nspred+endo_nbr*(endo_nbr+1)/2+endo_nbr*nspred^2+endo_nbr*nspred*exo_nbr+endo_nbr*exo_nbr^2+endo_nbr, totparam_nbr); dREDUCEDFORM = zeros(endo_nbr*nspred+endo_nbr*(endo_nbr+1)/2+endo_nbr*nspred^2+endo_nbr*nspred*exo_nbr+endo_nbr*exo_nbr^2+endo_nbr, totparam_nbr);
for j=1:totparam_nbr for j=1:totparam_nbr
dREDUCEDFORM(:,j) = [vec(oo.dr.derivs.dghx(:,:,j)); dREDUCEDFORM(:,j) = [vec(dr.derivs.dghx(:,:,j));
dyn_vech(oo.dr.derivs.dOm(:,:,j)); dyn_vech(dr.derivs.dOm(:,:,j));
vec(oo.dr.derivs.dghxx(:,:,j)); vec(dr.derivs.dghxx(:,:,j));
vec(oo.dr.derivs.dghxu(:,:,j)); vec(dr.derivs.dghxu(:,:,j));
vec(oo.dr.derivs.dghuu(:,:,j)); vec(dr.derivs.dghuu(:,:,j));
vec(oo.dr.derivs.dghs2(:,j))]; vec(dr.derivs.dghs2(:,j))];
end end
dREDUCEDFORM = [ [zeros(endo_nbr, stderrparam_nbr+corrparam_nbr) oo.dr.derivs.dYss]; dREDUCEDFORM ]; % add steady state dREDUCEDFORM = [ [zeros(endo_nbr, stderrparam_nbr+corrparam_nbr) dr.derivs.dYss]; dREDUCEDFORM ]; % add steady state
elseif order == 3 elseif order == 3
[~, g1, g2, g3 ] = feval([fname,'.dynamic'], yy0, oo.exo_steady_state', params, oo.dr.ys, 1); [~, g1, g2, g3 ] = feval([fname,'.dynamic'], yy0, exo_steady_state', params, dr.ys, 1);
%g1 is [endo_nbr by yy0ex0_nbr first derivative (wrt all dynamic variables) of dynamic model equations, i.e. df/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order %g1 is [endo_nbr by yy0ex0_nbr first derivative (wrt all dynamic variables) of dynamic model equations, i.e. df/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order
%g2 is [endo_nbr by yy0ex0_nbr^2] second derivative (wrt all dynamic variables) of dynamic model equations, i.e. d(df/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order %g2 is [endo_nbr by yy0ex0_nbr^2] second derivative (wrt all dynamic variables) of dynamic model equations, i.e. d(df/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order
DYNAMIC = [Yss; DYNAMIC = [Yss;
vec(g1); vec(g1);
vec(g2); vec(g2);
vec(g3)]; %add steady state and put rows of g1 and g2 in DR order vec(g3)]; %add steady state and put rows of g1 and g2 in DR order
dDYNAMIC = [oo.dr.derivs.dYss; dDYNAMIC = [dr.derivs.dYss;
reshape(oo.dr.derivs.dg1,size(oo.dr.derivs.dg1,1)*size(oo.dr.derivs.dg1,2),size(oo.dr.derivs.dg1,3)); %reshape dg1 in DR order reshape(dr.derivs.dg1,size(dr.derivs.dg1,1)*size(dr.derivs.dg1,2),size(dr.derivs.dg1,3)); %reshape dg1 in DR order
reshape(oo.dr.derivs.dg2,size(oo.dr.derivs.dg1,1)*size(oo.dr.derivs.dg1,2)^2,size(oo.dr.derivs.dg1,3)); reshape(dr.derivs.dg2,size(dr.derivs.dg1,1)*size(dr.derivs.dg1,2)^2,size(dr.derivs.dg1,3));
reshape(oo.dr.derivs.dg2,size(oo.dr.derivs.dg1,1)*size(oo.dr.derivs.dg1,2)^2,size(oo.dr.derivs.dg1,3))]; %reshape dg3 in DR order reshape(dr.derivs.dg2,size(dr.derivs.dg1,1)*size(dr.derivs.dg1,2)^2,size(dr.derivs.dg1,3))]; %reshape dg3 in DR order
REDUCEDFORM = [Yss; REDUCEDFORM = [Yss;
vec(oo.dr.ghx); vec(dr.ghx);
dyn_vech(oo.dr.ghu*Sigma_e*transpose(oo.dr.ghu)); dyn_vech(dr.ghu*Sigma_e*transpose(dr.ghu));
vec(oo.dr.ghxx); vec(oo.dr.ghxu); vec(oo.dr.ghuu); vec(oo.dr.ghs2); vec(dr.ghxx); vec(dr.ghxu); vec(dr.ghuu); vec(dr.ghs2);
vec(oo.dr.ghxxx); vec(oo.dr.ghxxu); vec(oo.dr.ghxuu); vec(oo.dr.ghuuu); vec(oo.dr.ghxss); vec(oo.dr.ghuss)]; %in DR order vec(dr.ghxxx); vec(dr.ghxxu); vec(dr.ghxuu); vec(dr.ghuuu); vec(dr.ghxss); vec(dr.ghuss)]; %in DR order
dREDUCEDFORM = zeros(size(REDUCEDFORM,1)-endo_nbr, totparam_nbr); dREDUCEDFORM = zeros(size(REDUCEDFORM,1)-endo_nbr, totparam_nbr);
for j=1:totparam_nbr for j=1:totparam_nbr
dREDUCEDFORM(:,j) = [vec(oo.dr.derivs.dghx(:,:,j)); dREDUCEDFORM(:,j) = [vec(dr.derivs.dghx(:,:,j));
dyn_vech(oo.dr.derivs.dOm(:,:,j)); dyn_vech(dr.derivs.dOm(:,:,j));
vec(oo.dr.derivs.dghxx(:,:,j)); vec(oo.dr.derivs.dghxu(:,:,j)); vec(oo.dr.derivs.dghuu(:,:,j)); vec(oo.dr.derivs.dghs2(:,j)) vec(dr.derivs.dghxx(:,:,j)); vec(dr.derivs.dghxu(:,:,j)); vec(dr.derivs.dghuu(:,:,j)); vec(dr.derivs.dghs2(:,j))
vec(oo.dr.derivs.dghxxx(:,:,j)); vec(oo.dr.derivs.dghxxu(:,:,j)); vec(oo.dr.derivs.dghxuu(:,:,j)); vec(oo.dr.derivs.dghuuu(:,:,j)); vec(oo.dr.derivs.dghxss(:,:,j)); vec(oo.dr.derivs.dghuss(:,:,j))]; vec(dr.derivs.dghxxx(:,:,j)); vec(dr.derivs.dghxxu(:,:,j)); vec(dr.derivs.dghxuu(:,:,j)); vec(dr.derivs.dghuuu(:,:,j)); vec(dr.derivs.dghxss(:,:,j)); vec(dr.derivs.dghuss(:,:,j))];
end end
dREDUCEDFORM = [ [zeros(endo_nbr, stderrparam_nbr+corrparam_nbr) oo.dr.derivs.dYss]; dREDUCEDFORM ]; % add steady state dREDUCEDFORM = [ [zeros(endo_nbr, stderrparam_nbr+corrparam_nbr) dr.derivs.dYss]; dREDUCEDFORM ]; % add steady state
end end
% Get (pruned) state space representation: % Get (pruned) state space representation:
pruned = pruned_state_space_system(M, options, oo.dr, indvobs, nlags, useautocorr, 1); pruned = pruned_SS.pruned_state_space_system(M_, options_, dr, indvobs, nlags, useautocorr, 1);
MEAN = pruned.E_y; MEAN = pruned.E_y;
dMEAN = pruned.dE_y; dMEAN = pruned.dE_y;
%storage for Jacobians used in dsge_likelihood.m for analytical Gradient and Hession of likelihood (only at order=1) %storage for Jacobians used in dsge_likelihood.m for analytical Gradient and Hession of likelihood (only at order=1)
derivatives_info = struct(); derivatives_info = struct();
if order == 1 if order == 1
dT = zeros(endo_nbr,endo_nbr,totparam_nbr); dT = zeros(endo_nbr,endo_nbr,totparam_nbr);
dT(:,(nstatic+1):(nstatic+nspred),:) = oo.dr.derivs.dghx; dT(:,(nstatic+1):(nstatic+nspred),:) = dr.derivs.dghx;
derivatives_info.DT = dT; derivatives_info.DT = dT;
derivatives_info.DOm = oo.dr.derivs.dOm; derivatives_info.DOm = dr.derivs.dOm;
derivatives_info.DYss = oo.dr.derivs.dYss; derivatives_info.DYss = dr.derivs.dYss;
end end
%% Compute dMOMENTS %% Compute dMOMENTS
...@@ -255,7 +258,7 @@ if ~no_identification_moments ...@@ -255,7 +258,7 @@ if ~no_identification_moments
if kronflag == -1 if kronflag == -1
%numerical derivative of autocovariogram %numerical derivative of autocovariogram
dMOMENTS = fjaco(str2func('identification_numerical_objective'), xparam1, 1, estim_params, M, oo, options, indpmodel, indpstderr, indpcorr, indvobs, useautocorr, nlags, grid_nbr); %[outputflag=1] dMOMENTS = identification.fjaco(str2func('identification.numerical_objective'), xparam1, 1, estim_params, M_, options_, indpmodel, indpstderr, indvobs, useautocorr, nlags, grid_nbr, dr, endo_steady_state, exo_steady_state, exo_det_steady_state); %[outputflag=1]
dMOMENTS = [dMEAN; dMOMENTS]; %add Jacobian of steady state of VAROBS variables dMOMENTS = [dMEAN; dMOMENTS]; %add Jacobian of steady state of VAROBS variables
else else
dMOMENTS = zeros(obs_nbr + obs_nbr*(obs_nbr+1)/2 + nlags*obs_nbr^2 , totparam_nbr); dMOMENTS = zeros(obs_nbr + obs_nbr*(obs_nbr+1)/2 + nlags*obs_nbr^2 , totparam_nbr);
...@@ -312,7 +315,7 @@ if ~no_identification_spectrum ...@@ -312,7 +315,7 @@ if ~no_identification_spectrum
IA = eye(size(pruned.A,1)); IA = eye(size(pruned.A,1));
if kronflag == -1 if kronflag == -1
%numerical derivative of spectral density %numerical derivative of spectral density
dOmega_tmp = fjaco(str2func('identification_numerical_objective'), xparam1, 2, estim_params, M, oo, options, indpmodel, indpstderr, indpcorr, indvobs, useautocorr, nlags, grid_nbr); %[outputflag=2] dOmega_tmp = identification.fjaco(str2func('identification.numerical_objective'), xparam1, 2, estim_params, M_, options_, indpmodel, indpstderr, indvobs, useautocorr, nlags, grid_nbr, dr, endo_steady_state, exo_steady_state, exo_det_steady_state); %[outputflag=2]
kk = 0; kk = 0;
for ig = 1:length(freqs) for ig = 1:length(freqs)
kk = kk+1; kk = kk+1;
...@@ -330,7 +333,7 @@ if ~no_identification_spectrum ...@@ -330,7 +333,7 @@ if ~no_identification_spectrum
dC = reshape(pruned.dC,size(pruned.dC,1)*size(pruned.dC,2),size(pruned.dC,3)); dC = reshape(pruned.dC,size(pruned.dC,1)*size(pruned.dC,2),size(pruned.dC,3));
dD = reshape(pruned.dD,size(pruned.dD,1)*size(pruned.dD,2),size(pruned.dD,3)); dD = reshape(pruned.dD,size(pruned.dD,1)*size(pruned.dD,2),size(pruned.dD,3));
dVarinov = reshape(pruned.dVarinov,size(pruned.dVarinov,1)*size(pruned.dVarinov,2),size(pruned.dVarinov,3)); dVarinov = reshape(pruned.dVarinov,size(pruned.dVarinov,1)*size(pruned.dVarinov,2),size(pruned.dVarinov,3));
K_obs_exo = commutation(obs_nbr,size(pruned.Varinov,1)); K_obs_exo = pruned_SS.commutation(obs_nbr,size(pruned.Varinov,1));
for ig=1:length(freqs) for ig=1:length(freqs)
z = tneg(ig); z = tneg(ig);
zIminusA = (z*IA - pruned.A); zIminusA = (z*IA - pruned.A);
...@@ -389,15 +392,15 @@ if ~no_identification_minimal ...@@ -389,15 +392,15 @@ if ~no_identification_minimal
dMINIMAL = []; dMINIMAL = [];
else else
% Derive and check minimal state vector of first-order % Derive and check minimal state vector of first-order
SYS.A = oo.dr.ghx(pruned.indx,:); SYS.A = dr.ghx(pruned.indx,:);
SYS.dA = oo.dr.derivs.dghx(pruned.indx,:,:); SYS.dA = dr.derivs.dghx(pruned.indx,:,:);
SYS.B = oo.dr.ghu(pruned.indx,:); SYS.B = dr.ghu(pruned.indx,:);
SYS.dB = oo.dr.derivs.dghu(pruned.indx,:,:); SYS.dB = dr.derivs.dghu(pruned.indx,:,:);
SYS.C = oo.dr.ghx(pruned.indy,:); SYS.C = dr.ghx(pruned.indy,:);
SYS.dC = oo.dr.derivs.dghx(pruned.indy,:,:); SYS.dC = dr.derivs.dghx(pruned.indy,:,:);
SYS.D = oo.dr.ghu(pruned.indy,:); SYS.D = dr.ghu(pruned.indy,:);
SYS.dD = oo.dr.derivs.dghu(pruned.indy,:,:); SYS.dD = dr.derivs.dghu(pruned.indy,:,:);
[CheckCO,minnx,SYS] = get_minimal_state_representation(SYS,1); [CheckCO,minnx,SYS] = identification.get_minimal_state_representation(SYS,1);
if CheckCO == 0 if CheckCO == 0
warning_KomunjerNg = 'WARNING: Komunjer and Ng (2011) failed:\n'; warning_KomunjerNg = 'WARNING: Komunjer and Ng (2011) failed:\n';
...@@ -415,12 +418,12 @@ if ~no_identification_minimal ...@@ -415,12 +418,12 @@ if ~no_identification_minimal
dminB = reshape(dminB,size(dminB,1)*size(dminB,2),size(dminB,3)); dminB = reshape(dminB,size(dminB,1)*size(dminB,2),size(dminB,3));
dminC = reshape(dminC,size(dminC,1)*size(dminC,2),size(dminC,3)); dminC = reshape(dminC,size(dminC,1)*size(dminC,2),size(dminC,3));
dminD = reshape(dminD,size(dminD,1)*size(dminD,2),size(dminD,3)); dminD = reshape(dminD,size(dminD,1)*size(dminD,2),size(dminD,3));
dvechSig = reshape(oo.dr.derivs.dSigma_e,exo_nbr*exo_nbr,totparam_nbr); dvechSig = reshape(dr.derivs.dSigma_e,exo_nbr*exo_nbr,totparam_nbr);
indvechSig= find(tril(ones(exo_nbr,exo_nbr))); indvechSig= find(tril(ones(exo_nbr,exo_nbr)));
dvechSig = dvechSig(indvechSig,:); dvechSig = dvechSig(indvechSig,:);
Inx = eye(minnx); Inx = eye(minnx);
Inu = eye(exo_nbr); Inu = eye(exo_nbr);
[~,Enu] = duplication(exo_nbr); [~,Enu] = pruned_SS.duplication(exo_nbr);
KomunjerNg_DL = [dminA; dminB; dminC; dminD; dvechSig]; KomunjerNg_DL = [dminA; dminB; dminC; dminD; dvechSig];
KomunjerNg_DT = [kron(transpose(minA),Inx) - kron(Inx,minA); KomunjerNg_DT = [kron(transpose(minA),Inx) - kron(Inx,minA);
kron(transpose(minB),Inx); kron(transpose(minB),Inx);
......
matlab/get_minimal_state_representation.m → matlab/+identification/get_minimal_state_representation.m
View file @ 4a32acc2
...@@ -53,7 +53,7 @@ function [CheckCO,minns,minSYS] = get_minimal_state_representation(SYS, derivs_f ...@@ -53,7 +53,7 @@ function [CheckCO,minns,minSYS] = get_minimal_state_representation(SYS, derivs_f
% Jacobian (wrt to all parameters) of measurement matrix minD % Jacobian (wrt to all parameters) of measurement matrix minD
% ------------------------------------------------------------------------- % -------------------------------------------------------------------------
% This function is called by % This function is called by
% * get_identification_jacobians.m (previously getJJ.m) % * identification.get_jacobians.m (previously getJJ.m)
% ------------------------------------------------------------------------- % -------------------------------------------------------------------------
% This function calls % This function calls
% * check_minimality (embedded) % * check_minimality (embedded)
...@@ -209,7 +209,7 @@ function [mSYS,U] = minrealold(SYS,tol) ...@@ -209,7 +209,7 @@ function [mSYS,U] = minrealold(SYS,tol)
a = SYS.A; a = SYS.A;
b = SYS.B; b = SYS.B;
c = SYS.C; c = SYS.C;
[ns,nu] = size(b); [ns,~] = size(b);
[am,bm,cm,U,k] = ControllabilityStaircaseRosenbrock(a,b,c,tol); [am,bm,cm,U,k] = ControllabilityStaircaseRosenbrock(a,b,c,tol);
kk = sum(k); kk = sum(k);
nu = ns - kk; nu = ns - kk;
...@@ -219,7 +219,7 @@ function [mSYS,U] = minrealold(SYS,tol) ...@@ -219,7 +219,7 @@ function [mSYS,U] = minrealold(SYS,tol)
cm = cm(:,nu+1:ns); cm = cm(:,nu+1:ns);
ns = ns - nu; ns = ns - nu;
if ns if ns
[am,bm,cm,t,k] = ObservabilityStaircaseRosenbrock(am,bm,cm,tol); [am,bm,cm,~,k] = ObservabilityStaircaseRosenbrock(am,bm,cm,tol);
kk = sum(k); kk = sum(k);
nu = ns - kk; nu = ns - kk;
nn = nn + nu; nn = nn + nu;
...@@ -242,8 +242,8 @@ end ...@@ -242,8 +242,8 @@ end
function [abar,bbar,cbar,t,k] = ControllabilityStaircaseRosenbrock(a, b, c, tol) function [abar,bbar,cbar,t,k] = ControllabilityStaircaseRosenbrock(a, b, c, tol)
% Controllability staircase algorithm of Rosenbrock, 1968 % Controllability staircase algorithm of Rosenbrock, 1968
[ra,ca] = size(a); [ra,~] = size(a);
[rb,cb] = size(b); [~,cb] = size(b);
ptjn1 = eye(ra); ptjn1 = eye(ra);
ajn1 = a; ajn1 = a;
bjn1 = b; bjn1 = b;
...@@ -255,8 +255,8 @@ function [abar,bbar,cbar,t,k] = ControllabilityStaircaseRosenbrock(a, b, c, tol) ...@@ -255,8 +255,8 @@ function [abar,bbar,cbar,t,k] = ControllabilityStaircaseRosenbrock(a, b, c, tol)
tol = ra*norm(a,1)*eps; tol = ra*norm(a,1)*eps;
end end
for jj = 1 : ra for jj = 1 : ra
[uj,sj,vj] = svd(bjn1); [uj,sj] = svd(bjn1);
[rsj,csj] = size(sj); [rsj,~] = size(sj);
%p = flip(eye(rsj),2); %p = flip(eye(rsj),2);
p = eye(rsj); p = eye(rsj);
p = p(:,end:-1:1); p = p(:,end:-1:1);
...@@ -264,7 +264,7 @@ function [abar,bbar,cbar,t,k] = ControllabilityStaircaseRosenbrock(a, b, c, tol) ...@@ -264,7 +264,7 @@ function [abar,bbar,cbar,t,k] = ControllabilityStaircaseRosenbrock(a, b, c, tol)
uj = uj*p; uj = uj*p;
bb = uj'*bjn1; bb = uj'*bjn1;
roj = rank(bb,tol); roj = rank(bb,tol);
[rbb,cbb] = size(bb); [rbb,~] = size(bb);
sigmaj = rbb - roj; sigmaj = rbb - roj;
sigmajn1 = sigmaj; sigmajn1 = sigmaj;
k(jj) = roj; k(jj) = roj;
......
matlab/get_perturbation_params_derivs.m → matlab/+identification/get_perturbation_params_derivs.m
View file @ 4a32acc2
function DERIVS = get_perturbation_params_derivs(M, options, estim_params, oo, 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, oo, 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
...@@ -16,24 +16,31 @@ function DERIVS = get_perturbation_params_derivs(M, options, estim_params, oo, i ...@@ -16,24 +16,31 @@ function DERIVS = get_perturbation_params_derivs(M, options, estim_params, oo, i
% %
% ========================================================================= % =========================================================================
% INPUTS % INPUTS
% M: [structure] storing the model information % M_: [structure] storing the model information
% options: [structure] storing the options % options_: [structure] storing the options
% estim_params: [structure] storing the estimation information % estim_params_: [structure] storing the estimation information
% oo: [structure] storing the solution results % dr [structure] Reduced form model.
% indpmodel: [modparam_nbr by 1] index of selected (estimated) parameters in M.params; % endo_steady_state [vector] steady state value for endogenous variables
% corresponds to model parameters (no stderr and no corr) in estimated_params block % exo_steady_state [vector] steady state value for exogenous variables
% exo_det_steady_state [vector] steady state value for exogenous deterministic variables
% indpmodel: [modparam_nbr by 1] index of selected (estimated) parameters in M_.params;
% corresponds to model parameters (no stderr and no corr)
% in estimated_params block
% indpstderr: [stderrparam_nbr by 1] index of selected (estimated) standard errors, % indpstderr: [stderrparam_nbr by 1] index of selected (estimated) standard errors,
% i.e. for all exogenous variables where 'stderr' is given in the estimated_params block % i.e. for all exogenous variables where 'stderr' is given
% in the estimated_params block
% indpcorr: [corrparam_nbr by 2] matrix of selected (estimated) correlations, % indpcorr: [corrparam_nbr by 2] matrix of selected (estimated) correlations,
% i.e. for all exogenous variables where 'corr' is given in the estimated_params block % i.e. for all exogenous variables where 'corr' is given in
% d2flag: [boolean] flag to compute second-order parameter derivatives of steady state and first-order Kalman transition matrices % the estimated_params block
% d2flag: [boolean] flag to compute second-order parameter derivatives of steady state
% and first-order Kalman transition matrices
% ------------------------------------------------------------------------- % -------------------------------------------------------------------------
% OUTPUTS % OUTPUTS
% DERIVS: Structure with the following fields: % DERIVS: Structure with the following fields:
% 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
...@@ -49,7 +56,7 @@ function DERIVS = get_perturbation_params_derivs(M, options, estim_params, oo, i ...@@ -49,7 +56,7 @@ function DERIVS = get_perturbation_params_derivs(M, options, estim_params, oo, i
% 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
...@@ -81,12 +88,13 @@ function DERIVS = get_perturbation_params_derivs(M, options, estim_params, oo, i ...@@ -81,12 +88,13 @@ function DERIVS = get_perturbation_params_derivs(M, options, estim_params, oo, i
% ------------------------------------------------------------------------- % -------------------------------------------------------------------------
% This function is called by % This function is called by
% * dsge_likelihood.m % * dsge_likelihood.m
% * get_identification_jacobians.m % * identification.get_jacobians.m
% ------------------------------------------------------------------------- % -------------------------------------------------------------------------
% This function calls % This function calls
% * [fname,'.dynamic'] % * [fname,'.dynamic']
% * [fname,'.dynamic_params_derivs'] % * [fname,'.dynamic_params_derivs']
% * [fname,'.static'] % * [fname,'.sparse.static_g1']
% * [fname,'.sparse.static_g2']
% * [fname,'.static_params_derivs'] % * [fname,'.static_params_derivs']
% * commutation % * commutation
% * dyn_vech % * dyn_vech
...@@ -102,7 +110,7 @@ function DERIVS = get_perturbation_params_derivs(M, options, estim_params, oo, i ...@@ -102,7 +110,7 @@ function DERIVS = get_perturbation_params_derivs(M, options, estim_params, oo, i
% * sylvester3a % * sylvester3a
% * get_perturbation_params_derivs_numerical_objective % * get_perturbation_params_derivs_numerical_objective
% ========================================================================= % =========================================================================
% Copyright © 2019-2020 Dynare Team % Copyright © 2019-2024 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
...@@ -119,61 +127,60 @@ function DERIVS = get_perturbation_params_derivs(M, options, estim_params, oo, i ...@@ -119,61 +127,60 @@ function DERIVS = get_perturbation_params_derivs(M, options, estim_params, oo, i
% 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/>.
% ========================================================================= % =========================================================================
% Get fields from M % Get fields from M_
Correlation_matrix = M.Correlation_matrix; Correlation_matrix = M_.Correlation_matrix;
dname = M.dname; dname = M_.dname;
dynamic_tmp_nbr = M.dynamic_tmp_nbr; dynamic_tmp_nbr = M_.dynamic_tmp_nbr;
endo_nbr = M.endo_nbr; endo_nbr = M_.endo_nbr;
exo_nbr = M.exo_nbr; exo_nbr = M_.exo_nbr;
exo_det_nbr = M.exo_det_nbr; exo_det_nbr = M_.exo_det_nbr;
fname = M.fname; fname = M_.fname;
lead_lag_incidence = M.lead_lag_incidence; lead_lag_incidence = M_.lead_lag_incidence;
nfwrd = M.nfwrd; nfwrd = M_.nfwrd;
npred = M.npred; npred = M_.npred;
nspred = M.nspred; nspred = M_.nspred;
nstatic = M.nstatic; nstatic = M_.nstatic;
params = M.params; params = M_.params;
param_nbr = M.param_nbr; param_nbr = M_.param_nbr;
Sigma_e = M.Sigma_e; Sigma_e = M_.Sigma_e;
% Get fields from options % Get fields from options_
analytic_derivation_mode = options.analytic_derivation_mode; analytic_derivation_mode = options_.analytic_derivation_mode;
% analytic_derivation_mode: select method to compute Jacobians, default is 0 % analytic_derivation_mode: select method to compute Jacobians, default is 0
% * 0: efficient sylvester equation method to compute analytical derivatives as in Ratto & Iskrev (2012) % * 0: efficient sylvester equation method to compute analytical derivatives as in Ratto & Iskrev (2012)
% * 1: kronecker products method to compute analytical derivatives as in Iskrev (2010), only for order=1 % * 1: kronecker products method to compute analytical derivatives as in Iskrev (2010), only for order=1
% * -1: numerical two-sided finite difference method to compute numerical derivatives of all output arguments using function get_perturbation_params_derivs_numerical_objective.m % * -1: numerical two-sided finite difference method to compute numerical derivatives of all output arguments using function get_perturbation_params_derivs_numerical_objective.m
% * -2: numerical two-sided finite difference method to compute numerically dYss, dg1, dg2, dg3, d2Yss and d2g1, the other output arguments are computed analytically as in kronflag=0 % * -2: numerical two-sided finite difference method to compute numerically dYss, dg1, dg2, dg3, d2Yss and d2g1, the other output arguments are computed analytically as in kronflag=0
gstep = options.gstep; gstep = options_.gstep;
order = options.order; order = options_.order;
if isempty(options.qz_criterium) if isempty(options_.qz_criterium)
% set default value for qz_criterium: if there are no unit roots one can use 1.0 % set default value for qz_criterium: if there are no unit roots one can use 1.0
% If they are possible, you may have have multiple unit roots and the accuracy % If they are possible, you may have have multiple unit roots and the accuracy
% decreases when computing the eigenvalues in lyapunov_symm. Hence, we normally use 1+1e-6 % decreases when computing the eigenvalues in lyapunov_symm. Hence, we normally use 1+1e-6
options = select_qz_criterium_value(options); options_ = select_qz_criterium_value(options_);
end end
qz_criterium = options.qz_criterium; qz_criterium = options_.qz_criterium;
threads_BC = options.threads.kronecker.sparse_hessian_times_B_kronecker_C; threads_BC = options_.threads.kronecker.sparse_hessian_times_B_kronecker_C;
% Get fields from oo % Get fields from dr
exo_steady_state = oo.exo_steady_state; ghx = dr.ghx;
ghx = oo.dr.ghx; ghu = dr.ghu;
ghu = oo.dr.ghu;
if order > 1 if order > 1
ghxx = oo.dr.ghxx; ghxx = dr.ghxx;
ghxu = oo.dr.ghxu; ghxu = dr.ghxu;
ghuu = oo.dr.ghuu; ghuu = dr.ghuu;
ghs2 = oo.dr.ghs2; ghs2 = dr.ghs2;
end end
if order > 2 if order > 2
ghxxx = oo.dr.ghxxx; ghxxx = dr.ghxxx;
ghxxu = oo.dr.ghxxu; ghxxu = dr.ghxxu;
ghxuu = oo.dr.ghxuu; ghxuu = dr.ghxuu;
ghuuu = oo.dr.ghuuu; ghuuu = dr.ghuuu;
ghxss = oo.dr.ghxss; ghxss = dr.ghxss;
ghuss = oo.dr.ghuss; ghuss = dr.ghuss;
end end
order_var = oo.dr.order_var; order_var = dr.order_var;
ys = oo.dr.ys; ys = dr.ys;
% Some checks % Some checks
if exo_det_nbr > 0 if exo_det_nbr > 0
...@@ -181,11 +188,11 @@ if exo_det_nbr > 0 ...@@ -181,11 +188,11 @@ if exo_det_nbr > 0
end end
if order > 1 && analytic_derivation_mode == 1 if order > 1 && analytic_derivation_mode == 1
%analytic derivatives using Kronecker products is implemented only at first-order, at higher order we reset to analytic derivatives with sylvester equations %analytic derivatives using Kronecker products is implemented only at first-order, at higher order we reset to analytic derivatives with sylvester equations
%options.analytic_derivation_mode = 0; fprintf('As order > 1, reset ''analytic_derivation_mode'' to 0\n'); %options_.analytic_derivation_mode = 0; fprintf('As order > 1, reset ''analytic_derivation_mode'' to 0\n');
analytic_derivation_mode = 0; fprintf('As order > 1, reset ''analytic_derivation_mode'' to 0\n'); analytic_derivation_mode = 0; fprintf('As order > 1, reset ''analytic_derivation_mode'' to 0\n');
end end
numerical_objective_fname = str2func('get_perturbation_params_derivs_numerical_objective'); numerical_objective_fname = str2func('identification.get_perturbation_params_derivs_numerical_objective');
idx_states = nstatic+(1:nspred); %index for state variables, in DR order idx_states = nstatic+(1:nspred); %index for state variables, in DR order
modparam_nbr = length(indpmodel); %number of selected model parameters modparam_nbr = length(indpmodel); %number of selected model parameters
stderrparam_nbr = length(indpstderr); %number of selected stderr parameters stderrparam_nbr = length(indpstderr); %number of selected stderr parameters
...@@ -221,35 +228,35 @@ else %normal models ...@@ -221,35 +228,35 @@ else %normal models
end end
if analytic_derivation_mode < 0 if analytic_derivation_mode < 0
%Create auxiliary estim_params blocks if not available for numerical derivatives, estim_params_model contains only model parameters %Create auxiliary estim_params_ blocks if not available for numerical derivatives, estim_params_model contains only model parameters
estim_params_model.np = length(indpmodel); estim_params_model.np = length(indpmodel);
estim_params_model.param_vals(:,1) = indpmodel; estim_params_model.param_vals(:,1) = indpmodel;
estim_params_model.nvx = 0; estim_params_model.ncx = 0; estim_params_model.nvn = 0; estim_params_model.ncn = 0; estim_params_model.nvx = 0; estim_params_model.ncx = 0; estim_params_model.nvn = 0; estim_params_model.ncn = 0;
modparam1 = get_all_parameters(estim_params_model, M); %get all selected model parameters modparam1 = get_all_parameters(estim_params_model, M_); %get all selected model parameters
if ~isempty(indpstderr) && isempty(estim_params.var_exo) %if there are stderr parameters but no estimated_params_block if ~isempty(indpstderr) && isempty(estim_params_.var_exo) %if there are stderr parameters but no estimated_params_block
%provide temporary necessary information for stderr parameters %provide temporary necessary information for stderr parameters
estim_params.nvx = length(indpstderr); estim_params_.nvx = length(indpstderr);
estim_params.var_exo(:,1) = indpstderr; estim_params_.var_exo(:,1) = indpstderr;
end end
if ~isempty(indpcorr) && isempty(estim_params.corrx) %if there are corr parameters but no estimated_params_block if ~isempty(indpcorr) && isempty(estim_params_.corrx) %if there are corr parameters but no estimated_params_block
%provide temporary necessary information for stderr parameters %provide temporary necessary information for stderr parameters
estim_params.ncx = size(indpcorr,1); estim_params_.ncx = size(indpcorr,1);
estim_params.corrx(:,1:2) = indpcorr; estim_params_.corrx(:,1:2) = indpcorr;
end end
if ~isfield(estim_params,'nvn') %stderr of measurement errors not yet if ~isfield(estim_params_,'nvn') %stderr of measurement errors not yet
estim_params.nvn = 0; estim_params_.nvn = 0;
estim_params.var_endo = []; estim_params_.var_endo = [];
end end
if ~isfield(estim_params,'ncn') %corr of measurement errors not yet if ~isfield(estim_params_,'ncn') %corr of measurement errors not yet
estim_params.ncn = 0; estim_params_.ncn = 0;
estim_params.corrn = []; estim_params_.corrn = [];
end end
if ~isempty(indpmodel) && isempty(estim_params.param_vals) %if there are model parameters but no estimated_params_block if ~isempty(indpmodel) && isempty(estim_params_.param_vals) %if there are model parameters but no estimated_params_block
%provide temporary necessary information for model parameters %provide temporary necessary information for model parameters
estim_params.np = length(indpmodel); estim_params_.np = length(indpmodel);
estim_params.param_vals(:,1) = indpmodel; estim_params_.param_vals(:,1) = indpmodel;
end end
xparam1 = get_all_parameters(estim_params, M); %get all selected stderr, corr, and model parameters in estimated_params block, stderr and corr come first, then model parameters xparam1 = get_all_parameters(estim_params_, M_); %get all selected stderr, corr, and model parameters in estimated_params block, stderr and corr come first, then model parameters
end end
if d2flag if d2flag
modparam_nbr2 = modparam_nbr*(modparam_nbr+1)/2; %number of unique entries of selected model parameters only in second-order derivative matrix modparam_nbr2 = modparam_nbr*(modparam_nbr+1)/2; %number of unique entries of selected model parameters only in second-order derivative matrix
...@@ -288,8 +295,8 @@ if analytic_derivation_mode == -1 ...@@ -288,8 +295,8 @@ 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 = fjaco(numerical_objective_fname, xparam1, 'perturbation_solution', estim_params, M, oo, options); 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);
ind_ghu = ind_ghx(end) + (1:endo_nbr*exo_nbr); ind_ghu = ind_ghx(end) + (1:endo_nbr*exo_nbr);
...@@ -320,7 +327,7 @@ if analytic_derivation_mode == -1 ...@@ -320,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
...@@ -342,16 +349,16 @@ if analytic_derivation_mode == -1 ...@@ -342,16 +349,16 @@ if analytic_derivation_mode == -1
end end
%Parameter Jacobian of dynamic model derivatives (wrt selected model parameters only) %Parameter Jacobian of dynamic model derivatives (wrt selected model parameters only)
dYss_g = fjaco(numerical_objective_fname, modparam1, 'dynamic_model', estim_params_model, M, oo, options); dYss_g = identification.fjaco(numerical_objective_fname, modparam1, 'dynamic_model', estim_params_model, M_, options_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state);
ind_Yss = 1:endo_nbr; ind_Yss = 1:endo_nbr;
if options.discretionary_policy || options.ramsey_policy if options_.discretionary_policy || options_.ramsey_policy
ind_g1 = ind_Yss(end) + (1:M.eq_nbr*yy0ex0_nbr); ind_g1 = ind_Yss(end) + (1:M_.eq_nbr*yy0ex0_nbr);
else else
ind_g1 = ind_Yss(end) + (1:endo_nbr*yy0ex0_nbr); ind_g1 = ind_Yss(end) + (1:endo_nbr*yy0ex0_nbr);
end end
DERIVS.dYss = dYss_g(ind_Yss, :); %in tensor notation, wrt selected model parameters only DERIVS.dYss = dYss_g(ind_Yss, :); %in tensor notation, wrt selected model parameters only
if options.discretionary_policy || options.ramsey_policy if options_.discretionary_policy || options_.ramsey_policy
DERIVS.dg1 = reshape(dYss_g(ind_g1,:),[M.eq_nbr, yy0ex0_nbr, modparam_nbr]); %in tensor notation, wrt selected model parameters only DERIVS.dg1 = reshape(dYss_g(ind_g1,:),[M_.eq_nbr, yy0ex0_nbr, modparam_nbr]); %in tensor notation, wrt selected model parameters only
else else
DERIVS.dg1 = reshape(dYss_g(ind_g1,:),[endo_nbr, yy0ex0_nbr, modparam_nbr]); %in tensor notation, wrt selected model parameters only DERIVS.dg1 = reshape(dYss_g(ind_g1,:),[endo_nbr, yy0ex0_nbr, modparam_nbr]); %in tensor notation, wrt selected model parameters only
end end
...@@ -365,11 +372,11 @@ if analytic_derivation_mode == -1 ...@@ -365,11 +372,11 @@ 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 = hessian_sparse(numerical_objective_fname, xparam1, gstep, 'Kalman_Transition', estim_params, M, oo, options); 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);
options.order = order; %make sure to set back options_.order = order; %make sure to set back
ind_KalmanA = ind_Yss(end) + (1:endo_nbr^2); ind_KalmanA = ind_Yss(end) + (1:endo_nbr^2);
DERIVS.d2KalmanA = d2Yss_KalmanA_Om(ind_KalmanA, indp2tottot2); %only unique elements DERIVS.d2KalmanA = d2Yss_KalmanA_Om(ind_KalmanA, indp2tottot2); %only unique elements
DERIVS.d2Om = d2Yss_KalmanA_Om(ind_KalmanA(end)+1:end , indp2tottot2); %only unique elements DERIVS.d2Om = d2Yss_KalmanA_Om(ind_KalmanA(end)+1:end , indp2tottot2); %only unique elements
...@@ -388,7 +395,7 @@ if analytic_derivation_mode == -2 ...@@ -388,7 +395,7 @@ if analytic_derivation_mode == -2
% The parameter derivatives of perturbation solution matrices are computed analytically below (analytic_derivation_mode=0) % The parameter derivatives of perturbation solution matrices are computed analytically below (analytic_derivation_mode=0)
if order == 3 if order == 3
[~, g1, g2, g3] = feval([fname,'.dynamic'], ys(I), exo_steady_state', params, ys, 1); [~, g1, g2, g3] = feval([fname,'.dynamic'], ys(I), exo_steady_state', params, ys, 1);
g3 = unfold_g3(g3, yy0ex0_nbr); g3 = identification.unfold_g3(g3, yy0ex0_nbr);
elseif order == 2 elseif order == 2
[~, g1, g2] = feval([fname,'.dynamic'], ys(I), exo_steady_state', params, ys, 1); [~, g1, g2] = feval([fname,'.dynamic'], ys(I), exo_steady_state', params, ys, 1);
elseif order == 1 elseif order == 1
...@@ -398,9 +405,9 @@ if analytic_derivation_mode == -2 ...@@ -398,9 +405,9 @@ if analytic_derivation_mode == -2
if d2flag if d2flag
% computation of d2Yss and d2g1 % computation of d2Yss and d2g1
% note that hessian_sparse does not take symmetry into account, i.e. compare hessian_sparse.m to hessian.m, but focuses already on unique values, which are duplicated below % note that hessian_sparse does not take symmetry into account, i.e. compare hessian_sparse.m to hessian.m, but focuses already on unique values, which are duplicated below
options.order = 1; %d2flag requires only first order options_.order = 1; %d2flag requires only first order
d2Yss_g1 = hessian_sparse(numerical_objective_fname, modparam1, gstep, 'dynamic_model', estim_params_model, M, oo, options); % d2flag requires only first-order d2Yss_g1 = identification.hessian_sparse(numerical_objective_fname, modparam1, gstep, 'dynamic_model', estim_params_model, M_, options_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state); % d2flag requires only first-order
options.order = order; %make sure to set back the order options_.order = order; %make sure to set back the order
d2Yss = reshape(full(d2Yss_g1(1:endo_nbr,:)), [endo_nbr modparam_nbr modparam_nbr]); %put into tensor notation d2Yss = reshape(full(d2Yss_g1(1:endo_nbr,:)), [endo_nbr modparam_nbr modparam_nbr]); %put into tensor notation
for j=1:endo_nbr for j=1:endo_nbr
d2Yss(j,:,:) = dyn_unvech(dyn_vech(d2Yss(j,:,:))); %add duplicate values to full hessian d2Yss(j,:,:) = dyn_unvech(dyn_vech(d2Yss(j,:,:))); %add duplicate values to full hessian
...@@ -425,7 +432,7 @@ if analytic_derivation_mode == -2 ...@@ -425,7 +432,7 @@ if analytic_derivation_mode == -2
end end
%Parameter Jacobian of dynamic model derivatives (wrt selected model parameters only) %Parameter Jacobian of dynamic model derivatives (wrt selected model parameters only)
dYss_g = fjaco(numerical_objective_fname, modparam1, 'dynamic_model', estim_params_model, M, oo, options); dYss_g = identification.fjaco(numerical_objective_fname, modparam1, 'dynamic_model', estim_params_model, M_, options_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state);
ind_Yss = 1:endo_nbr; ind_Yss = 1:endo_nbr;
ind_g1 = ind_Yss(end) + (1:endo_nbr*yy0ex0_nbr); ind_g1 = ind_Yss(end) + (1:endo_nbr*yy0ex0_nbr);
dYss = dYss_g(ind_Yss,:); %in tensor notation, wrt selected model parameters only dYss = dYss_g(ind_Yss,:); %in tensor notation, wrt selected model parameters only
...@@ -441,20 +448,23 @@ if analytic_derivation_mode == -2 ...@@ -441,20 +448,23 @@ if analytic_derivation_mode == -2
clear dYss_g clear dYss_g
elseif (analytic_derivation_mode == 0 || analytic_derivation_mode == 1) elseif (analytic_derivation_mode == 0 || analytic_derivation_mode == 1)
%% Analytical computation of Jacobian and Hessian (wrt selected model parameters) of steady state, i.e. dYss and d2Yss if ~exist(['+' fname filesep 'static_params_derivs.m'],'file')
[~, g1_static] = feval([fname,'.static'], ys, exo_steady_state', params); %g1_static is [endo_nbr by endo_nbr] first-derivative (wrt all endogenous variables) of static model equations f, i.e. df/dys, in declaration order
try
rp_static = feval([fname,'.static_params_derivs'], ys, exo_steady_state', params); %rp_static is [endo_nbr by param_nbr] first-derivative (wrt all model parameters) of static model equations f, i.e. df/dparams, in declaration order
catch
error('For analytical parameter derivatives ''static_params_derivs.m'' file is needed, this can be created by putting identification(order=%d) into your mod file.',order) error('For analytical parameter derivatives ''static_params_derivs.m'' file is needed, this can be created by putting identification(order=%d) into your mod file.',order)
end end
if ~exist(['+' fname filesep 'dynamic_params_derivs.m'],'file')
error('For analytical parameter derivatives ''dynamic_params_derivs.m'' file is needed, this can be created by putting identification(order=%d) into your mod file.',order)
end
%% Analytical computation of Jacobian and Hessian (wrt selected model parameters) of steady state, i.e. dYss and d2Yss
[g1_static, T_order_static, T_static] = feval([fname,'.sparse.static_g1'], ys, exo_steady_state', params, M_.static_g1_sparse_rowval, M_.static_g1_sparse_colval, M_.static_g1_sparse_colptr); %g1_static is [endo_nbr by endo_nbr] first-derivative (wrt all endogenous variables) of static model equations f, i.e. df/dys, in declaration order
rp_static = feval([fname,'.static_params_derivs'], ys, exo_steady_state', params); %rp_static is [endo_nbr by param_nbr] first-derivative (wrt all model parameters) of static model equations f, i.e. df/dparams, in declaration order
dys = -g1_static\rp_static; %use implicit function theorem (equation 5 of Ratto and Iskrev (2012) to compute [endo_nbr by param_nbr] first-derivative (wrt all model parameters) of steady state for all endogenous variables analytically, note that dys is in declaration order dys = -g1_static\rp_static; %use implicit function theorem (equation 5 of Ratto and Iskrev (2012) to compute [endo_nbr by param_nbr] first-derivative (wrt all model parameters) of steady state for all endogenous variables analytically, note that dys is in declaration order
d2ys = zeros(endo_nbr, param_nbr, param_nbr); %initialize in tensor notation, note that d2ys is only needed for d2flag, i.e. for g1pp d2ys = zeros(endo_nbr, param_nbr, param_nbr); %initialize in tensor notation, note that d2ys is only needed for d2flag, i.e. for g1pp
if d2flag if d2flag
[~, ~, g2_static] = feval([fname,'.static'], ys, exo_steady_state', params); %g2_static is [endo_nbr by endo_nbr^2] second derivative (wrt all endogenous variables) of static model equations f, i.e. d(df/dys)/dys, in declaration order g2_static_v = feval([fname,'.sparse.static_g2'], ys, exo_steady_state', params, T_order_static, T_static);
g2_static = build_two_dim_hessian(M_.static_g2_sparse_indices, g2_static_v, endo_nbr, endo_nbr); %g2_static is [endo_nbr by endo_nbr^2] second derivative (wrt all endogenous variables) of static model equations f, i.e. d(df/dys)/dys, in declaration order
if order < 3 if order < 3
[~, g1, g2, g3] = feval([fname,'.dynamic'], ys(I), exo_steady_state', params, ys, 1); %note that g3 does not contain symmetric elements [~, g1, g2, g3] = feval([fname,'.dynamic'], ys(I), exo_steady_state', params, ys, 1); %note that g3 does not contain symmetric elements
g3 = unfold_g3(g3, yy0ex0_nbr); %add symmetric elements to g3 g3 = identification.unfold_g3(g3, yy0ex0_nbr); %add symmetric elements to g3
else else
T = NaN(sum(dynamic_tmp_nbr(1:5))); T = NaN(sum(dynamic_tmp_nbr(1:5)));
T = feval([fname, '.dynamic_g4_tt'], T, ys(I), exo_steady_state', params, ys, 1); T = feval([fname, '.dynamic_g4_tt'], T, ys(I), exo_steady_state', params, ys, 1);
...@@ -462,20 +472,16 @@ elseif (analytic_derivation_mode == 0 || analytic_derivation_mode == 1) ...@@ -462,20 +472,16 @@ elseif (analytic_derivation_mode == 0 || analytic_derivation_mode == 1)
g2 = feval([fname, '.dynamic_g2'], T, ys(I), exo_steady_state', params, ys, 1, false); %g2 is [endo_nbr by yy0ex0_nbr^2] second derivative (wrt all dynamic variables) of dynamic model equations, i.e. d(df/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order g2 = feval([fname, '.dynamic_g2'], T, ys(I), exo_steady_state', params, ys, 1, false); %g2 is [endo_nbr by yy0ex0_nbr^2] second derivative (wrt all dynamic variables) of dynamic model equations, i.e. d(df/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order
g3 = feval([fname, '.dynamic_g3'], T, ys(I), exo_steady_state', params, ys, 1, false); %note that g3 does not contain symmetric elements g3 = feval([fname, '.dynamic_g3'], T, ys(I), exo_steady_state', params, ys, 1, false); %note that g3 does not contain symmetric elements
g4 = feval([fname, '.dynamic_g4'], T, ys(I), exo_steady_state', params, ys, 1, false); %note that g4 does not contain symmetric elements g4 = feval([fname, '.dynamic_g4'], T, ys(I), exo_steady_state', params, ys, 1, false); %note that g4 does not contain symmetric elements
g3 = unfold_g3(g3, yy0ex0_nbr); %add symmetric elements to g3, %g3 is [endo_nbr by yy0ex0_nbr^3] third-derivative (wrt all dynamic variables) of dynamic model equations, i.e. (d(df/dyy0ex0)/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order g3 = identification.unfold_g3(g3, yy0ex0_nbr); %add symmetric elements to g3, %g3 is [endo_nbr by yy0ex0_nbr^3] third-derivative (wrt all dynamic variables) of dynamic model equations, i.e. (d(df/dyy0ex0)/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order
g4 = unfold_g4(g4, yy0ex0_nbr); %add symmetric elements to g4, %g4 is [endo_nbr by yy0ex0_nbr^4] fourth-derivative (wrt all dynamic variables) of dynamic model equations, i.e. ((d(df/dyy0ex0)/dyy0ex0)/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order g4 = identification.unfold_g4(g4, yy0ex0_nbr); %add symmetric elements to g4, %g4 is [endo_nbr by yy0ex0_nbr^4] fourth-derivative (wrt all dynamic variables) of dynamic model equations, i.e. ((d(df/dyy0ex0)/dyy0ex0)/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order
end end
%g1 is [endo_nbr by yy0ex0_nbr first derivative (wrt all dynamic variables) of dynamic model equations, i.e. df/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order %g1 is [endo_nbr by yy0ex0_nbr first derivative (wrt all dynamic variables) of dynamic model equations, i.e. df/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order
%g2 is [endo_nbr by yy0ex0_nbr^2] second derivative (wrt all dynamic variables) of dynamic model equations, i.e. d(df/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order %g2 is [endo_nbr by yy0ex0_nbr^2] second derivative (wrt all dynamic variables) of dynamic model equations, i.e. d(df/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order
%g3 is [endo_nbr by yy0ex0_nbr^3] third-derivative (wrt all dynamic variables) of dynamic model equations, i.e. (d(df/dyy0ex0)/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order %g3 is [endo_nbr by yy0ex0_nbr^3] third-derivative (wrt all dynamic variables) of dynamic model equations, i.e. (d(df/dyy0ex0)/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order
try
[~, g1p_static, rpp_static] = feval([fname,'.static_params_derivs'], ys, exo_steady_state', params); [~, g1p_static, rpp_static] = feval([fname,'.static_params_derivs'], ys, exo_steady_state', params);
%g1p_static is [endo_nbr by endo_nbr by param_nbr] first derivative (wrt all model parameters) of first-derivative (wrt all endogenous variables) of static model equations f, i.e. (df/dys)/dparams, in declaration order %g1p_static is [endo_nbr by endo_nbr by param_nbr] first derivative (wrt all model parameters) of first-derivative (wrt all endogenous variables) of static model equations f, i.e. (df/dys)/dparams, in declaration order
%rpp_static is [#second_order_residual_terms by 4] and contains nonzero values and corresponding indices of second derivatives (wrt all model parameters) of static model equations f, i.e. d(df/dparams)/dparams, in declaration order, where %rpp_static is [#second_order_residual_terms by 4] and contains nonzero values and corresponding indices of second derivatives (wrt all model parameters) of static model equations f, i.e. d(df/dparams)/dparams, in declaration order, where
% column 1 contains equation number; column 2 contains first parameter; column 3 contains second parameter; column 4 contains value of derivative % column 1 contains equation number; column 2 contains first parameter; column 3 contains second parameter; column 4 contains value of derivative
catch
error('For analytical parameter derivatives ''static_params_derivs.m'' file is needed, this can be created by putting identification(order=%d) into your mod file.',order)
end
rpp_static = get_all_resid_2nd_derivs(rpp_static, endo_nbr, param_nbr); %make full matrix out of nonzero values and corresponding indices rpp_static = get_all_resid_2nd_derivs(rpp_static, endo_nbr, param_nbr); %make full matrix out of nonzero values and corresponding indices
%rpp_static is [endo_nbr by param_nbr by param_nbr] second derivatives (wrt all model parameters) of static model equations, i.e. d(df/dparams)/dparams, in declaration order %rpp_static is [endo_nbr by param_nbr by param_nbr] second derivatives (wrt all model parameters) of static model equations, i.e. d(df/dparams)/dparams, in declaration order
if isempty(find(g2_static)) if isempty(find(g2_static))
...@@ -501,7 +507,7 @@ elseif (analytic_derivation_mode == 0 || analytic_derivation_mode == 1) ...@@ -501,7 +507,7 @@ elseif (analytic_derivation_mode == 0 || analytic_derivation_mode == 1)
end end
%handling of steady state for nonstationary variables %handling of steady state for nonstationary variables
if any(any(isnan(dys))) if any(any(isnan(dys)))
[U,T] = schur(g1_static); [U,T] = schur(full(g1_static));
e1 = abs(ordeig(T)) < qz_criterium-1; e1 = abs(ordeig(T)) < qz_criterium-1;
k = sum(e1); % Number of non stationary variables. k = sum(e1); % Number of non stationary variables.
% Number of stationary variables: n = length(e1)-k % Number of stationary variables: n = length(e1)-k
...@@ -519,58 +525,42 @@ elseif (analytic_derivation_mode == 0 || analytic_derivation_mode == 1) ...@@ -519,58 +525,42 @@ elseif (analytic_derivation_mode == 0 || analytic_derivation_mode == 1)
end end
if d2flag if d2flag
try
if order < 3 if order < 3
[~, g1p, ~, g1pp, g2p] = feval([fname,'.dynamic_params_derivs'], ys(I), exo_steady_state', params, ys, 1, dys, d2ys); [~, g1p, ~, g1pp, g2p] = feval([fname,'.dynamic_params_derivs'], ys(I), exo_steady_state', params, ys, 1, dys, d2ys);
else else
[~, g1p, ~, g1pp, g2p, g3p] = feval([fname,'.dynamic_params_derivs'], ys(I), exo_steady_state', params, ys, 1, dys, d2ys); [~, g1p, ~, g1pp, g2p, g3p] = feval([fname,'.dynamic_params_derivs'], ys(I), exo_steady_state', params, ys, 1, dys, d2ys);
end end
catch
error('For analytical parameter derivatives ''dynamic_params_derivs.m'' file is needed, this can be created by putting identification(order=%d) into your mod file.',order)
end
%g1pp are nonzero values and corresponding indices of second-derivatives (wrt all model parameters) of first-derivative (wrt all dynamic variables) of dynamic model equations, i.e. d(d(df/dyy0ex0)/dparam)/dparam, rows are in declaration order, first column in declaration order %g1pp are nonzero values and corresponding indices of second-derivatives (wrt all model parameters) of first-derivative (wrt all dynamic variables) of dynamic model equations, i.e. d(d(df/dyy0ex0)/dparam)/dparam, rows are in declaration order, first column in declaration order
d2Yss = d2ys(order_var,indpmodel,indpmodel); %[endo_nbr by mod_param_nbr by mod_param_nbr], put into DR order and focus only on selected model parameters d2Yss = d2ys(order_var,indpmodel,indpmodel); %[endo_nbr by mod_param_nbr by mod_param_nbr], put into DR order and focus only on selected model parameters
else else
if order == 1 if order == 1
try
[~, g1p] = feval([fname,'.dynamic_params_derivs'], ys(I), exo_steady_state', params, ys, 1, dys, d2ys); [~, g1p] = feval([fname,'.dynamic_params_derivs'], ys(I), exo_steady_state', params, ys, 1, dys, d2ys);
%g1p is [endo_nbr by yy0ex0_nbr by param_nbr] first-derivative (wrt all model parameters) of first-derivative (wrt all dynamic variables) of dynamic model equations, i.e. d(df/dyy0ex0)/dparam, rows are in declaration order, column in lead_lag_incidence order %g1p is [endo_nbr by yy0ex0_nbr by param_nbr] first-derivative (wrt all model parameters) of first-derivative (wrt all dynamic variables) of dynamic model equations, i.e. d(df/dyy0ex0)/dparam, rows are in declaration order, column in lead_lag_incidence order
catch
error('For analytical parameter derivatives ''dynamic_params_derivs.m'' file is needed, this can be created by putting identification(order=%d) into your mod file.',order)
end
[~, g1, g2 ] = feval([fname,'.dynamic'], ys(I), exo_steady_state', params, ys, 1); [~, g1, g2 ] = feval([fname,'.dynamic'], ys(I), exo_steady_state', params, ys, 1);
%g1 is [endo_nbr by yy0ex0_nbr first derivative (wrt all dynamic variables) of dynamic model equations, i.e. df/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order %g1 is [endo_nbr by yy0ex0_nbr first derivative (wrt all dynamic variables) of dynamic model equations, i.e. df/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order
%g2 is [endo_nbr by yy0ex0_nbr^2] second derivatives (wrt all dynamic variables) of dynamic model equations, i.e. d(df/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order %g2 is [endo_nbr by yy0ex0_nbr^2] second derivatives (wrt all dynamic variables) of dynamic model equations, i.e. d(df/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order
elseif order == 2 elseif order == 2
try
[~, g1p, ~, ~, g2p] = feval([fname,'.dynamic_params_derivs'], ys(I), exo_steady_state', params, ys, 1, dys, d2ys); [~, g1p, ~, ~, g2p] = feval([fname,'.dynamic_params_derivs'], ys(I), exo_steady_state', params, ys, 1, dys, d2ys);
%g1p is [endo_nbr by yy0ex0_nbr by param_nbr] first-derivative (wrt all model parameters) of first-derivative (wrt all dynamic variables) of dynamic model equations, i.e. d(df/dyy0ex0)/dparam, rows are in declaration order, column in lead_lag_incidence order %g1p is [endo_nbr by yy0ex0_nbr by param_nbr] first-derivative (wrt all model parameters) of first-derivative (wrt all dynamic variables) of dynamic model equations, i.e. d(df/dyy0ex0)/dparam, rows are in declaration order, column in lead_lag_incidence order
%g2p are nonzero values and corresponding indices of first-derivative (wrt all model parameters) of second-derivatives (wrt all dynamic variables) of dynamic model equations, i.e. d(d(df/dyy0ex0)/dyy0ex0)/dparam, rows are in declaration order, first and second column in declaration order %g2p are nonzero values and corresponding indices of first-derivative (wrt all model parameters) of second-derivatives (wrt all dynamic variables) of dynamic model equations, i.e. d(d(df/dyy0ex0)/dyy0ex0)/dparam, rows are in declaration order, first and second column in declaration order
catch
error('For analytical parameter derivatives ''dynamic_params_derivs.m'' file is needed, this can be created by putting identification(order=%d) into your mod file.',order)
end
[~, g1, g2, g3] = feval([fname,'.dynamic'], ys(I), exo_steady_state', params, ys, 1); %note that g3 does not contain symmetric elements [~, g1, g2, g3] = feval([fname,'.dynamic'], ys(I), exo_steady_state', params, ys, 1); %note that g3 does not contain symmetric elements
g3 = unfold_g3(g3, yy0ex0_nbr); %add symmetric elements to g3 g3 = identification.unfold_g3(g3, yy0ex0_nbr); %add symmetric elements to g3
%g1 is [endo_nbr by yy0ex0_nbr first derivative (wrt all dynamic variables) of dynamic model equations, i.e. df/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order %g1 is [endo_nbr by yy0ex0_nbr first derivative (wrt all dynamic variables) of dynamic model equations, i.e. df/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order
%g2 is [endo_nbr by yy0ex0_nbr^2] second derivative (wrt all dynamic variables) of dynamic model equations, i.e. d(df/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order %g2 is [endo_nbr by yy0ex0_nbr^2] second derivative (wrt all dynamic variables) of dynamic model equations, i.e. d(df/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order
%g3 is [endo_nbr by yy0ex0_nbr^3] third-derivative (wrt all dynamic variables) of dynamic model equations, i.e. (d(df/dyy0ex0)/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order %g3 is [endo_nbr by yy0ex0_nbr^3] third-derivative (wrt all dynamic variables) of dynamic model equations, i.e. (d(df/dyy0ex0)/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order
elseif order == 3 elseif order == 3
try
[~, g1p, ~, ~, g2p, g3p] = feval([fname,'.dynamic_params_derivs'], ys(I), exo_steady_state', params, ys, 1, dys, d2ys); [~, g1p, ~, ~, g2p, g3p] = feval([fname,'.dynamic_params_derivs'], ys(I), exo_steady_state', params, ys, 1, dys, d2ys);
%g1p is [endo_nbr by yy0ex0_nbr by param_nbr] first-derivative (wrt all model parameters) of first-derivative (wrt all dynamic variables) of dynamic model equations, i.e. d(df/dyy0ex0)/dparam, rows are in declaration order, column in lead_lag_incidence order %g1p is [endo_nbr by yy0ex0_nbr by param_nbr] first-derivative (wrt all model parameters) of first-derivative (wrt all dynamic variables) of dynamic model equations, i.e. d(df/dyy0ex0)/dparam, rows are in declaration order, column in lead_lag_incidence order
%g2p are nonzero values and corresponding indices of first-derivative (wrt all model parameters) of second-derivatives (wrt all dynamic variables) of dynamic model equations, i.e. d(d(df/dyy0ex0)/dyy0ex0)/dparam, rows are in declaration order, first and second column in declaration order %g2p are nonzero values and corresponding indices of first-derivative (wrt all model parameters) of second-derivatives (wrt all dynamic variables) of dynamic model equations, i.e. d(d(df/dyy0ex0)/dyy0ex0)/dparam, rows are in declaration order, first and second column in declaration order
%g3p are nonzero values and corresponding indices of first-derivative (wrt all model parameters) of third-derivatives (wrt all dynamic variables) of dynamic model equations, i.e. d(d(d(df/dyy0ex0)/dyy0ex0)/dyy0ex0)/dparam, rows are in declaration order, first, second and third column in declaration order %g3p are nonzero values and corresponding indices of first-derivative (wrt all model parameters) of third-derivatives (wrt all dynamic variables) of dynamic model equations, i.e. d(d(d(df/dyy0ex0)/dyy0ex0)/dyy0ex0)/dparam, rows are in declaration order, first, second and third column in declaration order
catch
error('For analytical parameter derivatives ''dynamic_params_derivs.m'' file is needed, this can be created by putting identification(order=%d) into your mod file.',order)
end
T = NaN(sum(dynamic_tmp_nbr(1:5))); T = NaN(sum(dynamic_tmp_nbr(1:5)));
T = feval([fname, '.dynamic_g4_tt'], T, ys(I), exo_steady_state', params, ys, 1); T = feval([fname, '.dynamic_g4_tt'], T, ys(I), exo_steady_state', params, ys, 1);
g1 = feval([fname, '.dynamic_g1'], T, ys(I), exo_steady_state', params, ys, 1, false); %g1 is [endo_nbr by yy0ex0_nbr first derivative (wrt all dynamic variables) of dynamic model equations, i.e. df/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order g1 = feval([fname, '.dynamic_g1'], T, ys(I), exo_steady_state', params, ys, 1, false); %g1 is [endo_nbr by yy0ex0_nbr first derivative (wrt all dynamic variables) of dynamic model equations, i.e. df/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order
g2 = feval([fname, '.dynamic_g2'], T, ys(I), exo_steady_state', params, ys, 1, false); %g2 is [endo_nbr by yy0ex0_nbr^2] second derivative (wrt all dynamic variables) of dynamic model equations, i.e. d(df/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order g2 = feval([fname, '.dynamic_g2'], T, ys(I), exo_steady_state', params, ys, 1, false); %g2 is [endo_nbr by yy0ex0_nbr^2] second derivative (wrt all dynamic variables) of dynamic model equations, i.e. d(df/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order
g3 = feval([fname, '.dynamic_g3'], T, ys(I), exo_steady_state', params, ys, 1, false); %note that g3 does not contain symmetric elements g3 = feval([fname, '.dynamic_g3'], T, ys(I), exo_steady_state', params, ys, 1, false); %note that g3 does not contain symmetric elements
g4 = feval([fname, '.dynamic_g4'], T, ys(I), exo_steady_state', params, ys, 1, false); %note that g4 does not contain symmetric elements g4 = feval([fname, '.dynamic_g4'], T, ys(I), exo_steady_state', params, ys, 1, false); %note that g4 does not contain symmetric elements
g3 = unfold_g3(g3, yy0ex0_nbr); %add symmetric elements to g3, %g3 is [endo_nbr by yy0ex0_nbr^3] third-derivative (wrt all dynamic variables) of dynamic model equations, i.e. (d(df/dyy0ex0)/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order g3 = identification.unfold_g3(g3, yy0ex0_nbr); %add symmetric elements to g3, %g3 is [endo_nbr by yy0ex0_nbr^3] third-derivative (wrt all dynamic variables) of dynamic model equations, i.e. (d(df/dyy0ex0)/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order
g4 = unfold_g4(g4, yy0ex0_nbr); %add symmetric elements to g4, %g4 is [endo_nbr by yy0ex0_nbr^4] fourth-derivative (wrt all dynamic variables) of dynamic model equations, i.e. ((d(df/dyy0ex0)/dyy0ex0)/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order g4 = identification.unfold_g4(g4, yy0ex0_nbr); %add symmetric elements to g4, %g4 is [endo_nbr by yy0ex0_nbr^4] fourth-derivative (wrt all dynamic variables) of dynamic model equations, i.e. ((d(df/dyy0ex0)/dyy0ex0)/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order
end end
end end
% Parameter Jacobian of steady state in different orderings, note dys is in declaration order % Parameter Jacobian of steady state in different orderings, note dys is in declaration order
...@@ -592,7 +582,7 @@ elseif (analytic_derivation_mode == 0 || analytic_derivation_mode == 1) ...@@ -592,7 +582,7 @@ elseif (analytic_derivation_mode == 0 || analytic_derivation_mode == 1)
[II, JJ] = ind2sub([yy0ex0_nbr yy0ex0_nbr], find(g2(j,:))); %g2 is [endo_nbr by yy0ex0_nbr^2] [II, JJ] = ind2sub([yy0ex0_nbr yy0ex0_nbr], find(g2(j,:))); %g2 is [endo_nbr by yy0ex0_nbr^2]
for i = 1:yy0ex0_nbr for i = 1:yy0ex0_nbr
is = find(II==i); is = find(II==i);
is = is(find(JJ(is)<=yy0_nbr)); %focus only on oo.dr.ys(I) derivatives as exogenous variables are 0 in a stochastic context is = is(find(JJ(is)<=yy0_nbr)); %focus only on dr.ys(I) derivatives as exogenous variables are 0 in a stochastic context
if ~isempty(is) if ~isempty(is)
tmp_g2 = full(g2(j,find(g2(j,:)))); tmp_g2 = full(g2(j,find(g2(j,:))));
dg1(j,i,:) = tmp_g2(is)*dyy0(JJ(is),:); %put into tensor notation dg1(j,i,:) = tmp_g2(is)*dyy0(JJ(is),:); %put into tensor notation
...@@ -795,7 +785,7 @@ if analytic_derivation_mode == 1 ...@@ -795,7 +785,7 @@ if analytic_derivation_mode == 1
dghu = [zeros(endo_nbr*exo_nbr, stderrparam_nbr+corrparam_nbr) dghu]; dghu = [zeros(endo_nbr*exo_nbr, stderrparam_nbr+corrparam_nbr) dghu];
% Compute dOm = dvec(ghu*Sigma_e*ghu') from expressions 34 in Iskrev (2010) Appendix A % Compute dOm = dvec(ghu*Sigma_e*ghu') from expressions 34 in Iskrev (2010) Appendix A
dOm = kron(I_endo,ghu*Sigma_e)*(commutation(endo_nbr, exo_nbr)*dghu)... dOm = kron(I_endo,ghu*Sigma_e)*(pruned_SS.commutation(endo_nbr, exo_nbr)*dghu)...
+ kron(ghu,ghu)*reshape(dSigma_e, exo_nbr^2, totparam_nbr) + kron(ghu*Sigma_e,I_endo)*dghu; + kron(ghu,ghu)*reshape(dSigma_e, exo_nbr^2, totparam_nbr) + kron(ghu*Sigma_e,I_endo)*dghu;
% Put into tensor notation % Put into tensor notation
...@@ -1540,7 +1530,7 @@ function g22 = get_2nd_deriv_mat(gpp,i,j,npar) ...@@ -1540,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);
...@@ -1553,7 +1543,7 @@ return ...@@ -1553,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
...@@ -1580,7 +1570,7 @@ return ...@@ -1580,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/get_perturbation_params_derivs_numerical_objective.m → matlab/+identification/get_perturbation_params_derivs_numerical_objective.m
View file @ 4a32acc2
function [out,info] = get_perturbation_params_derivs_numerical_objective(params, outputflag, estim_params, M, oo, options) function [out,info] = get_perturbation_params_derivs_numerical_objective(params, outputflag, estim_params, M_, options_, dr, endo_steady_state, exo_steady_state, exo_det_steady_state)
%function [out,info] = get_perturbation_params_derivs_numerical_objective(params, outputflag, estim_params, M, oo, options) %function [out,info] = get_perturbation_params_derivs_numerical_objective(params, outputflag, estim_params, M_, options_, dr, steady_state, exo_steady_state, exo_det_steady_state)
% ------------------------------------------------------------------------- % -------------------------------------------------------------------------
% Objective function to compute numerically the Jacobians used for get_perturbation_params_derivs % Objective function to compute numerically the Jacobians used for get_perturbation_params_derivs
% ========================================================================= % =========================================================================
...@@ -8,9 +8,12 @@ function [out,info] = get_perturbation_params_derivs_numerical_objective(params, ...@@ -8,9 +8,12 @@ function [out,info] = get_perturbation_params_derivs_numerical_objective(params,
% stderr parameters come first, corr parameters second, model parameters third % stderr parameters come first, corr parameters second, model parameters third
% outputflag: [string] flag which objective to compute (see below) % outputflag: [string] flag which objective to compute (see below)
% estim_params: [structure] storing the estimation information % estim_params: [structure] storing the estimation information
% M: [structure] storing the model information % M_: [structure] storing the model information
% oo: [structure] storing the solution results % options_: [structure] storing the options
% options: [structure] storing the options % dr [structure] Reduced form model.
% endo_steady_state [vector] steady state value for endogenous variables
% exo_steady_state [vector] steady state value for exogenous variables
% exo_det_steady_state [vector] steady state value for exogenous deterministic variables
% ------------------------------------------------------------------------- % -------------------------------------------------------------------------
% %
% OUTPUT % OUTPUT
...@@ -27,7 +30,7 @@ function [out,info] = get_perturbation_params_derivs_numerical_objective(params, ...@@ -27,7 +30,7 @@ function [out,info] = get_perturbation_params_derivs_numerical_objective(params,
% * get_perturbation_params_derivs.m (previously getH.m) % * get_perturbation_params_derivs.m (previously getH.m)
% ------------------------------------------------------------------------- % -------------------------------------------------------------------------
% This function calls % This function calls
% * [M.fname,'.dynamic'] % * [M_.fname,'.dynamic']
% * resol % * resol
% * dyn_vech % * dyn_vech
% ========================================================================= % =========================================================================
...@@ -50,58 +53,58 @@ function [out,info] = get_perturbation_params_derivs_numerical_objective(params, ...@@ -50,58 +53,58 @@ function [out,info] = get_perturbation_params_derivs_numerical_objective(params,
% ========================================================================= % =========================================================================
%% Update stderr, corr and model parameters and compute perturbation approximation and steady state with updated parameters %% Update stderr, corr and model parameters and compute perturbation approximation and steady state with updated parameters
M = set_all_parameters(params,estim_params,M); M_ = set_all_parameters(params,estim_params,M_);
[oo.dr,info,M.params] = compute_decision_rules(M,options,oo.dr, oo.steady_state, oo.exo_steady_state, oo.exo_det_steady_state); [dr,info,M_.params] = compute_decision_rules(M_,options_,dr,endo_steady_state,exo_steady_state,exo_det_steady_state);
Sigma_e = M.Sigma_e; Sigma_e = M_.Sigma_e;
if info(1) > 0 if info(1) > 0
% there are errors in the solution algorithm % there are errors in the solution algorithm
out = []; out = [];
return return
else else
ys = oo.dr.ys; %steady state of model variables in declaration order ys = dr.ys; %steady state of model variables in declaration order
ghx = oo.dr.ghx; ghu = oo.dr.ghu; ghx = dr.ghx; ghu = dr.ghu;
if options.order > 1 if options_.order > 1
ghxx = oo.dr.ghxx; ghxu = oo.dr.ghxu; ghuu = oo.dr.ghuu; ghs2 = oo.dr.ghs2; ghxx = dr.ghxx; ghxu = dr.ghxu; ghuu = dr.ghuu; ghs2 = dr.ghs2;
end end
if options.order > 2 if options_.order > 2
ghxxx = oo.dr.ghxxx; ghxxu = oo.dr.ghxxu; ghxuu = oo.dr.ghxuu; ghxss = oo.dr.ghxss; ghuuu = oo.dr.ghuuu; ghuss = oo.dr.ghuss; ghxxx = dr.ghxxx; ghxxu = dr.ghxxu; ghxuu = dr.ghxuu; ghxss = dr.ghxss; ghuuu = dr.ghuuu; ghuss = dr.ghuss;
end end
end end
Yss = ys(oo.dr.order_var); %steady state of model variables in DR order Yss = ys(dr.order_var); %steady state of model variables in DR order
%% out = [vec(Sigma_e);vec(ghx);vec(ghu);vec(ghxx);vec(ghxu);vec(ghuu);vec(ghs2);vec(ghxxx);vec(ghxxu);vec(ghxuu);vec(ghuuu);vec(ghxss);vec(ghuss)] %% out = [vec(Sigma_e);vec(ghx);vec(ghu);vec(ghxx);vec(ghxu);vec(ghuu);vec(ghs2);vec(ghxxx);vec(ghxxu);vec(ghxuu);vec(ghuuu);vec(ghxss);vec(ghuss)]
if strcmp(outputflag,'perturbation_solution') if strcmp(outputflag,'perturbation_solution')
out = [Sigma_e(:); ghx(:); ghu(:)]; out = [Sigma_e(:); ghx(:); ghu(:)];
if options.order > 1 if options_.order > 1
out = [out; ghxx(:); ghxu(:); ghuu(:); ghs2(:);]; out = [out; ghxx(:); ghxu(:); ghuu(:); ghs2(:);];
end end
if options.order > 2 if options_.order > 2
out = [out; ghxxx(:); ghxxu(:); ghxuu(:); ghuuu(:); ghxss(:); ghuss(:)]; out = [out; ghxxx(:); ghxxu(:); ghxuu(:); ghuuu(:); ghxss(:); ghuss(:)];
end end
end end
%% out = [Yss; vec(g1); vec(g2); vec(g3)]; of all endogenous variables, in DR order %% out = [Yss; vec(g1); vec(g2); vec(g3)]; of all endogenous variables, in DR order
if strcmp(outputflag,'dynamic_model') if strcmp(outputflag,'dynamic_model')
[I,~] = find(M.lead_lag_incidence'); %I is used to evaluate dynamic model files [I,~] = find(M_.lead_lag_incidence'); %I is used to evaluate dynamic model files
if options.order == 1 if options_.order == 1
[~, g1] = feval([M.fname,'.dynamic'], ys(I), oo.exo_steady_state', M.params, ys, 1); [~, g1] = feval([M_.fname,'.dynamic'], ys(I), exo_steady_state', M_.params, ys, 1);
out = [Yss; g1(:)]; out = [Yss; g1(:)];
elseif options.order == 2 elseif options_.order == 2
[~, g1, g2] = feval([M.fname,'.dynamic'], ys(I), oo.exo_steady_state', M.params, ys, 1); [~, g1, g2] = feval([M_.fname,'.dynamic'], ys(I), exo_steady_state', M_.params, ys, 1);
out = [Yss; g1(:); g2(:)]; out = [Yss; g1(:); g2(:)];
elseif options.order == 3 elseif options_.order == 3
[~, g1, g2, g3] = feval([M.fname,'.dynamic'], ys(I), oo.exo_steady_state', M.params, ys, 1); [~, g1, g2, g3] = feval([M_.fname,'.dynamic'], ys(I), exo_steady_state', M_.params, ys, 1);
g3 = unfold_g3(g3, length(ys(I))+M.exo_nbr); g3 = identification.unfold_g3(g3, length(ys(I))+M_.exo_nbr);
out = [Yss; g1(:); g2(:); g3(:)]; out = [Yss; g1(:); g2(:); g3(:)];
end end
end end
%% out = [Yss; vec(KalmanA); dyn_vech(KalmanB*Sigma_e*KalmanB')]; in DR order, where A and B are Kalman transition matrices %% out = [Yss; vec(KalmanA); dyn_vech(KalmanB*Sigma_e*KalmanB')]; in DR order, where A and B are Kalman transition matrices
if strcmp(outputflag,'Kalman_Transition') if strcmp(outputflag,'Kalman_Transition')
if options.order == 1 if options_.order == 1
KalmanA = zeros(M.endo_nbr,M.endo_nbr); KalmanA = zeros(M_.endo_nbr,M_.endo_nbr);
KalmanA(:,M.nstatic+(1:M.nspred)) = ghx; KalmanA(:,M_.nstatic+(1:M_.nspred)) = ghx;
Om = ghu*Sigma_e*transpose(ghu); Om = ghu*Sigma_e*transpose(ghu);
out = [Yss; KalmanA(:); dyn_vech(Om)]; out = [Yss; KalmanA(:); dyn_vech(Om)];
else else
......