From d064db34ab2fb134808daa524c6bcc5a3ad12b21 Mon Sep 17 00:00:00 2001 From: Michel Juillard Date: Sun, 20 May 2012 08:47:54 +0200 Subject: [PATCH] Cleaning homotopy1. Removing globals from homotopy3. Now it returns last successful values when it fails --- matlab/homotopy1.m | 5 +-- matlab/homotopy3.m | 77 ++++++++++++++++++++++++++++------------------ matlab/steady.m | 12 ++++---- 3 files changed, 54 insertions(+), 40 deletions(-) diff --git a/matlab/homotopy1.m b/matlab/homotopy1.m index 4b7009905..1d5d203ca 100644 --- a/matlab/homotopy1.m +++ b/matlab/homotopy1.m @@ -1,4 +1,4 @@ -function [M,oo,info,last_values,ip,ix,ixd] = homotopy1(values, step_nbr, M, options, oo) +function [M,oo,info,ip,ix,ixd] = homotopy1(values, step_nbr, M, options, oo) % function homotopy1(values, step_nbr) % % Implements homotopy (mode 1) for steady-state computation. @@ -23,7 +23,6 @@ function [M,oo,info,last_values,ip,ix,ixd] = homotopy1(values, step_nbr, M, opti % OUTPUTS % M struct of model parameters % oo struct of outputs -% last_values last set of values fir which a solution was found % ip index of parameters % ix index of exogenous variables % ixp index of exogenous deterministic variables @@ -48,7 +47,6 @@ function [M,oo,info,last_values,ip,ix,ixd] = homotopy1(values, step_nbr, M, opti % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . -last_values = []; nv = size(values, 1); ip = find(values(:,1) == 4); % Parameters @@ -96,7 +94,6 @@ for i=1:step_nbr+1 M.params = old_params; oo.exo_steady_state = old_exo; oo.exo_det_steady_state = old_exo_det; - last_values = points(:,i-1); break end end diff --git a/matlab/homotopy3.m b/matlab/homotopy3.m index 5359680bf..5f77c0677 100644 --- a/matlab/homotopy3.m +++ b/matlab/homotopy3.m @@ -1,4 +1,4 @@ -function homotopy3(values, step_nbr) +function [M,oo,info,ip,ix,ixd]=homotopy3(values, step_nbr, M, options, oo) % function homotopy3(values, step_nbr) % % Implements homotopy (mode 3) for steady-state computation. @@ -18,9 +18,17 @@ function homotopy3(values, step_nbr) % Column 3 can contain NaNs, in which case previous % initialization of variable will be used as initial value. % step_nbr: maximum number of steps to try before aborting +% M struct of model parameters +% options struct of options +% oo struct of outputs % % OUTPUTS -% none +% M struct of model parameters +% oo struct of outputs +% info return status 0: OK, 1: failed +% ip index of parameters +% ix index of exogenous variables +% ixp index of exogenous deterministic variables % % SPECIAL REQUIREMENTS % none @@ -42,8 +50,7 @@ function homotopy3(values, step_nbr) % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . -global M_ oo_ options_ - +info = []; tol = 1e-8; nv = size(values,1); @@ -58,41 +65,43 @@ end % Construct vector of starting values, using previously initialized values % when initial value has not been given in homotopy_setup block -oldvalues = values(:,3); -ipn = find(values(:,1) == 4 & isnan(oldvalues)); -oldvalues(ipn) = M_.params(values(ipn, 2)); -ixn = find(values(:,1) == 1 & isnan(oldvalues)); -oldvalues(ixn) = oo_.exo_steady_state(values(ixn, 2)); -ixdn = find(values(:,1) == 2 & isnan(oldvalues)); -oldvalues(ixdn) = oo_.exo_det_steady_state(values(ixdn, 2)); +last_values = values(:,3); +ipn = find(values(:,1) == 4 & isnan(last_values)); +last_values(ipn) = M.params(values(ipn, 2)); +ixn = find(values(:,1) == 1 & isnan(last_values)); +last_values(ixn) = oo.exo_steady_state(values(ixn, 2)); +ixdn = find(values(:,1) == 2 & isnan(last_values)); +last_values(ixdn) = oo.exo_det_steady_state(values(ixdn, 2)); targetvalues = values(:,4); -if min(abs(targetvalues - oldvalues)) < tol +if min(abs(targetvalues - last_values)) < tol error('HOMOTOPY mode 3: distance between initial and final values should be at least %e for all variables', tol) end -iplus = find(targetvalues > oldvalues); -iminus = find(targetvalues < oldvalues); +iplus = find(targetvalues > last_values); +iminus = find(targetvalues < last_values); -curvalues = oldvalues; -inc = (targetvalues-oldvalues)/2; +curvalues = last_values; +inc = (targetvalues-last_values)/2; kplus = []; kminus = []; +last_values = []; disp('HOMOTOPY mode 3: launching solver at initial point...') iter = 1; while iter < step_nbr - M_.params(values(ip,2)) = curvalues(ip); - oo_.exo_steady_state(values(ix,2)) = curvalues(ix); - oo_.exo_det_steady_state(values(ixd,2)) = curvalues(ixd); + M.params(values(ip,2)) = curvalues(ip); + oo.exo_steady_state(values(ix,2)) = curvalues(ix); + oo.exo_det_steady_state(values(ixd,2)) = curvalues(ixd); - old_ss = oo_.steady_state; + old_ss = oo.steady_state; - try - oo_.steady_state = steady_(M_,options_,oo_); - + [steady_state,params,info] = steady_(M,options,oo); + if info(1) == 0 + oo.steady_state = steady_state; + M.params = params; if length([kplus; kminus]) == nv return end @@ -101,25 +110,33 @@ while iter < step_nbr else disp('HOMOTOPY mode 3: successful step, now multiplying increment by 2...') end - oldvalues = curvalues; + last_values = curvalues; + old_params = params; + old_exo_steady_state = oo.exo_steady_state; + old_exo_det_steady_state = oo.exo_det_steady_state; inc = 2*inc; - catch E - disp(E.message) + elseif iter == 1 + error('HOMOTOPY mode 3: can''t solve the model at 1st iteration') + else disp('HOMOTOPY mode 3: failed step, now dividing increment by 2...') inc = inc/2; - oo_.steady_state = old_ss; + oo.steady_state = old_ss; end - curvalues = oldvalues + inc; + curvalues = last_values + inc; kplus = find(curvalues(iplus) >= targetvalues(iplus)); curvalues(iplus(kplus)) = targetvalues(iplus(kplus)); kminus = find(curvalues(iminus) <= targetvalues(iminus)); curvalues(iminus(kminus)) = targetvalues(iminus(kminus)); if max(abs(inc)) < tol - error('HOMOTOPY mode 3: failed, increment has become too small') + disp('HOMOTOPY mode 3: failed, increment has become too small') + M.params = old_params; + oo.exo_steady_state = old_exo_steady_state; + oo.exo_det_steady_state = old_exo_det_steady_state; + return end iter = iter + 1; end -error('HOMOTOPY mode 3: failed, maximum iterations reached') +disp('HOMOTOPY mode 3: failed, maximum iterations reached') diff --git a/matlab/steady.m b/matlab/steady.m index ea1e7508b..7bd3c800d 100644 --- a/matlab/steady.m +++ b/matlab/steady.m @@ -46,29 +46,29 @@ M_.Sigma_e = zeros(size(Sigma_e)); info = 0; switch options_.homotopy_mode case 1 - [M_,oo_,info,last_values,ip,ix,ixd] = homotopy1(options_.homotopy_values, options_.homotopy_steps,M_,options_,oo_); + [M_,oo_,info,ip,ix,ixd] = homotopy1(options_.homotopy_values,options_.homotopy_steps,M_,options_,oo_); case 2 homotopy2(options_.homotopy_values, options_.homotopy_steps); case 3 - homotopy3(options_.homotopy_values, options_.homotopy_steps); + [M_,oo_,info,ip,ix,ixd] = homotopy3(options_.homotopy_values,options_.homotopy_steps,M_,options_,oo_); end if info(1) hv = options_.homotopy_values; disp(' ') - disp('WARNING: homotopy step was not comleted') + disp('WARNING: homotopy step was not completed') disp('The last values for which a solution was found are:') for i=1:length(ip) disp(sprintf('%12s %12.6f',M_.param_names(hv(ip(i),2),:), ... - last_values(ip(i)))) + M_.params(hv(ip(i))))) end for i=1:length(ix) disp(sprintf('%12s %12.6f',M_.exo_names(hv(ix(i),2),:), ... - last_values(ix(i)))) + oo_.exo_steady_state(hv(ix(i))))) end for i=1:length(ixd) disp(sprintf('%12s %12.6f',M_.exo_det_names(hv(ixd(i),2),:), ... - last_values(ixd(i)))) + oo_.exo_det_steady_state(hv(ixd(i))))) end if options_.homotopy_force_continue