Skip to content
Snippets Groups Projects

Compare revisions

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

Source

Select target project
No results found
Select Git revision
Loading items

Target

Select target project
  • giovanma/dynare
  • giorgiomas/dynare
  • Vermandel/dynare
  • Dynare/dynare
  • normann/dynare
  • MichelJuillard/dynare
  • wmutschl/dynare
  • FerhatMihoubi/dynare
  • sebastien/dynare
  • lnsongxf/dynare
  • rattoma/dynare
  • CIMERS/dynare
  • FredericKarame/dynare
  • SumuduK/dynare
  • MinjeJeon/dynare
  • camilomrch/dynare
  • DoraK/dynare
  • avtishin/dynare
  • selma/dynare
  • claudio_olguin/dynare
  • jeffjiang07/dynare
  • EthanSystem/dynare
  • stepan-a/dynare
  • wjgatt/dynare
  • JohannesPfeifer/dynare
  • gboehl/dynare
  • chskcau/dynare-doc-fixes
27 results
Select Git revision
Loading items
Show changes
Showing
with 261 additions and 295 deletions
File moved
function P=lyapunov_solver(T,R,Q,DynareOptions) function P=lyapunov_solver(T,R,Q,options_)
% function P=lyapunov_solver(T,R,Q,DynareOptions) % function P=lyapunov_solver(T,R,Q,options_)
% Solves the Lyapunov equation P-T*P*T' = R*Q*R' arising in a state-space % Solves the Lyapunov equation P-T*P*T' = R*Q*R' arising in a state-space
% system, where P is the variance of the states % system, where P is the variance of the states
% %
...@@ -7,7 +7,7 @@ function P=lyapunov_solver(T,R,Q,DynareOptions) ...@@ -7,7 +7,7 @@ function P=lyapunov_solver(T,R,Q,DynareOptions)
% T [double] n*n matrix. % T [double] n*n matrix.
% R [double] n*m matrix. % R [double] n*m matrix.
% Q [double] m*m matrix. % Q [double] m*m matrix.
% DynareOptions [structure] Dynare options % options_ [structure] Dynare options
% %
% Outputs % Outputs
% P [double] n*n matrix. % P [double] n*n matrix.
...@@ -15,11 +15,11 @@ function P=lyapunov_solver(T,R,Q,DynareOptions) ...@@ -15,11 +15,11 @@ function P=lyapunov_solver(T,R,Q,DynareOptions)
% Algorithms % Algorithms
% Default, if none of the other algorithms is selected: % Default, if none of the other algorithms is selected:
% Reordered Schur decomposition (Bartels-Stewart algorithm) % Reordered Schur decomposition (Bartels-Stewart algorithm)
% DynareOptions.lyapunov_fp == true % options_.lyapunov_fp == true
% iteration-based fixed point algorithm % iteration-based fixed point algorithm
% DynareOptions.lyapunov_db == true % options_.lyapunov_db == true
% doubling algorithm % doubling algorithm
% DynareOptions.lyapunov_srs == true % options_.lyapunov_srs == true
% Square-root solver for discrete-time Lyapunov equations (requires Matlab System Control toolbox % Square-root solver for discrete-time Lyapunov equations (requires Matlab System Control toolbox
% or Octave control package) % or Octave control package)
...@@ -40,14 +40,14 @@ function P=lyapunov_solver(T,R,Q,DynareOptions) ...@@ -40,14 +40,14 @@ function P=lyapunov_solver(T,R,Q,DynareOptions)
% 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/>.
if DynareOptions.lyapunov_fp if options_.lyapunov_fp
P = lyapunov_symm(T,R*Q*R',DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, 3, DynareOptions.debug); P = lyapunov_symm(T,R*Q*R',options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold, 3, options_.debug);
elseif DynareOptions.lyapunov_db elseif options_.lyapunov_db
[P, errorflag] = disclyap_fast(T,R*Q*R',DynareOptions.lyapunov_doubling_tol); [P, errorflag] = disclyap_fast(T,R*Q*R',options_.lyapunov_doubling_tol);
if errorflag %use Schur-based method if errorflag %use Schur-based method
P = lyapunov_symm(T,R*Q*R',DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, [], DynareOptions.debug); P = lyapunov_symm(T,R*Q*R',options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold, [], options_.debug);
end end
elseif DynareOptions.lyapunov_srs elseif options_.lyapunov_srs
% works only with Matlab System Control toolbox or Octave control package, % works only with Matlab System Control toolbox or Octave control package,
if isoctave if isoctave
if ~user_has_octave_forge_package('control') if ~user_has_octave_forge_package('control')
...@@ -62,7 +62,7 @@ elseif DynareOptions.lyapunov_srs ...@@ -62,7 +62,7 @@ elseif DynareOptions.lyapunov_srs
R_P = dlyapchol(T,chol_Q); R_P = dlyapchol(T,chol_Q);
P = R_P' * R_P; P = R_P' * R_P;
else else
P = lyapunov_symm(T,R*Q*R',DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, [], DynareOptions.debug); P = lyapunov_symm(T,R*Q*R',options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold, [], options_.debug);
end end
return % --*-- Unit tests --*-- return % --*-- Unit tests --*--
...@@ -92,7 +92,7 @@ tmp2=randn(m_large,m_large); ...@@ -92,7 +92,7 @@ tmp2=randn(m_large,m_large);
Q_large=tmp2*tmp2'; Q_large=tmp2*tmp2';
R_large=randn(n_large,m_large); R_large=randn(n_large,m_large);
% DynareOptions.lyapunov_fp == 1 % options_.lyapunov_fp == 1
options_.lyapunov_fp = true; options_.lyapunov_fp = true;
try try
Pstar1_small = lyapunov_solver(T_small,R_small,Q_small,options_); Pstar1_small = lyapunov_solver(T_small,R_small,Q_small,options_);
...@@ -102,7 +102,7 @@ catch ...@@ -102,7 +102,7 @@ catch
t(1) = 0; t(1) = 0;
end end
% Dynareoptions.lyapunov_db == 1 % options_.lyapunov_db == 1
options_.lyapunov_fp = false; options_.lyapunov_fp = false;
options_.lyapunov_db = true; options_.lyapunov_db = true;
try try
...@@ -113,7 +113,7 @@ catch ...@@ -113,7 +113,7 @@ catch
t(2) = 0; t(2) = 0;
end end
% Dynareoptions.lyapunov_srs == 1 % options_.lyapunov_srs == 1
if (isoctave && user_has_octave_forge_package('control')) || (~isoctave && user_has_matlab_license('control_toolbox')) if (isoctave && user_has_octave_forge_package('control')) || (~isoctave && user_has_matlab_license('control_toolbox'))
options_.lyapunov_db = false; options_.lyapunov_db = false;
options_.lyapunov_srs = true; options_.lyapunov_srs = true;
......
File moved
function [res,jac,domerr] = mcpath_function(func,z,jacflag,varargin)
domerr = 0;
if jacflag
[res,jac] = func(z,varargin{:});
else
res = func(z,varargin{:});
end
function [xparams, logpost, options_]=metropolis_draw(init,options_,estim_params_,M_)
% function [xparams, logpost, options_]=metropolis_draw(init,options_,estim_params_,M_)
% Builds draws from metropolis
%
% INPUTS:
% init: scalar equal to
% 1: first call to store the required information
% on files, lines, and chains to be read
% in persistent variables to make them available
% for future calls
% 0: load a parameter draw
% Additional Inputs required for initialization
% options_ [structure] Matlab's structure describing the options (initialized by dynare, see @ref{options_}).
% M_ [structure] Matlab's structure describing the Model (initialized by dynare, see @ref{M_}).
% estim_params_ [structure] Matlab's structure describing the estimated_parameters (initialized by dynare, see @ref{estim_params_}).
%
% OUTPUTS:
% xparams: if init==0: vector of estimated parameters
% if init==1: error flaog
% logpost: log of posterior density
% options_: [structure] Matlab's structure describing the options (initialized by dynare, see @ref{options_}).
%
% SPECIAL REQUIREMENTS
%
% Requires CutSample to be run before in order to set up mh_history-file
% Copyright © 2003-2023 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/>.
persistent mh_nblck NumberOfDraws BaseName FirstLine FirstMhFile MAX_nruns
xparams = 0;
logpost = 0;
if init
%get number of parameters
nvx = estim_params_.nvx;
nvn = estim_params_.nvn;
ncx = estim_params_.ncx;
ncn = estim_params_.ncn;
np = estim_params_.np ;
npar = nvx+nvn+ncx+ncn+np;
%get path of metropolis files
MetropolisFolder = CheckPath('metropolis',M_.dname);
FileName = M_.fname;
BaseName = [MetropolisFolder filesep FileName];
%load mh_history-file with info on what to load
record=load_last_mh_history_file(MetropolisFolder, FileName);
FirstMhFile = record.KeepedDraws.FirstMhFile;
FirstLine = record.KeepedDraws.FirstLine;
TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8); %number of parameters plus posterior plus ?
mh_nblck = options_.mh_nblck;
% set sub_draws option if empty
if isempty(options_.sub_draws)
options_.sub_draws = min(options_.posterior_max_subsample_draws, ceil(.25*NumberOfDraws));
else
if options_.sub_draws>NumberOfDraws*mh_nblck
skipline()
disp(['The value of option sub_draws (' num2str(options_.sub_draws) ') is greater than the number of available draws in the MCMC (' num2str(NumberOfDraws*mh_nblck) ')!'])
disp('You can either change the value of sub_draws, reduce the value of mh_drop, or run another mcmc (with the load_mh_file option).')
skipline()
xparams = 1; % xparams is interpreted as an error flag
end
end
return
else %not initialization, return one draw
%get random draw from random chain
ChainNumber = ceil(rand*mh_nblck);
DrawNumber = ceil(rand*NumberOfDraws);
if DrawNumber <= MAX_nruns-FirstLine+1 %draw in first file, needs to account for first line
MhFilNumber = FirstMhFile;
MhLine = FirstLine+DrawNumber-1;
else %draw in other file
DrawNumber = DrawNumber-(MAX_nruns-FirstLine+1);
MhFilNumber = FirstMhFile+ceil(DrawNumber/MAX_nruns);
MhLine = DrawNumber-(MhFilNumber-FirstMhFile-1)*MAX_nruns;
end
%load parameters and posterior
load( [ BaseName '_mh' int2str(MhFilNumber) '_blck' int2str(ChainNumber) '.mat' ],'x2','logpo2');
xparams = x2(MhLine,:);
logpost= logpo2(MhLine);
end
\ No newline at end of file
% [dynpp_derivs, dyn_derivs] = k_order_perturbation(dr,DynareModel,DynareOptions) % [dynpp_derivs, dyn_derivs] = k_order_perturbation(dr,M_,options_)
% computes a k-th order perturbation solution % computes a k-th order perturbation solution
% %
% INPUTS % INPUTS
% dr: struct describing the reduced form solution of the model. % dr: struct describing the reduced form solution of the model.
% DynareModel: struct jobs's parameters % M_: struct jobs's parameters
% DynareOptions: struct job's options % options_: struct job's options
% %
% OUTPUTS % OUTPUTS
% dynpp_derivs struct Derivatives of the decision rule in Dynare++ format. % dynpp_derivs struct Derivatives of the decision rule in Dynare++ format.
...@@ -25,7 +25,7 @@ ...@@ -25,7 +25,7 @@
% dynare/mex/sources/k_order_perturbation.cc and it uses code provided by % dynare/mex/sources/k_order_perturbation.cc and it uses code provided by
% dynare++ % dynare++
% Copyright © 2013-2021 Dynare Team % Copyright © 2013-2023 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
......
% W = k_order_welfare(dr, DynareModel, DynareOptions) % W = k_order_welfare(dr, M_, options_)
% computes a k-th order approximation of welfare % computes a k-th order approximation of welfare
% %
% INPUTS % INPUTS
% dr: struct describing the reduced form solution of the model. % dr: struct describing the reduced form solution of the model.
% DynareModel: struct jobs's parameters % M_: struct jobs's parameters
% DynareOptions: struct job's options % options_: struct job's options
% %
% OUTPUTS % OUTPUTS
% %
......
...@@ -12,7 +12,7 @@ function tf = contains(string, pattern, varargin) ...@@ -12,7 +12,7 @@ function tf = contains(string, pattern, varargin)
% OUTPUT % OUTPUT
% - tf [logical] % - tf [logical]
% %
% Copyright © 2019 Dynare Team % Copyright © 2019-2024 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
...@@ -62,7 +62,7 @@ end ...@@ -62,7 +62,7 @@ end
tf = false(size(string)); tf = false(size(string));
for ii = 1:numel(pattern) for ii = 1:numel(pattern)
idx = regexp(string, pattern{ii}); idx = regexp(string, regexptranslate('escape', pattern{ii}));
for jj = 1:numel(string) for jj = 1:numel(string)
tf(jj) = tf(jj) || ~isempty(idx{jj}); tf(jj) = tf(jj) || ~isempty(idx{jj});
end end
......
function [c, ia, ib] = intersect_stable(a, b)
% Crude implementation of intersect(…, 'stable'), which is missing before
% Octave 6.2 and buggy before 8.4 (#60347)
% Copyright (C) 2019-2023 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/>.
if nargin ~= 2
error('intersect_stable: needs exactly 2 input arguments');
end
if isnumeric (a) && isnumeric (b)
c = [];
elseif iscell (b)
c = {};
else
c = '';
end
ia = [];
ib = [];
if isempty (a) || isempty (b)
return
else
isrowvec = isrow (a) && isrow (b);
for i = 1:numel(a)
if iscellstr(c)
idx = strcmp(a(i), b);
else
idx = a(i) == b;
end
if any(idx) && ~ismember(a(i), c)
c = [c(:); a(i)];
if nargout > 1
ia = [ia, i];
ib = [ib, find(idx)];
end
end
end
%% Adjust output orientation for MATLAB compatibility
if isrowvec
c = c.';
end
end
end
%!test
%! a = [3 4 1 5];
%! b = [2 4 9 1 6];
%! [c,ia,ib]=intersect_stable(a,b);
%! assert(c, [4 1])
%! assert(ia, [2 3])
%! assert(ib, [2 4])
%! assert(a(ia), c)
%! assert(b(ib), c)
%!test
%! a = [3 4 1 5]';
%! b = [2 4 9 1 6]';
%! [c,ia,ib]=intersect_stable(a,b);
%! assert(c, [4 1]')
%! assert(ia, [2 3])
%! assert(ib, [2 4])
%! assert(a(ia), c)
%! assert(b(ib), c)
%!test
%! a = { 'defun', 'mapcar', 'let', 'eval-when'};
%! b = { 'setf', 'let', 'list', 'cdr', 'defun'};
%! [c,ia,ib]=intersect_stable(a,b);
%! assert(c, { 'defun', 'let' })
%! assert(ia, [1 3])
%! assert(ib, [5 2])
%! assert(a(ia), c)
%! assert(b(ib), c)
function [X, info] = cycle_reduction(A0, A1, A2, cvg_tol, ch) function [X, info] = cycle_reduction(A0, A1, A2, cvg_tol, max_it, ch)
%@info: %@info:
%! @deftypefn {Function File} {[@var{X}, @var{info}] =} cycle_reduction (@var{A0},@var{A1},@var{A2},@var{cvg_tol},@var{ch}) %! @deftypefn {Function File} {[@var{X}, @var{info}] =} cycle_reduction (@var{A0},@var{A1},@var{A2},@var{cvg_tol},@var{ch})
...@@ -60,14 +60,13 @@ function [X, info] = cycle_reduction(A0, A1, A2, cvg_tol, ch) ...@@ -60,14 +60,13 @@ function [X, info] = cycle_reduction(A0, A1, A2, cvg_tol, ch)
% 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/>.
max_it = 300;
it = 0; it = 0;
info = 0; info = 0;
X = []; X = [];
crit = Inf; crit = Inf;
A0_0 = A0; A0_0 = A0;
Ahat1 = A1; Ahat1 = A1;
if (nargin == 5 && ~isempty(ch) ) if (nargin == 6 && ~isempty(ch) )
A1_0 = A1; A1_0 = A1;
A2_0 = A2; A2_0 = A2;
end end
...@@ -94,7 +93,7 @@ while cont ...@@ -94,7 +93,7 @@ while cont
return return
elseif isnan(crit) elseif isnan(crit)
info(1) = 402; info(1) = 402;
info(2) = log(norm(A1,1)) info(2) = log(norm(A1,1));
return return
end end
it = it + 1; it = it + 1;
...@@ -102,12 +101,12 @@ end ...@@ -102,12 +101,12 @@ end
X = -Ahat1\A0_0; X = -Ahat1\A0_0;
if (nargin == 5 && ~isempty(ch) ) if (nargin == 6 && ~isempty(ch) )
%check the solution %check the solution
res = norm(A0_0 + A1_0 * X + A2_0 * X * X, 1); res = norm(A0_0 + A1_0 * X + A2_0 * X * X, 1);
if (res > cvg_tol) if (res > cvg_tol)
info(1) = 403 info(1) = 403;
info(2) = log(res) info(2) = log(res);
dprintf('The norm of the residual is %s whereas the tolerance criterion is %s', num2str(res), num2str(cvg_tol)); dprintf('The norm of the residual is %s whereas the tolerance criterion is %s', num2str(res), num2str(cvg_tol));
end end
end end
...@@ -128,7 +127,7 @@ C = diag(15*ones(n,1)); C = C - diag(5*ones(n-1,1),-1); C = C - diag(5*ones(n-1, ...@@ -128,7 +127,7 @@ C = diag(15*ones(n,1)); C = C - diag(5*ones(n-1,1),-1); C = C - diag(5*ones(n-1,
% Solve the equation with the cycle reduction algorithm % Solve the equation with the cycle reduction algorithm
try try
tic; X1 = cycle_reduction(C,B,A,1e-7); elapsedtime = toc; tic; X1 = cycle_reduction(C,B,A,1e-7,100); elapsedtime = toc;
disp(['Elapsed time for cycle reduction algorithm is: ' num2str(elapsedtime) ' (n=' int2str(n) ').']) disp(['Elapsed time for cycle reduction algorithm is: ' num2str(elapsedtime) ' (n=' int2str(n) ').'])
t(1) = 1; t(1) = 1;
catch catch
......
...@@ -39,10 +39,6 @@ if nargin > 5 || nargin < 2 || nargout > 7 || nargout == 0 ...@@ -39,10 +39,6 @@ if nargin > 5 || nargin < 2 || nargout > 7 || nargout == 0
error('MJDGGES: takes 2, 3 or 4 input arguments and between 1 and 7 output arguments.') error('MJDGGES: takes 2, 3 or 4 input arguments and between 1 and 7 output arguments.')
end end
if isoctave && octave_ver_less_than('7')
error('Octave version 7 or higher is required (Octave 6 lacks the ordqz function)')
end
[me, ne] = size(e); [me, ne] = size(e);
[md, nd] = size(d); [md, nd] = size(d);
if ~isreal(e) || ~isreal(d) || me ~= ne || md ~= nd || me ~= nd if ~isreal(e) || ~isreal(d) || me ~= ne || md ~= nd || me ~= nd
...@@ -64,7 +60,7 @@ try ...@@ -64,7 +60,7 @@ try
eigval = ordeig(ss, tt); eigval = ordeig(ss, tt);
select = abs(eigval) < qz_criterium; select = abs(eigval) < qz_criterium;
sdim = sum(select); sdim = sum(select);
[ss, tt, qq, zz] = ordqz(ss, tt, qq, zz, select); [ss, tt, ~, zz] = ordqz(ss, tt, qq, zz, select);
eigval = ordeig(ss, tt); eigval = ordeig(ss, tt);
catch catch
info = 1; % Not as precise as lapack's info! info = 1; % Not as precise as lapack's info!
......
...@@ -9,7 +9,7 @@ function g = best_1978(a ,b) ...@@ -9,7 +9,7 @@ function g = best_1978(a ,b)
% OUTPUTS % OUTPUTS
% - g [double] n*1 vector, gamma variates. % - g [double] n*1 vector, gamma variates.
% Copyright © 2006-2020 Dynare Team % Copyright © 2006-2023 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
...@@ -44,12 +44,7 @@ while mm ...@@ -44,12 +44,7 @@ while mm
X(index) = bb(index)+Y(index); % x X(index) = bb(index)+Y(index); % x
id1 = index(X(index)<0); % Reject. id1 = index(X(index)<0); % Reject.
id2 = setdiff(index, id1); id2 = setdiff(index, id1);
if numel(id2) ~= 0 || isoctave || ~matlab_ver_less_than('9.1')
% If id2=[], LHS of the .* has size [0,0], while RHS is [0,1].
% Since there is no automatic broadcast in MATLAB < R2016b, skip the
% statement in that case.
Z(id2) = 64.0*(W(id2).^3).*(rand(length(id2),1).^2); % d Z(id2) = 64.0*(W(id2).^3).*(rand(length(id2),1).^2); % d
end
id3 = id2(Z(id2)>1.0-2.0*Y(id2).*Y(id2)./X(id2)); % Reject. id3 = id2(Z(id2)>1.0-2.0*Y(id2).*Y(id2)./X(id2)); % Reject.
id4 = id3(log(Z(id3))>2.0*(bb(id3).*log(X(id3)./bb(id3))-Y(id3))); % Reject. id4 = id3(log(Z(id3))>2.0*(bb(id3).*log(X(id3)./bb(id3))-Y(id3))); % Reject.
index = [id1, id4]; index = [id1, id4];
......
...@@ -9,7 +9,7 @@ function g = knuth(a, b) ...@@ -9,7 +9,7 @@ function g = knuth(a, b)
% OUTPUTS % OUTPUTS
% - g [double] n*1 vector, gamma variates. % - g [double] n*1 vector, gamma variates.
% Copyright © 2006-2020 Dynare Team % Copyright © 2006-2023 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
...@@ -38,15 +38,7 @@ while mm ...@@ -38,15 +38,7 @@ while mm
X(index) = Y(index).*bb(index) + a(index) - 1; X(index) = Y(index).*bb(index) + a(index) - 1;
id1 = index(X(index)<=0); % Rejected draws. id1 = index(X(index)<=0); % Rejected draws.
id2 = setdiff(index, id1); id2 = setdiff(index, id1);
if numel(id2) == 0 && ~isoctave && matlab_ver_less_than('9.1')
% The LHS of the > comparison in the "else" branch has size = [0, 1]
% while the RHS has size = [0, 0].
% Automatic broadcasting was only introduced in MATLAB R2016b, so we
% must handle this case separately to avoid an error.
id3 = [];
else
id3 = id2(rand(length(id2), 1)>(1+Y(id2).*Y(id2)).*exp((a(id2)-1).*(log(X(id2))-log(a(id2)-1))-bb(id2).*Y(id2))); % Rejected draws. id3 = id2(rand(length(id2), 1)>(1+Y(id2).*Y(id2)).*exp((a(id2)-1).*(log(X(id2))-log(a(id2)-1))-bb(id2).*Y(id2))); % Rejected draws.
end
index = [id1, id3]; index = [id1, id3];
mm = length(index); mm = length(index);
end end
......
...@@ -132,7 +132,7 @@ return % --*-- Unit tests --*-- ...@@ -132,7 +132,7 @@ return % --*-- Unit tests --*--
if ~isoctave && ~user_has_matlab_license('statistics_toolbox') if ~isoctave && ~user_has_matlab_license('statistics_toolbox')
method = struct('small', 'Weibull-rejection', 'large', 'Knuth'); method = struct('small', 'Weibull-rejection', 'large', 'Knuth');
n = 1000000; n = 1000000;
m = 100; m = 1;
a = 0.1; a = 0.1;
b = 1.0; b = 1.0;
try try
...@@ -171,7 +171,7 @@ end ...@@ -171,7 +171,7 @@ end
if ~isoctave && ~user_has_matlab_license('statistics_toolbox') if ~isoctave && ~user_has_matlab_license('statistics_toolbox')
method = struct('small', 'Johnk', 'large', 'Knuth'); method = struct('small', 'Johnk', 'large', 'Knuth');
n = 1000000; n = 1000000;
m = 100; m = 1;
a = 0.1; a = 0.1;
b = 1.0; b = 1.0;
try try
...@@ -210,7 +210,7 @@ end ...@@ -210,7 +210,7 @@ end
if ~isoctave && ~user_has_matlab_license('statistics_toolbox') if ~isoctave && ~user_has_matlab_license('statistics_toolbox')
method = struct('small', 'Berman', 'large', 'Knuth'); method = struct('small', 'Berman', 'large', 'Knuth');
n = 1000000; n = 1000000;
m = 100; m = 3;
a = 0.1; a = 0.1;
b = 1.0; b = 1.0;
try try
...@@ -249,7 +249,7 @@ end ...@@ -249,7 +249,7 @@ end
if ~isoctave && ~user_has_matlab_license('statistics_toolbox') if ~isoctave && ~user_has_matlab_license('statistics_toolbox')
method = struct('small', 'Ahrens-Dieter', 'large', 'Knuth'); method = struct('small', 'Ahrens-Dieter', 'large', 'Knuth');
n = 1000000; n = 1000000;
m = 100; m = 1;
a = 0.1; a = 0.1;
b = 1.0; b = 1.0;
try try
...@@ -288,7 +288,7 @@ end ...@@ -288,7 +288,7 @@ end
if ~isoctave && ~user_has_matlab_license('statistics_toolbox') if ~isoctave && ~user_has_matlab_license('statistics_toolbox')
method = struct('small', 'Best', 'large', 'Knuth'); method = struct('small', 'Best', 'large', 'Knuth');
n = 1000000; n = 1000000;
m = 100; m = 1;
a = 0.1; a = 0.1;
b = 1.0; b = 1.0;
try try
...@@ -327,7 +327,7 @@ end ...@@ -327,7 +327,7 @@ end
if ~isoctave && ~user_has_matlab_license('statistics_toolbox') if ~isoctave && ~user_has_matlab_license('statistics_toolbox')
method = struct('small', 'Weibull-rejection', 'large', 'Knuth'); method = struct('small', 'Weibull-rejection', 'large', 'Knuth');
n = 1000000; n = 1000000;
m = 100; m = 3;
a = 1.5; a = 1.5;
b = 1.0; b = 1.0;
try try
...@@ -366,7 +366,7 @@ end ...@@ -366,7 +366,7 @@ end
if ~isoctave && ~user_has_matlab_license('statistics_toolbox') if ~isoctave && ~user_has_matlab_license('statistics_toolbox')
method = struct('small', 'Weibull-rejection', 'large', 'Cheng'); method = struct('small', 'Weibull-rejection', 'large', 'Cheng');
n = 1000000; n = 1000000;
m = 100; m = 1;
a = 1.5; a = 1.5;
b = 1.0; b = 1.0;
try try
...@@ -405,7 +405,7 @@ end ...@@ -405,7 +405,7 @@ end
if ~isoctave && ~user_has_matlab_license('statistics_toolbox') if ~isoctave && ~user_has_matlab_license('statistics_toolbox')
method = struct('small', 'Weibull-rejection', 'large', 'Best'); method = struct('small', 'Weibull-rejection', 'large', 'Best');
n = 1000000; n = 1000000;
m = 100; m = 20;
a = 1.5; a = 1.5;
b = 1.0; b = 1.0;
try try
......
function model_diagnostics(M,options,oo) function model_diagnostics(M_,options_,oo_)
% function model_diagnostics(M,options,oo) % function model_diagnostics(M_,options_,oo_)
% computes various diagnostics on the model % computes various diagnostics on the model
% INPUTS % INPUTS
% M [matlab structure] Definition of the model. % M_ [matlab structure] Definition of the model.
% options [matlab structure] Global options. % options_ [matlab structure] Global options.
% oo [matlab structure] Results % oo_ [matlab structure] Results
% %
% OUTPUTS % OUTPUTS
% none % none
...@@ -16,7 +16,7 @@ function model_diagnostics(M,options,oo) ...@@ -16,7 +16,7 @@ function model_diagnostics(M,options,oo)
% none. % none.
% %
% Copyright © 1996-2023 Dynare Team % Copyright © 1996-2024 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
...@@ -33,22 +33,22 @@ function model_diagnostics(M,options,oo) ...@@ -33,22 +33,22 @@ function model_diagnostics(M,options,oo)
% You should have received a copy of the GNU General Public License % You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>. % along with Dynare. If not, see <https://www.gnu.org/licenses/>.
endo_names = M.endo_names; endo_names = M_.endo_names;
lead_lag_incidence = M.lead_lag_incidence; lead_lag_incidence = M_.lead_lag_incidence;
maximum_endo_lag = M.maximum_endo_lag; maximum_endo_lag = M_.maximum_endo_lag;
if options.ramsey_policy if options_.ramsey_policy
%test whether specification matches %test whether specification matches
inst_nbr = size(options.instruments,1); inst_nbr = size(options_.instruments,1);
if inst_nbr~=0 if inst_nbr~=0
implied_inst_nbr = M.ramsey_orig_endo_nbr - M.ramsey_orig_eq_nbr; implied_inst_nbr = M_.ramsey_orig_endo_nbr - M_.ramsey_orig_eq_nbr;
if inst_nbr>implied_inst_nbr if inst_nbr>implied_inst_nbr
warning('You have specified more steady state instruments than there are omitted equations. While there are use cases for this setup, it is rather unusual. Check whether this is desired.') warning('You have specified more steady state instruments than there are omitted equations. While there are use cases for this setup, it is rather unusual. Check whether this is desired.')
elseif inst_nbr<implied_inst_nbr elseif inst_nbr<implied_inst_nbr
warning('You have specified fewer steady state instruments than there are omitted equations. While there are use cases for this setup, it is rather unusual. Check whether this is desired.') warning('You have specified fewer steady state instruments than there are omitted equations. While there are use cases for this setup, it is rather unusual. Check whether this is desired.')
end end
else else
if options.steadystate_flag if options_.steadystate_flag
warning('You have specified a steady state file, but not provided steady state instruments. In this case, you typically need to make sure to provide all steady state values, including the ones for the planner''s instrument(s).') warning('You have specified a steady state file, but not provided steady state instruments. In this case, you typically need to make sure to provide all steady state values, including the ones for the planner''s instrument(s).')
end end
end end
...@@ -57,21 +57,21 @@ end ...@@ -57,21 +57,21 @@ end
problem_dummy=0; problem_dummy=0;
%naming conflict in steady state file %naming conflict in steady state file
if options.steadystate_flag == 1 if options_.steadystate_flag == 1
if strmatch('ys',M.endo_names,'exact') if strmatch('ys',M_.endo_names,'exact')
disp(['MODEL_DIAGNOSTICS: using the name ys for an endogenous variable will typically conflict with the internal naming in user-defined steady state files.']) disp('MODEL_DIAGNOSTICS: using the name ys for an endogenous variable will typically conflict with the internal naming in user-defined steady state files.')
problem_dummy=1; problem_dummy=1;
end end
if strmatch('ys',M.param_names,'exact') if strmatch('ys',M_.param_names,'exact')
disp(['MODEL_DIAGNOSTICS: using the name ys for a parameter will typically conflict with the internal naming in user-defined steady state files.']) disp('MODEL_DIAGNOSTICS: using the name ys for a parameter will typically conflict with the internal naming in user-defined steady state files.')
problem_dummy=1; problem_dummy=1;
end end
if strmatch('M_',M.endo_names,'exact') if strmatch('M_',M_.endo_names,'exact')
disp(['MODEL_DIAGNOSTICS: using the name M_ for an endogenous variable will typically conflict with the internal naming in user-defined steady state files.']) disp('MODEL_DIAGNOSTICS: using the name M_ for an endogenous variable will typically conflict with the internal naming in user-defined steady state files.')
problem_dummy=1; problem_dummy=1;
end end
if strmatch('M_',M.param_names,'exact') if strmatch('M_',M_.param_names,'exact')
disp(['MODEL_DIAGNOSTICS: using the name M_ for a parameter will typically conflict with the internal naming in user-defined steady state files.']) disp('MODEL_DIAGNOSTICS: using the name M_ for a parameter will typically conflict with the internal naming in user-defined steady state files.')
problem_dummy=1; problem_dummy=1;
end end
end end
...@@ -94,34 +94,40 @@ end ...@@ -94,34 +94,40 @@ end
% %
info = 0; info = 0;
if M.exo_nbr == 0 if M_.exo_nbr == 0
oo.exo_steady_state = [] ; oo_.exo_steady_state = [] ;
end end
info=test_for_deep_parameters_calibration(M); info=test_for_deep_parameters_calibration(M_);
if info if info
problem_dummy=1; problem_dummy=1;
end end
% check if ys is steady state % check if ys is steady state
options.debug=true; %locally set debug option to true options_.debug=true; %locally set debug option to true
if options.logged_steady_state %if steady state was previously logged, undo this if options_.logged_steady_state %if steady state was previously logged, undo this
oo.dr.ys=exp(oo.dr.ys); oo_.dr.ys=exp(oo_.dr.ys);
oo.steady_state=exp(oo.steady_state); oo_.steady_state=exp(oo_.steady_state);
options.logged_steady_state=0; options_.logged_steady_state=0;
end end
[dr.ys,M.params,check1]=evaluate_steady_state(oo.steady_state,[oo.exo_steady_state; oo.exo_det_steady_state],M,options,options.steadystate.nocheck); [dr.ys,M_.params,check1]=evaluate_steady_state(oo_.steady_state,[oo_.exo_steady_state; oo_.exo_det_steady_state],M_,options_,~options_.steadystate.nocheck);
if isfield(M_,'occbin')
if any(oo_.exo_steady_state)
disp('MODEL_DIAGNOSTICS: OccBin was detected in conjunction with a non-zero steady state of the exogenous variables. That will usually create issues.')
problem_dummy=1;
end
end
% testing for problem % testing for problem
if check1(1) if check1(1)
problem_dummy=1; problem_dummy=1;
disp('MODEL_DIAGNOSTICS: The steady state cannot be computed') disp('MODEL_DIAGNOSTICS: The steady state cannot be computed')
if any(isnan(dr.ys)) if any(isnan(dr.ys))
disp(['MODEL_DIAGNOSTICS: Steady state contains NaNs']) disp('MODEL_DIAGNOSTICS: Steady state contains NaNs')
end end
if any(isinf(dr.ys)) if any(isinf(dr.ys))
disp(['MODEL_DIAGNOSTICS: Steady state contains Inf']) disp('MODEL_DIAGNOSTICS: Steady state contains Inf')
end end
return return
end end
...@@ -137,35 +143,36 @@ end ...@@ -137,35 +143,36 @@ end
% singular Jacobian of static model % singular Jacobian of static model
% %
singularity_problem = 0; singularity_problem = 0;
if ~options.block if ~options_.block
nb = 1; nb = 1;
else else
nb = length(M.block_structure_stat.block); nb = length(M_.block_structure_stat.block);
end end
exo = [oo.exo_steady_state; oo.exo_det_steady_state]; exo = [oo_.exo_steady_state; oo_.exo_det_steady_state];
for b=1:nb for b=1:nb
if options.bytecode if options_.bytecode
if nb == 1 if nb == 1
[res, jacob] = bytecode(dr.ys, exo, M.params, dr.ys, 1, exo, ... [~, jacob] = bytecode(M_, options_, dr.ys, exo, M_.params, dr.ys, 1, exo, ...
'evaluate', 'static'); 'evaluate', 'static');
else else
[res, jacob] = bytecode(dr.ys, exo, M.params, dr.ys, 1, exo, ... [~, jacob] = bytecode(M_, options_, dr.ys, exo, M_.params, dr.ys, 1, exo, ...
'evaluate', 'static', 'block_decomposed', ['block=' ... 'evaluate', 'static', 'block_decomposed', ['block=' ...
int2str(b)]); int2str(b)]);
end end
n_vars_jacob=size(jacob,2);
else else
if options.block if options_.block
T = NaN(M.block_structure_stat.tmp_nbr, 1); T = NaN(M_.block_structure_stat.tmp_nbr, 1);
fh_static = str2func(sprintf('%s.sparse.block.static_%d', M.fname, b)); fh_static = str2func(sprintf('%s.sparse.block.static_%d', M_.fname, b));
[~, ~,~, jacob] = fh_static(dr.ys, exo, M.params, M.block_structure_stat.block(b).g1_sparse_rowval, ... [~, ~,~, jacob] = fh_static(dr.ys, exo, M_.params, M_.block_structure_stat.block(b).g1_sparse_rowval, ...
M.block_structure_stat.block(b).g1_sparse_colval, ... M_.block_structure_stat.block(b).g1_sparse_colval, ...
M.block_structure_stat.block(b).g1_sparse_colptr, T); M_.block_structure_stat.block(b).g1_sparse_colptr, T);
n_vars_jacob=size(jacob,2); n_vars_jacob=size(jacob,2);
else else
[res, T_order, T] = feval([M.fname '.sparse.static_resid'], dr.ys, exo, M.params); [~, T_order, T] = feval([M_.fname '.sparse.static_resid'], dr.ys, exo, M_.params);
jacob = feval([M.fname '.sparse.static_g1'], dr.ys, exo, M.params, M.static_g1_sparse_rowval, M.static_g1_sparse_colval, M.static_g1_sparse_colptr, T_order, T); jacob = feval([M_.fname '.sparse.static_g1'], dr.ys, exo, M_.params, M_.static_g1_sparse_rowval, M_.static_g1_sparse_colval, M_.static_g1_sparse_colptr, T_order, T);
n_vars_jacob=M.endo_nbr; n_vars_jacob=M_.endo_nbr;
end end
jacob=full(jacob); jacob=full(jacob);
end end
...@@ -173,19 +180,19 @@ for b=1:nb ...@@ -173,19 +180,19 @@ for b=1:nb
problem_dummy=1; problem_dummy=1;
[infrow,infcol]=find(isinf(jacob) | isnan(jacob)); [infrow,infcol]=find(isinf(jacob) | isnan(jacob));
fprintf('\nMODEL_DIAGNOSTICS: The Jacobian of the static model contains Inf or NaN. The problem arises from: \n\n') fprintf('\nMODEL_DIAGNOSTICS: The Jacobian of the static model contains Inf or NaN. The problem arises from: \n\n')
display_problematic_vars_Jacobian(infrow,infcol,M,dr.ys,'static','MODEL_DIAGNOSTICS: ') display_problematic_vars_Jacobian(infrow,infcol,M_,dr.ys,'static','MODEL_DIAGNOSTICS: ')
end end
if any(any(~isreal(jacob))) if any(any(~isreal(jacob)))
problem_dummy=1; problem_dummy=1;
[imagrow,imagcol]=find(abs(imag(jacob))>1e-15); [imagrow,imagcol]=find(abs(imag(jacob))>1e-15);
fprintf('\nMODEL_DIAGNOSTICS: The Jacobian of the static model contains imaginary parts. The problem arises from: \n\n') fprintf('\nMODEL_DIAGNOSTICS: The Jacobian of the static model contains imaginary parts. The problem arises from: \n\n')
display_problematic_vars_Jacobian(imagrow,imagcol,M,dr.ys,'static','MODEL_DIAGNOSTICS: ') display_problematic_vars_Jacobian(imagrow,imagcol,M_,dr.ys,'static','MODEL_DIAGNOSTICS: ')
end end
try try
if (~isoctave && matlab_ver_less_than('9.12')) || isempty(options.jacobian_tolerance) if (~isoctave && matlab_ver_less_than('9.12')) || isempty(options_.jacobian_tolerance)
rank_jacob = rank(jacob); %can sometimes fail rank_jacob = rank(jacob); %can sometimes fail
else else
rank_jacob = rank(jacob,options.jacobian_tolerance); %can sometimes fail rank_jacob = rank(jacob,options_.jacobian_tolerance); %can sometimes fail
end end
catch catch
rank_jacob=size(jacob,1); rank_jacob=size(jacob,1);
...@@ -197,48 +204,48 @@ for b=1:nb ...@@ -197,48 +204,48 @@ for b=1:nb
'singular']) 'singular'])
disp(['MODEL_DIAGNOSTICS: there is ' num2str(n_vars_jacob-rank_jacob) ... disp(['MODEL_DIAGNOSTICS: there is ' num2str(n_vars_jacob-rank_jacob) ...
' collinear relationships between the variables and the equations']) ' collinear relationships between the variables and the equations'])
if (~isoctave && matlab_ver_less_than('9.12')) || isempty(options.jacobian_tolerance) if (~isoctave && matlab_ver_less_than('9.12')) || isempty(options_.jacobian_tolerance)
ncol = null(jacob); ncol = null(jacob);
else else
ncol = null(jacob,options.jacobian_tolerance); %can sometimes fail ncol = null(jacob,options_.jacobian_tolerance); %can sometimes fail
end end
n_rel = size(ncol,2); n_rel = size(ncol,2);
for i = 1:n_rel for i = 1:n_rel
if n_rel > 1 if n_rel > 1
disp(['Relation ' int2str(i)]) disp(['Relation ' int2str(i)])
end end
disp('Colinear variables:') disp('Collinear variables:')
for j=1:10 for j=1:10
k = find(abs(ncol(:,i)) > 10^-j); k = find(abs(ncol(:,i)) > 10^-j);
if max(abs(jacob(:,k)*ncol(k,i))) < 1e-6 if max(abs(jacob(:,k)*ncol(k,i))) < 1e-6
break break
end end
end end
if options.block && ~options.bytecode if options_.block && ~options_.bytecode
fprintf('%s\n',endo_names{M.block_structure_stat.block(b).variable(k)}) fprintf('%s\n',endo_names{M_.block_structure_stat.block(b).variable(k)})
else else
fprintf('%s\n',endo_names{k}) fprintf('%s\n',endo_names{k})
end end
end end
if (~isoctave && matlab_ver_less_than('9.12')) || isempty(options.jacobian_tolerance) if (~isoctave && matlab_ver_less_than('9.12')) || isempty(options_.jacobian_tolerance)
neq = null(jacob'); %can sometimes fail neq = null(jacob'); %can sometimes fail
else else
neq = null(jacob',options.jacobian_tolerance); %can sometimes fail neq = null(jacob',options_.jacobian_tolerance); %can sometimes fail
end end
n_rel = size(neq,2); n_rel = size(neq,2);
for i = 1:n_rel for i = 1:n_rel
if n_rel > 1 if n_rel > 1
disp(['Relation ' int2str(i)]) disp(['Relation ' int2str(i)])
end end
disp('Colinear equations') disp('Collinear equations')
for j=1:10 for j=1:10
k = find(abs(neq(:,i)) > 10^-j); k = find(abs(neq(:,i)) > 10^-j);
if max(abs(jacob(k,:)'*neq(k,i))) < 1e-6 if max(abs(jacob(k,:)'*neq(k,i))) < 1e-6
break break
end end
end end
if options.block && ~options.bytecode if options_.block && ~options_.bytecode
disp(M.block_structure_stat.block(b).equation(k)) disp(M_.block_structure_stat.block(b).equation(k))
else else
disp(k') disp(k')
end end
...@@ -248,9 +255,9 @@ end ...@@ -248,9 +255,9 @@ end
if singularity_problem if singularity_problem
try try
options_check=options; options_check=options_;
options_check.noprint=1; options_check.noprint=1;
[eigenvalues_] = check(M, options_check, oo); [eigenvalues_] = check(M_, options_check, oo_);
if any(abs(abs(eigenvalues_)-1)<1e-6) if any(abs(abs(eigenvalues_)-1)<1e-6)
fprintf('MODEL_DIAGNOSTICS: The singularity seems to be (partly) caused by the presence of a unit root\n') fprintf('MODEL_DIAGNOSTICS: The singularity seems to be (partly) caused by the presence of a unit root\n')
fprintf('MODEL_DIAGNOSTICS: as the absolute value of one eigenvalue is in the range of +-1e-6 to 1.\n') fprintf('MODEL_DIAGNOSTICS: as the absolute value of one eigenvalue is in the range of +-1e-6 to 1.\n')
...@@ -265,51 +272,55 @@ if singularity_problem ...@@ -265,51 +272,55 @@ if singularity_problem
end end
%%check dynamic Jacobian %%check dynamic Jacobian
klen = M.maximum_lag + M.maximum_lead + 1; dyn_endo_ss = repmat(dr.ys, 3, 1);
exo_simul = [repmat(oo.exo_steady_state',klen,1) repmat(oo.exo_det_steady_state',klen,1)];
iyv = M.lead_lag_incidence';
iyv = iyv(:);
iyr0 = find(iyv) ;
it_ = M.maximum_lag + 1;
z = repmat(dr.ys,1,klen);
if options.order == 1 if options_.order == 1
if (options.bytecode) if (options_.bytecode)
[~, loc_dr] = bytecode('dynamic','evaluate', z,exo_simul, ... [~, loc_dr] = bytecode('dynamic','evaluate', M_, options_, z, exo_simul, ...
M.params, dr.ys, 1); M_.params, dr.ys, 1);
jacobia_ = [loc_dr.g1 loc_dr.g1_x loc_dr.g1_xd]; % TODO: simplify the following once bytecode MEX has been updated to sparse format
else g1 = zeros(M_.endo_nbr, 3*M_.endo_nbr+M_.exo_nbr+M_.exo_det_nbr);
[~,jacobia_] = feval([M.fname '.dynamic'],z(iyr0),exo_simul, ... if M_.maximum_endo_lag > 0
M.params, dr.ys, it_); g1(:, find(M_.lead_lag_incidence(M_.maximum_endo_lag, :))) = loc_dr.g1(:, 1:M_.nspred);
end end
elseif options.order >= 2 [~,icurr] = find(M_.lead_lag_incidence(M_.maximum_endo_lag+1, :));
if (options.bytecode) g1(:, M_.endo_nbr + icurr) = loc_dr.g1(:, M_.nspred+(1:length(icurr)));
[~, loc_dr] = bytecode('dynamic','evaluate', z,exo_simul, ... if M_.maximum_endo_lead > 0
M.params, dr.ys, 1); g1(:, 2*M_.endo_nbr + find(M_.lead_lag_incidence(M_.maximum_endo_lag+2, :))) = loc_dr.g1(:, M_.nspred+M_.endo_nbr+(1:M_.nsfwrd));
jacobia_ = [loc_dr.g1 loc_dr.g1_x]; end
g1(:, 3*M_.endo_nbr+(1:M_.exo_nbr)) = loc_dr.g1_x;
g1(:, 3*M_.endo_nbr+M_.exo_nbr+(1:M_.exo_det_nbr)) = loc_dr.g1_xd;
g1 = sparse(g1);
else else
[~,jacobia_,hessian1] = feval([M.fname '.dynamic'],z(iyr0),... g1 = feval([M_.fname '.sparse.dynamic_g1'], dyn_endo_ss, exo, M_.params, dr.ys, ...
exo_simul, ... M_.dynamic_g1_sparse_rowval, M_.dynamic_g1_sparse_colval, ...
M.params, dr.ys, it_); M_.dynamic_g1_sparse_colptr);
end
elseif options_.order >= 2
[g1, T_order, T] = feval([M_.fname '.sparse.dynamic_g1'], dyn_endo_ss, exo, M_.params, ...
dr.ys, M_.dynamic_g1_sparse_rowval, M_.dynamic_g1_sparse_colval, ...
M_.dynamic_g1_sparse_colptr);
if exist(['+' M_.fname '/+sparse/dynamic_g2.m'],"file") || exist(['+' M_.fname '/+sparse/dynamic_g2.' mexext],"file")
g2_v = feval([M_.fname '.sparse.dynamic_g2'], dyn_endo_ss, exo, M_.params, dr.ys, T_order, T);
end end
end end
if any(any(isinf(jacobia_) | isnan(jacobia_))) if any(any(isinf(g1) | isnan(g1)))
problem_dummy=1; problem_dummy=1;
[infrow,infcol]=find(isinf(jacobia_) | isnan(jacobia_)); [infrow,infcol]=find(isinf(g1) | isnan(g1));
fprintf('\nMODEL_DIAGNOSTICS: The Jacobian of the dynamic model contains Inf or NaN. The problem arises from: \n\n') fprintf('\nMODEL_DIAGNOSTICS: The Jacobian of the dynamic model contains Inf or NaN. The problem arises from: \n\n')
display_problematic_vars_Jacobian(infrow,infcol,M,dr.ys,'dynamic','MODEL_DIAGNOSTICS: ') display_problematic_vars_Jacobian(infrow,infcol,M_,dr.ys,'dynamic','MODEL_DIAGNOSTICS: ')
end end
if any(any(~isreal(jacobia_))) if any(any(~isreal(g1)))
[imagrow,imagcol]=find(abs(imag(jacobia_))>1e-15); [imagrow,imagcol]=find(abs(imag(g1))>1e-15);
if ~isempty(imagrow) if ~isempty(imagrow)
problem_dummy=1; problem_dummy=1;
fprintf('\nMODEL_DIAGNOSTICS: The Jacobian of the dynamic model contains imaginary parts. The problem arises from: \n\n') fprintf('\nMODEL_DIAGNOSTICS: The Jacobian of the dynamic model contains imaginary parts. The problem arises from: \n\n')
display_problematic_vars_Jacobian(imagrow,imagcol,M,dr.ys,'dynamic','MODEL_DIAGNOSTICS: ') display_problematic_vars_Jacobian(imagrow,imagcol,M_,dr.ys,'dynamic','MODEL_DIAGNOSTICS: ')
end end
end end
if exist('hessian1','var') if exist('g2_v','var')
if any(any(isinf(hessian1) | isnan(hessian1))) if any(any(isinf(g2_v) | isnan(g2_v)))
problem_dummy=1; problem_dummy=1;
fprintf('\nMODEL_DIAGNOSTICS: The Hessian of the dynamic model contains Inf or NaN.\n') fprintf('\nMODEL_DIAGNOSTICS: The Hessian of the dynamic model contains Inf or NaN.\n')
end end
......
...@@ -279,22 +279,11 @@ function print_line(names,var_index,lead_lag,M_) ...@@ -279,22 +279,11 @@ function print_line(names,var_index,lead_lag,M_)
end end
else else
aux_index=find([M_.aux_vars(:).endo_index]==var_index); aux_index=find([M_.aux_vars(:).endo_index]==var_index);
aux_type=M_.aux_vars(aux_index).type;
if ~isfield(M_.aux_vars(aux_index),'orig_lead_lag') || isempty(M_.aux_vars(aux_index).orig_lead_lag)
if ismember(aux_type,[1,3])
str = subst_auxvar(var_index, -1, M_);
elseif ismember(aux_type,[0,2])
str = subst_auxvar(var_index, 1, M_);
else
if lead_lag==0 if lead_lag==0
str = subst_auxvar(var_index, [], M_); str = subst_auxvar(var_index, [], M_);
else else
str = subst_auxvar(var_index, lead_lag, M_); str = subst_auxvar(var_index, lead_lag, M_);
end end
end
else
str = subst_auxvar(var_index, M_.aux_vars(aux_index).orig_lead_lag, M_);
end
aux_orig_expression=M_.aux_vars(aux_index).orig_expr; aux_orig_expression=M_.aux_vars(aux_index).orig_expr;
if isempty(aux_orig_expression) if isempty(aux_orig_expression)
fprintf('%s\n',str); fprintf('%s\n',str);
......
function [endogenousvariables, exogenousvariables] = model_inversion(constraints, ... function [endogenousvariables, exogenousvariables] = model_inversion(constraints, ...
exogenousvariables, ... exogenousvariables, ...
initialconditions, DynareModel, DynareOptions, DynareOutput) initialconditions, M_, options_, oo_)
% function [endogenousvariables, exogenousvariables] = model_inversion(constraints, ... % function [endogenousvariables, exogenousvariables] = model_inversion(constraints, ...
% exogenousvariables, ... % exogenousvariables, ...
% initialconditions, DynareModel, DynareOptions, DynareOutput) % initialconditions, M_, options_, oo_)
% INPUTS % INPUTS
% - constraints [dseries] with N constrained endogenous variables from t1 to t2. % - constraints [dseries] with N constrained endogenous variables from t1 to t2.
% - exogenousvariables [dseries] with Q exogenous variables. % - exogenousvariables [dseries] with Q exogenous variables.
% - initialconditions [dseries] with M endogenous variables starting before t1 (M initialcond must contain at least the state variables). % - initialconditions [dseries] with M endogenous variables starting before t1 (M initialcond must contain at least the state variables).
% - DynareModel [struct] M_, Dynare global structure containing informations related to the model. % - M_ [struct] Dynare global structure containing informations related to the model.
% - DynareOptions [struct] options_, Dynare global structure containing all the options. % - options_ [struct] Dynare global structure containing all the options.
% - DynareOutput [struct] oo_, Dynare global structure containing all the options. % - oo_ [struct] Dynare global structure containing all the options.
% %
% OUTPUTS % OUTPUTS
% - endogenousvariables [dseries] % - endogenousvariables [dseries]
...@@ -18,7 +18,7 @@ function [endogenousvariables, exogenousvariables] = model_inversion(constraints ...@@ -18,7 +18,7 @@ function [endogenousvariables, exogenousvariables] = model_inversion(constraints
% %
% REMARKS % REMARKS
% Copyright © 2018-2021 Dynare Team % Copyright © 2018-2023 Dynare Team
% %
% This file is part of Dynare. % This file is part of Dynare.
% %
...@@ -51,7 +51,7 @@ if ~isempty(initialconditions) && ~isdseries(initialconditions) ...@@ -51,7 +51,7 @@ if ~isempty(initialconditions) && ~isdseries(initialconditions)
error('model_inversion: Third input argument must be a dseries object!') error('model_inversion: Third input argument must be a dseries object!')
end end
if ~isstruct(DynareModel) if ~isstruct(M_)
error('model_inversion: Last input argument must be structures (M_)!') error('model_inversion: Last input argument must be structures (M_)!')
end end
...@@ -71,29 +71,29 @@ if exogenousvariables.vobs>constraints.vobs ...@@ -71,29 +71,29 @@ if exogenousvariables.vobs>constraints.vobs
observed_exogenous_variables_flag = true; observed_exogenous_variables_flag = true;
end end
if DynareModel.maximum_lag if M_.maximum_lag
% Add auxiliary variables in initialconditions object. % Add auxiliary variables in initialconditions object.
initialconditions = checkdatabase(initialconditions, DynareModel, true, false); initialconditions = checkdatabase(initialconditions, M_, true, false);
end end
% Get the list of endogenous and exogenous variables. % Get the list of endogenous and exogenous variables.
endo_names = DynareModel.endo_names; endo_names = M_.endo_names;
exo_names = DynareModel.exo_names; exo_names = M_.exo_names;
exogenousvariables = exogenousvariables{exo_names{:}}; exogenousvariables = exogenousvariables{exo_names{:}};
% Use specidalized routine if the model is backward looking. % Use specidalized routine if the model is backward looking.
if ~DynareModel.maximum_lead if ~M_.maximum_lead
if DynareModel.maximum_lag if M_.maximum_lag
[endogenousvariables, exogenousvariables] = ... [endogenousvariables, exogenousvariables] = ...
backward_model_inversion(constraints, exogenousvariables, initialconditions, ... backward_model_inversion(constraints, exogenousvariables, initialconditions, ...
endo_names, exo_names, freeinnovations, ... endo_names, exo_names, freeinnovations, ...
DynareModel, DynareOptions, DynareOutput); M_, options_, oo_);
else else
[endogenousvariables, exogenousvariables] = ... [endogenousvariables, exogenousvariables] = ...
static_model_inversion(constraints, exogenousvariables, ... static_model_inversion(constraints, exogenousvariables, ...
endo_names, exo_names, freeinnovations, ... endo_names, exo_names, freeinnovations, ...
DynareModel, DynareOptions, DynareOutput); M_, options_, oo_);
end end
return return
end end
......