Merge branch 'dynare-complex_resid'

Ref. !2155
kalman-mex
Sébastien Villemot 2023-09-20 10:09:07 +02:00
commit 44f307ce45
No known key found for this signature in database
GPG Key ID: 2CECE9350ECEBE4A
7 changed files with 69 additions and 45 deletions

View File

@ -52,6 +52,7 @@ options_.lyapunov_complex_threshold = 1e-15;
options_.solve_algo = 4;
options_.solve_tolf = eps^(1/3);
options_.solve_tolx = eps^(2/3);
options_.solve_randomize_initial_guess = true;
options_.trust_region_initial_step_bound_factor = 1;
options_.dr_display_tol=1e-6;
options_.minimal_workspace = false;

View File

@ -51,7 +51,7 @@ if any(imag(oo_.steady_state))
end
if options_.steadystate_flag
[oo_.steady_state,M_.params,info] = ...
[oo_.steady_state,M_.params] = ...
evaluate_steady_state(oo_.steady_state,[oo_.exo_steady_state; oo_.exo_det_steady_state],M_,options_,false);
end
@ -91,6 +91,8 @@ if nargout == 0
first_eq = 1;
last_eq = M_.orig_endo_nbr;
end
disp_format_tags_real=sprintf('Equation number %%%uu: %%-%us: %%14.6f\n',length(num2str(M_.eq_nbr)),size(strvcat(tags(:,3)),2)+1);
disp_format_tags_complex=sprintf('Equation number %%%uu: %%-%us: real: %%14.6f, imaginary: %%g\n',length(num2str(M_.eq_nbr)),size(strvcat(tags(:,3)),2)+1);
for i=first_eq:last_eq
if abs(z(i)) < options_.solve_tolf/100
tmp = 0;
@ -104,11 +106,17 @@ if nargout == 0
end
if ~(non_zero && tmp == 0)
if ~istag || length(ind) == 0
disp(['Equation number ' int2str(i) ' : ' num2str(tmp)])
if ~isreal(z)
fprintf('Equation number %u: %g (imaginary part: %g)\n', i, real(tmp), imag(tmp));
else
fprintf('Equation number %u: %g\n', i, tmp);
end
else
t1 = tg( ind , 2 );
s = cell2mat(t1);
disp( ['Equation number ', int2str(i) ,' : ', num2str(tmp) ,' : ' s])
if ~isreal(z)
fprintf(disp_format_tags_complex, i, tg{ind , 2}, real(tmp), imag(tmp) );
else
fprintf(disp_format_tags_real, i, tg{ind , 2}, tmp);
end
end
end
end

View File

@ -59,7 +59,6 @@ else
end
% checking initial values
% TODO We should have an option to deactivate the randomization.
if jacobian_flag
[fvec, fjac] = feval(f, x, varargin{:});
wrong_initial_guess_flag = false;
@ -70,32 +69,38 @@ if jacobian_flag
errorcode = -11;
return;
end
disp_verbose('Randomize initial guess...', options.verbosity)
% Let's try random numbers for the variables initialized with the default value.
wrong_initial_guess_flag = true;
% First try with positive numbers.
tentative_number = 0;
while wrong_initial_guess_flag && tentative_number<=in0*10
tentative_number = tentative_number+1;
x(idx) = rand(in0, 1)*10;
[fvec, fjac] = feval(f, x, varargin{:});
wrong_initial_guess_flag = ~all(isfinite(fvec)) || any(isinf(fjac(:))) || any(isnan((fjac(:))));
end
% If all previous attempts failed, try with real numbers.
tentative_number = 0;
while wrong_initial_guess_flag && tentative_number<=in0*10
tentative_number = tentative_number+1;
x(idx) = randn(in0, 1)*10;
[fvec, fjac] = feval(f, x, varargin{:});
wrong_initial_guess_flag = ~all(isfinite(fvec)) || any(isinf(fjac(:))) || any(isnan((fjac(:))));
end
% Last tentative, ff all previous attempts failed, try with negative numbers.
tentative_number = 0;
while wrong_initial_guess_flag && tentative_number<=in0*10
tentative_number = tentative_number+1;
x(idx) = -rand(in0, 1)*10;
[fvec, fjac] = feval(f, x, varargin{:});
wrong_initial_guess_flag = ~all(isfinite(fvec)) || any(isinf(fjac(:))) || any(isnan((fjac(:))));
if options.solve_randomize_initial_guess
if any(~isreal(fvec)) || any(~isreal(fjac(:)))
disp_verbose('dynare_solve: starting value results in complex values. Randomize initial guess...', options.verbosity)
else
disp_verbose('dynare_solve: starting value results in nonfinite/NaN value. Randomize initial guess...', options.verbosity)
end
% Let's try random numbers for the variables initialized with the default value.
wrong_initial_guess_flag = true;
% First try with positive numbers.
tentative_number = 0;
while wrong_initial_guess_flag && tentative_number<=in0*10
tentative_number = tentative_number+1;
x(idx) = rand(in0, 1)*10;
[fvec, fjac] = feval(f, x, varargin{:});
wrong_initial_guess_flag = ~all(isfinite(fvec)) || any(isinf(fjac(:))) || any(isnan((fjac(:)))) || any(~isreal(fvec)) || any(~isreal(fjac(:)));
end
% If all previous attempts failed, try with real numbers.
tentative_number = 0;
while wrong_initial_guess_flag && tentative_number<=in0*10
tentative_number = tentative_number+1;
x(idx) = randn(in0, 1)*10;
[fvec, fjac] = feval(f, x, varargin{:});
wrong_initial_guess_flag = ~all(isfinite(fvec)) || any(isinf(fjac(:))) || any(isnan((fjac(:)))) || any(~isreal(fvec)) || any(~isreal(fjac(:)));
end
% Last tentative, ff all previous attempts failed, try with negative numbers.
tentative_number = 0;
while wrong_initial_guess_flag && tentative_number<=in0*10
tentative_number = tentative_number+1;
x(idx) = -rand(in0, 1)*10;
[fvec, fjac] = feval(f, x, varargin{:});
wrong_initial_guess_flag = ~all(isfinite(fvec)) || any(isinf(fjac(:))) || any(isnan((fjac(:)))) || any(~isreal(fvec)) || any(~isreal(fjac(:)));
end
end
end
else

View File

@ -293,22 +293,19 @@ elseif ~options.bytecode && ~options.block
options.mcppath.lb = lb;
options.mcppath.ub = ub;
end
[ys,check,fvec] = dynare_solve(@static_mcp_problem,...
[ys,check,fvec, fjac, errorcode] = dynare_solve(@static_mcp_problem,...
ys_init,...
options.steady.maxit, options.solve_tolf, options.solve_tolx, ...
options, exo_ss, params,...
M.endo_nbr, static_resid, static_g1, ...
M.static_g1_sparse_rowval, M.static_g1_sparse_colval, M.static_g1_sparse_colptr, eq_index);
else
[ys, check] = dynare_solve(@static_problem, ys_init, ...
[ys, check, fvec, fjac, errorcode] = dynare_solve(@static_problem, ys_init, ...
options.steady.maxit, options.solve_tolf, options.solve_tolx, ...
options, exo_ss, params, M.endo_nbr, static_resid, static_g1, ...
M.static_g1_sparse_rowval, M.static_g1_sparse_colval, M.static_g1_sparse_colptr);
end
if check && options.debug
[ys, check, fvec, fjac, errorcode] = dynare_solve(@static_problem, ys_init, ...
options.steady.maxit, options.solve_tolf, options.solve_tolx, ...
options, exo_ss, params, M.endo_nbr, static_resid, static_g1, M.static_g1_sparse_rowval, M.static_g1_sparse_colval, M.static_g1_sparse_colptr);
dprintf('Nonlinear solver routine returned errorcode=%i.', errorcode)
skipline()
[infrow,infcol]=find(isinf(fjac) | isnan(fjac));
@ -480,10 +477,14 @@ if M.static_and_dynamic_models_differ
end
if ~isreal(ys)
info(1) = 21;
info(2) = sum(imag(ys).^2);
ys = real(ys);
if sum(imag(ys).^2) < 1e-7
ys=real(ys);
else
info(1) = 21;
info(2) = sum(imag(ys).^2);
ys = real(ys);
return
end
end
if ~isempty(find(isnan(ys)))

View File

@ -67,9 +67,9 @@ switch info(1)
message = 'The steadystate file did not compute the steady state';
case 20
if DynareOptions.linear
message = sprintf('Impossible to find the steady state (the sum of square residuals of the static equations is %5.4f). Either the model doesn''t have a steady state or there are an infinity of steady states Check whether your model is truly linear or whether there is a mistake in linearization.', info(2));
message = sprintf('Impossible to find the steady state (the sum of squared residuals of the static equations is %5.4f). Either the model doesn''t have a steady state or there are an infinity of steady states. Check whether your model is truly linear or whether there is a mistake in linearization.', info(2));
else
message = sprintf('Impossible to find the steady state (the sum of square residuals of the static equations is %5.4f). Either the model doesn''t have a steady state, there are an infinity of steady states, or the guess values are too far from the solution', info(2));
message = sprintf('Impossible to find the steady state (the sum of squared residuals of the static equations is %5.4f). Either the model doesn''t have a steady state, there are an infinity of steady states, or the guess values are too far from the solution', info(2));
end
case 21
message = sprintf('The steady state is complex (the sum of square residuals of imaginary parts of the steady state is %5.4f)', info(2));

View File

@ -94,9 +94,18 @@ else
end
end
if options_.debug
fprintf('\nThe steady state computation failed. It terminated with the following values:\n')
fprintf('\nsteady: The steady state computation failed. It terminated with the following values:\n')
if ~isreal(oo_.steady_state)
format_string=sprintf('%%-%us= %%g%%+gi\n',size(strvcat(M_.endo_names),2)+1);
else
format_string=sprintf('%%-%us= %%14.6f\n',size(strvcat(M_.endo_names),2)+1);
end
for i=1:M_.orig_endo_nbr
fprintf('%s \t\t %g\n', M_.endo_names{i}, oo_.steady_state(i));
if ~isreal(oo_.steady_state)
fprintf(format_string, M_.endo_names{i}, real(oo_.steady_state(i)),imag(oo_.steady_state(i)));
else
fprintf(format_string, M_.endo_names{i}, oo_.steady_state(i));
end
end
end
print_info(info,options_.noprint, options_);

@ -1 +1 @@
Subproject commit d95c2f2bd81a4e8fc56552d19e62ac37b8d08a81
Subproject commit bdb5cdacf22019c75db38dcdfe696f2978ab31f9