Skip to content
GitLab
Commit a40cfedd authored by Sébastien Villemot's avatar Sébastien Villemot
Browse files

Merge branch 'identification_higher_order' into 'master' 

Identification Toolbox: improvements, bugfixes, license, preparation for getting rid of globals

Closes #1694 and #1689

See merge request Dynare/dynare!1697
parents d3e90a8d 3d51ee9e
Expand all Hide whitespace changes
Inline Side-by-side
license.txt
View file @ a40cfedd
......@@ -100,8 +100,32 @@ Copyright: 1997 Tom Minka <minka@microsoft.com>
License: GPL-3+
Files: matlab/allVL1.m
Copyright: 2009 Bruno Luong
License: BSD-2-clause
Copyright: 2007-2010 Bruno Luong <brunoluong@yahoo.com>
2019 Dynare Team
License: GPL-3+
Comment: The original author gave authorization to change
the license from BSD-2-clause to GPL-3+ and redistribute
it under GPL-3+ with Dynare.
Files: matlab/uperm.m
Copyright: 2014 Bruno Luong <brunoluong@yahoo.com>
2019 Dynare Team
License: GPL-3+
Comment: The original author gave authorization to change
the license from BSD-2-clause to GPL-3+ and redistribute
it under GPL-3+ with Dynare.
Files: matlab/prodmom.m matlab/bivmom.m
Copyright: 2008-2015 Raymond Kan <kan@chass.utoronto.ca>
2019-2020 Dynare Team
License: GPL-3+
Comment: The author gave authorization to redistribute
these functions under GPL-3+ with Dynare and would
appreciate acknowledgement of the source by citation
of the following paper:
Kan, R.: "From moments of sum to moments of product."
Journal of Multivariate Analysis, 2008, vol. 99, issue 3,
pages 542-554.
Files: matlab/gsa/Morris_Measure_Groups.m
matlab/gsa/Sampling_Function_2.m
......
matlab/allVL1.m
View file @ a40cfedd
......@@ -29,7 +29,7 @@ function v = allVL1(n, L1, L1ops, MaxNbSol)
% This function can be used to generate all orders of all
% multivariable polynomials of degree p in R^n:
% Order = allVL1(n, p)
% Author: Bruno Luong
% Author: Bruno Luong (brunoluong@yahoo.com)
% Original, 30/nov/2007
% Version 1.1, 30/apr/2008: Add H1 line as suggested by John D'Errico
% 1.2, 17/may/2009: Possibility to get the number of permutations
......@@ -38,31 +38,25 @@ function v = allVL1(n, L1, L1ops, MaxNbSol)
% 1.4, 18/Dec/2010: + non-recursive engine
% Retrieved from https://www.mathworks.com/matlabcentral/fileexchange/17818-all-permutations-of-integers-with-sum-criteria
% Copyright (c) 2009 Bruno Luong
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:
%
% * Redistributions of source code must retain the above copyright
% notice, this list of conditions and the following disclaimer.
% * Redistributions in binary form must reproduce the above copyright
% notice, this list of conditions and the following disclaimer in
% the documentation and/or other materials provided with the distribution
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.
% =========================================================================
% Copyright (C) 2007-2010 Bruno Luong <brunoluong@yahoo.com>
% Copyright (C) 2020 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 <http://www.gnu.org/licenses/>.
% =========================================================================
global MaxCounter;
if nargin<3 || isempty(L1ops)
......
matlab/bivmom.m
View file @ a40cfedd
......@@ -10,6 +10,29 @@
% of the equation.
% Usage: bivmom(p,rho)
%
% Retrieved from http://www-2.rotman.utoronto.ca/~kan/papers/prodmom.zip
% This function is part of replication codes of the following paper:
% Kan, R.: "From moments of sum to moments of product." Journal of
% Multivariate Analysis, 2008, vol. 99, issue 3, pages 542-554.
% =========================================================================
% Copyright (C) 2008-2015 Raymond Kan <kan@chass.utoronto.ca>
% Copyright (C) 2019-2020 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 <http://www.gnu.org/licenses/>.
% =========================================================================
function [y,dy] = bivmom(p,rho)
s1 = p(1);
s2 = p(2);
......
matlab/disp_identification.m
View file @ a40cfedd
......@@ -71,7 +71,7 @@ if options_ident.normalize_jacobians == 1
else
fprintf(' Normalize Jacobians: No\n');
end
fprintf(' Tolerance level for rank computations: %.0d\n',options_ident.tol_rank);
fprintf(' Tolerance level for rank computations: %s\n',num2str(options_ident.tol_rank));
fprintf(' Tolerance level for selecting nonzero columns: %.0d\n',options_ident.tol_deriv);
fprintf(' Tolerance level for selecting nonzero singular values: %.0d\n',options_ident.tol_sv);
......
matlab/dynare_identification.m
View file @ a40cfedd
This diff is collapsed.
matlab/fjaco.m
View file @ a40cfedd
......@@ -10,7 +10,7 @@ function fjac = fjaco(f,x,varargin)
% OUTPUT
% fjac : finite difference Jacobian
%
% Copyright (C) 2010-2017,2019 Dynare Team
% Copyright (C) 2010-2017,2019-2020 Dynare Team
%
% This file is part of Dynare.
%
......
matlab/get_identification_jacobians.m
View file @ a40cfedd
This diff is collapsed.
matlab/get_minimal_state_representation.m
View file @ a40cfedd
function [CheckCO,minnx,minA,minB,minC,minD,dminA,dminB,dminC,dminD] = get_minimal_state_representation(A,B,C,D,dA,dB,dC,dD)
% [CheckCO,minnx,minA,minB,minC,minD,dminA,dminB,dminC,dminD] = get_minimal_state_representation(A,B,C,D,dA,dB,dC,dD)
% Derives and checks the minimal state representation of the ABCD
% representation of a state space model
function [CheckCO,minns,minSYS] = get_minimal_state_representation(SYS, derivs_flag)
% Derives and checks the minimal state representation
% Let x = A*x(-1) + B*u and y = C*x(-1) + D*u be a linear state space
% system, then this function computes the following representation
% xmin = minA*xmin(-1) + minB*u and and y=minC*xmin(-1) + minD*u
%
% -------------------------------------------------------------------------
% INPUTS
% A: [endo_nbr by endo_nbr] Transition matrix from Kalman filter
% for all endogenous declared variables, in DR order
% B: [endo_nbr by exo_nbr] Transition matrix from Kalman filter
% mapping shocks today to endogenous variables today, in DR order
% C: [obs_nbr by endo_nbr] Measurement matrix from Kalman filter
% linking control/observable variables to states, in DR order
% D: [obs_nbr by exo_nbr] Measurement matrix from Kalman filter
% mapping shocks today to controls/observables today, in DR order
% dA: [endo_nbr by endo_nbr by totparam_nbr] in DR order
% SYS [structure]
% with the following necessary fields:
% A: [nspred by nspred] in DR order
% Transition matrix for all state variables
% B: [nspred by exo_nbr] in DR order
% Transition matrix mapping shocks today to states today
% C: [varobs_nbr by nspred] in DR order
% Measurement matrix linking control/observable variables to states
% D: [varobs_nbr by exo_nbr] in DR order
% Measurement matrix mapping shocks today to controls/observables today
% and optional fields:
% dA: [nspred by nspred by totparam_nbr] in DR order
% Jacobian (wrt to all parameters) of transition matrix A
% dB: [endo_nbr by exo_nbr by totparam_nbr] in DR order
% dB: [nspred by exo_nbr by totparam_nbr] in DR order
% Jacobian (wrt to all parameters) of transition matrix B
% dC: [obs_nbr by endo_nbr by totparam_nbr] in DR order
% dC: [varobs_nbr by nspred by totparam_nbr] in DR order
% Jacobian (wrt to all parameters) of measurement matrix C
% dD: [obs_nbr by exo_nbr by totparam_nbr] in DR order
% dD: [varobs_nbr by exo_nbr by totparam_nbr] in DR order
% Jacobian (wrt to all parameters) of measurement matrix D
% derivs_flag [scalar]
% (optional) indicator whether to output parameter derivatives
% -------------------------------------------------------------------------
% OUTPUTS
% CheckCO: [scalar] indicator, equals to 1 if minimal state representation is found
% minnx: [scalar] length of minimal state vector
% minA: [minnx by minnx] Transition matrix A for evolution of minimal state vector
% minB: [minnx by exo_nbr] Transition matrix B for evolution of minimal state vector
% minC: [obs_nbr by minnx] Measurement matrix C for evolution of controls, depending on minimal state vector only
% minD: [obs_nbr by minnx] Measurement matrix D for evolution of controls, depending on minimal state vector only
% dminA: [minnx by minnx by totparam_nbr] in DR order
% CheckCO: [scalar]
% equals to 1 if minimal state representation is found
% minns: [scalar]
% length of minimal state vector
% SYS [structure]
% with the following fields:
% minA: [minns by minns] in DR-order
% transition matrix A for evolution of minimal state vector
% minB: [minns by exo_nbr] in DR-order
% transition matrix B for evolution of minimal state vector
% minC: [varobs_nbr by minns] in DR-order
% measurement matrix C for evolution of controls, depending on minimal state vector only
% minD: [varobs_nbr by minns] in DR-order
% measurement matrix D for evolution of controls, depending on minimal state vector only
% dminA: [minns by minns by totparam_nbr] in DR order
% Jacobian (wrt to all parameters) of transition matrix minA
% dminB: [minnx by exo_nbr by totparam_nbr] in DR order
% dminB: [minns by exo_nbr by totparam_nbr] in DR order
% Jacobian (wrt to all parameters) of transition matrix minB
% dminC: [obs_nbr by minnx by totparam_nbr] in DR order
% dminC: [varobs_nbr by minns by totparam_nbr] in DR order
% Jacobian (wrt to all parameters) of measurement matrix minC
% dminD: [obs_nbr by exo_nbr by totparam_nbr] in DR order
% dminD: [varobs_nbr by u_nbr by totparam_nbr] in DR order
% Jacobian (wrt to all parameters) of measurement matrix minD
% -------------------------------------------------------------------------
% This function is called by
......@@ -42,8 +57,9 @@ function [CheckCO,minnx,minA,minB,minC,minD,dminA,dminB,dminC,dminD] = get_minim
% -------------------------------------------------------------------------
% This function calls
% * check_minimality (embedded)
% * minrealold (embedded)
% =========================================================================
% Copyright (C) 2019 Dynare Team
% Copyright (C) 2019-2020 Dynare Team
%
% This file is part of Dynare.
%
......@@ -60,102 +76,123 @@ function [CheckCO,minnx,minA,minB,minC,minD,dminA,dminB,dminC,dminD] = get_minim
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
% =========================================================================
nx = size(A,2);
ny = size(C,1);
nu = size(B,2);
if nargin == 1
derivs_flag = 0;
end
realsmall = 1e-7;
[nspred,exo_nbr] = size(SYS.B);
varobs_nbr = size(SYS.C,1);
% Check controllability and observability conditions for full state vector
CheckCO = check_minimality(A,B,C);
if CheckCO == 1 % If model is already minimal
minnx = nx;
minA = A;
minB = B;
minC = C;
minD = D;
if nargout > 6
dminA = dA;
dminB = dB;
dminC = dC;
dminD = dD;
end
CheckCO = check_minimality(SYS.A,SYS.B,SYS.C);
if CheckCO == 1 % If model is already minimal, we are finished
minns = nspred;
minSYS = SYS;
else
%Model is not minimal
realsmall = 1e-7;
% create indices for unnecessary states
endogstateindex = find(abs(sum(A,1))<realsmall);
exogstateindex = find(abs(sum(A,1))>realsmall);
minnx = length(exogstateindex);
% remove unnecessary states from solution matrices
A_1 = A(endogstateindex,exogstateindex);
A_2 = A(exogstateindex,exogstateindex);
B_2 = B(exogstateindex,:);
C_1 = C(:,endogstateindex);
C_2 = C(:,exogstateindex);
D = D;
ind_A1 = endogstateindex;
ind_A2 = exogstateindex;
% minimal representation
minA = A_2;
minB = B_2;
minC = C_2;
minD = D;
% Check controllability and observability conditions
CheckCO = check_minimality(minA,minB,minC);
if CheckCO ~=1
j=1;
while (CheckCO==0 && j<minnx)
combis = nchoosek(1:minnx,j);
i=1;
while i <= size(combis,1)
ind_A2 = exogstateindex;
ind_A1 = [endogstateindex ind_A2(combis(j,:))];
ind_A2(combis(j,:)) = [];
% remove unnecessary states from solution matrices
A_1 = A(ind_A1,ind_A2);
A_2 = A(ind_A2,ind_A2);
B_2 = B(ind_A2,:);
C_1 = C(:,ind_A1);
C_2 = C(:,ind_A2);
D = D;
% minimal representation
minA = A_2;
minB = B_2;
minC = C_2;
minD = D;
% Check controllability and observability conditions
CheckCO = check_minimality(minA,minB,minC);
if CheckCO == 1
minnx = length(ind_A2);
break;
%Model is not minimal
try
minreal_flag = 1;
% In future we will use SLICOT TB01PD.f mex file [to do @wmutschl], currently use workaround
[mSYS,U] = minrealold(SYS,realsmall);
minns = size(mSYS.A,1);
CheckCO = check_minimality(mSYS.A,mSYS.B,mSYS.C);
if CheckCO
minSYS.A = mSYS.A;
minSYS.B = mSYS.B;
minSYS.C = mSYS.C;
minSYS.D = mSYS.D;
if derivs_flag
totparam_nbr = size(SYS.dA,3);
minSYS.dA = zeros(minns,minns,totparam_nbr);
minSYS.dB = zeros(minns,exo_nbr,totparam_nbr);
minSYS.dC = zeros(varobs_nbr,minns,totparam_nbr);
% Note that orthogonal matrix U is such that (U*dA*U',U*dB,dC*U') is a Kalman decomposition of (dA,dB,dC) %
for jp=1:totparam_nbr
dA_tmp = U*SYS.dA(:,:,jp)*U';
dB_tmp = U*SYS.dB(:,:,jp);
dC_tmp = SYS.dC(:,:,jp)*U';
minSYS.dA(:,:,jp) = dA_tmp(1:minns,1:minns);
minSYS.dB(:,:,jp) = dB_tmp(1:minns,:);
minSYS.dC(:,:,jp) = dC_tmp(:,1:minns);
end
i=i+1;
minSYS.dD = SYS.dD;
end
j=j+1;
end
catch
minreal_flag = 0; % if something went wrong use below procedure
end
if nargout > 6
dminA = dA(ind_A2,ind_A2,:);
dminB = dB(ind_A2,:,:);
dminC = dC(:,ind_A2,:);
dminD = dD;
if minreal_flag == 0
fprintf('Use naive brute-force approach to find minimal state space system.\n These computations may be inaccurate/wrong as ''minreal'' is not available, please raise an issue on GitLab or the forum\n')
% create indices for unnecessary states
exogstateindex = find(abs(sum(SYS.A,1))>realsmall);
minns = length(exogstateindex);
% remove unnecessary states from solution matrices
A_2 = SYS.A(exogstateindex,exogstateindex);
B_2 = SYS.B(exogstateindex,:);
C_2 = SYS.C(:,exogstateindex);
D = SYS.D;
ind_A2 = exogstateindex;
% minimal representation
minSYS.A = A_2;
minSYS.B = B_2;
minSYS.C = C_2;
minSYS.D = D;
% Check controllability and observability conditions
CheckCO = check_minimality(minSYS.A,minSYS.B,minSYS.C);
if CheckCO ~=1
% do brute-force search
j=1;
while (CheckCO==0 && j<minns)
combis = nchoosek(1:minns,j);
i=1;
while i <= size(combis,1)
ind_A2 = exogstateindex;
ind_A2(combis(j,:)) = [];
% remove unnecessary states from solution matrices
A_2 = SYS.A(ind_A2,ind_A2);
B_2 = SYS.B(ind_A2,:);
C_2 = SYS.C(:,ind_A2);
D = SYS.D;
% minimal representation
minSYS.A = A_2;
minSYS.B = B_2;
minSYS.C = C_2;
minSYS.D = D;
% Check controllability and observability conditions
CheckCO = check_minimality(minSYS.A,minSYS.B,minSYS.C);
if CheckCO == 1
minns = length(ind_A2);
break;
end
i=i+1;
end
j=j+1;
end
end
if derivs_flag
minSYS.dA = SYS.dA(ind_A2,ind_A2,:);
minSYS.dB = SYS.dB(ind_A2,:,:);
minSYS.dC = SYS.dC(:,ind_A2,:);
minSYS.dD = SYS.dD;
end
end
end
function CheckCO = check_minimality(A,B,C)
function CheckCO = check_minimality(a,b,c)
%% This function computes the controllability and the observability matrices of the ABCD system and checks if the system is minimal
%
% Inputs: Solution matrices A,B,C of ABCD representation of a DSGE model
% Outputs: CheckCO: equals 1, if both rank conditions for observability and controllability are fullfilled
n = size(A,1);
CC = B; % Initialize controllability matrix
OO = C; % Initialize observability matrix
n = size(a,1);
CC = b; % Initialize controllability matrix
OO = c; % Initialize observability matrix
if n >= 2
for indn = 1:1:n-1
CC = [CC, (A^indn)*B]; % Set up controllability matrix
OO = [OO; C*(A^indn)]; % Set up observability matrix
CC = [CC, (a^indn)*b]; % Set up controllability matrix
OO = [OO; c*(a^indn)]; % Set up observability matrix
end
end
CheckC = (rank(full(CC))==n); % Check rank of controllability matrix
......@@ -163,4 +200,87 @@ CheckO = (rank(full(OO))==n); % Check rank of observability matrix
CheckCO = CheckC&CheckO; % equals 1 if minimal state
end % check_minimality end
function [mSYS,U] = minrealold(SYS,tol)
% This is a temporary replacement for minreal, will be replaced by a mex file from SLICOT TB01PD.f soon
a = SYS.A;
b = SYS.B;
c = SYS.C;
[ns,nu] = size(b);
[am,bm,cm,U,k] = ControllabilityStaircaseRosenbrock(a,b,c,tol);
kk = sum(k);
nu = ns - kk;
nn = nu;
am = am(nu+1:ns,nu+1:ns);
bm = bm(nu+1:ns,:);
cm = cm(:,nu+1:ns);
ns = ns - nu;
if ns
[am,bm,cm,t,k] = ObservabilityStaircaseRosenbrock(am,bm,cm,tol);
kk = sum(k);
nu = ns - kk;
nn = nn + nu;
am = am(nu+1:ns,nu+1:ns);
bm = bm(nu+1:ns,:);
cm = cm(:,nu+1:ns);
end
mSYS.A = am;
mSYS.B = bm;
mSYS.C = cm;
mSYS.D = SYS.D;
end
function [abar,bbar,cbar,t,k] = ObservabilityStaircaseRosenbrock(a,b,c,tol)
%Observability staircase form
[aa,bb,cc,t,k] = ControllabilityStaircaseRosenbrock(a',c',b',tol);
abar = aa'; bbar = cc'; cbar = bb';
end
function [abar,bbar,cbar,t,k] = ControllabilityStaircaseRosenbrock(a, b, c, tol)
% Controllability staircase algorithm of Rosenbrock, 1968
[ra,ca] = size(a);
[rb,cb] = size(b);
ptjn1 = eye(ra);
ajn1 = a;
bjn1 = b;
rojn1 = cb;
deltajn1 = 0;
sigmajn1 = ra;
k = zeros(1,ra);
if nargin == 3
tol = ra*norm(a,1)*eps;
end
for jj = 1 : ra
[uj,sj,vj] = svd(bjn1);
[rsj,csj] = size(sj);
p = flip(eye(rsj),2);
p = permute(p,[2 1 3:ndims(eye(rsj))]);
uj = uj*p;
bb = uj'*bjn1;
roj = rank(bb,tol);
[rbb,cbb] = size(bb);
sigmaj = rbb - roj;
sigmajn1 = sigmaj;
k(jj) = roj;
if roj == 0, break, end
if sigmaj == 0, break, end
abxy = uj' * ajn1 * uj;
aj = abxy(1:sigmaj,1:sigmaj);
bj = abxy(1:sigmaj,sigmaj+1:sigmaj+roj);
ajn1 = aj;
bjn1 = bj;
[ruj,cuj] = size(uj);
ptj = ptjn1 * ...
[uj zeros(ruj,deltajn1); ...
zeros(deltajn1,cuj) eye(deltajn1)];
ptjn1 = ptj;
deltaj = deltajn1 + roj;
deltajn1 = deltaj;
end
t = ptjn1';
abar = t * a * t';
bbar = t * b;
cbar = c * t';
end
end % Main function end
matlab/get_perturbation_params_derivs.m
View file @ a40cfedd
This diff is collapsed.
matlab/get_perturbation_params_derivs_numerical_objective.m
View file @ a40cfedd
......@@ -13,22 +13,25 @@ function [out,info] = get_perturbation_params_derivs_numerical_objective(params,
% options: [structure] storing the options
% -------------------------------------------------------------------------
%
% OUTPUT (dependent on outputflag and order of approximation):
% - 'perturbation_solution': 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)]
% - 'dynamic_model': out = [Yss; vec(g1); vec(g2); vec(g3)]
% - 'Kalman_Transition': out = [Yss; vec(KalmanA); dyn_vech(KalmanB*Sigma_e*KalmanB')];
% all in DR-order
% OUTPUT
% out (dependent on outputflag and order of approximation):
% - 'perturbation_solution': out = out1 = [vec(Sigma_e);vec(ghx);vec(ghu)]; (order==1)
% out = out2 = [out1;vec(ghxx);vec(ghxu);vec(ghuu);vec(ghs2)]; (order==2)
% out = out3 = [out1;out2;vec(ghxxx);vec(ghxxu);vec(ghxuu);vec(ghuuu);vec(ghxss);vec(ghuss)]; (order==3)
% - 'dynamic_model': out = [Yss; vec(g1); vec(g2); vec(g3)]
% - 'Kalman_Transition': out = [Yss; vec(KalmanA); dyn_vech(KalmanB*Sigma_e*KalmanB')];
% all in DR-order
% info [integer] output from resol
% -------------------------------------------------------------------------
% This function is called by
% * get_solution_params_deriv.m (previously getH.m)
% * get_identification_jacobians.m (previously getJJ.m)
% * get_perturbation_params_derivs.m (previously getH.m)
% -------------------------------------------------------------------------
% This function calls
% * [M.fname,'.dynamic']
% * resol
% * dyn_vech
% =========================================================================
% Copyright (C) 2019 Dynare Team
% Copyright (C) 2019-2020 Dynare Team
%
% This file is part of Dynare.
%
......@@ -48,8 +51,8 @@ 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
M = set_all_parameters(params,estim_params,M);
Sigma_e = M.Sigma_e;
[~,info,M,options,oo] = resol(0,M,options,oo);
Sigma_e = M.Sigma_e;
if info(1) > 0
% there are errors in the solution algorithm
......@@ -60,11 +63,9 @@ else
ghx = oo.dr.ghx; ghu = oo.dr.ghu;
if options.order > 1
ghxx = oo.dr.ghxx; ghxu = oo.dr.ghxu; ghuu = oo.dr.ghuu; ghs2 = oo.dr.ghs2;
%ghxs = oo.dr.ghxs; ghus = oo.dr.ghus; %these are zero due to certainty equivalence and Gaussian shocks
end
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;
%ghxxs = oo.dr.ghxxs; ghxus = oo.dr.ghxus; ghuus = oo.dr.ghuus; ghsss = oo.dr.ghsss; %these are zero due to certainty equivalence and Gaussian shocks
end
end
Yss = ys(oo.dr.order_var); %steady state of model variables in DR order
......
matlab/identification_analysis.m
View file @ a40cfedd
......@@ -28,7 +28,7 @@ function [ide_moments, ide_spectrum, ide_minimal, ide_hess, ide_reducedform, ide
% 1: prior exists. Identification is checked for stderr, corr and model parameters as declared in estimated_params block
% 0: prior does not exist. Identification is checked for all stderr and model parameters, but no corr parameters
% * init [integer]
% flag for initialization of persistent vars. This is needed if one may wants to make more calls to identification in the same dynare mod file
% flag for initialization of persistent vars. This is needed if one may want to make more calls to identification in the same mod file