function homotopy3(values, step_nbr)
Implements homotopy (mode 3) for steady-state computation.
Tries first the most extreme values. If it fails to compute the steady
state, the interval between initial and desired values is divided by two
for each parameter. Every time that it is impossible to find a steady
state, the previous interval is divided by two. When one succeed to find
a steady state, the previous interval is multiplied by two.
INPUTS
values: a matrix with 4 columns, representing the content of
homotopy_setup block, with one variable per line.
Column 1 is variable type (1 for exogenous, 2 for
exogenous deterministic, 4 for parameters)
Column 2 is symbol integer identifier.
Column 3 is initial value, and column 4 is final value.
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
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
none0001 function [M,oo,info,ip,ix,ixd]=homotopy3(values, step_nbr, M, options, oo) 0002 % function homotopy3(values, step_nbr) 0003 % 0004 % Implements homotopy (mode 3) for steady-state computation. 0005 % Tries first the most extreme values. If it fails to compute the steady 0006 % state, the interval between initial and desired values is divided by two 0007 % for each parameter. Every time that it is impossible to find a steady 0008 % state, the previous interval is divided by two. When one succeed to find 0009 % a steady state, the previous interval is multiplied by two. 0010 % 0011 % INPUTS 0012 % values: a matrix with 4 columns, representing the content of 0013 % homotopy_setup block, with one variable per line. 0014 % Column 1 is variable type (1 for exogenous, 2 for 0015 % exogenous deterministic, 4 for parameters) 0016 % Column 2 is symbol integer identifier. 0017 % Column 3 is initial value, and column 4 is final value. 0018 % Column 3 can contain NaNs, in which case previous 0019 % initialization of variable will be used as initial value. 0020 % step_nbr: maximum number of steps to try before aborting 0021 % M struct of model parameters 0022 % options struct of options 0023 % oo struct of outputs 0024 % 0025 % OUTPUTS 0026 % M struct of model parameters 0027 % oo struct of outputs 0028 % info return status 0: OK, 1: failed 0029 % ip index of parameters 0030 % ix index of exogenous variables 0031 % ixp index of exogenous deterministic variables 0032 % 0033 % SPECIAL REQUIREMENTS 0034 % none 0035 0036 % Copyright (C) 2008-2009 Dynare Team 0037 % 0038 % This file is part of Dynare. 0039 % 0040 % Dynare is free software: you can redistribute it and/or modify 0041 % it under the terms of the GNU General Public License as published by 0042 % the Free Software Foundation, either version 3 of the License, or 0043 % (at your option) any later version. 0044 % 0045 % Dynare is distributed in the hope that it will be useful, 0046 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0047 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0048 % GNU General Public License for more details. 0049 % 0050 % You should have received a copy of the GNU General Public License 0051 % along with Dynare. If not, see <http://www.gnu.org/licenses/>. 0052 0053 info = []; 0054 tol = 1e-8; 0055 0056 nv = size(values,1); 0057 0058 ip = find(values(:,1) == 4); % Parameters 0059 ix = find(values(:,1) == 1); % Exogenous 0060 ixd = find(values(:,1) == 2); % Exogenous deterministic 0061 0062 if length([ip; ix; ixd]) ~= nv 0063 error('HOMOTOPY mode 3: incorrect variable types specified') 0064 end 0065 0066 % Construct vector of starting values, using previously initialized values 0067 % when initial value has not been given in homotopy_setup block 0068 last_values = values(:,3); 0069 ipn = find(values(:,1) == 4 & isnan(last_values)); 0070 last_values(ipn) = M.params(values(ipn, 2)); 0071 ixn = find(values(:,1) == 1 & isnan(last_values)); 0072 last_values(ixn) = oo.exo_steady_state(values(ixn, 2)); 0073 ixdn = find(values(:,1) == 2 & isnan(last_values)); 0074 last_values(ixdn) = oo.exo_det_steady_state(values(ixdn, 2)); 0075 0076 targetvalues = values(:,4); 0077 0078 %if min(abs(targetvalues - last_values)) < tol 0079 % error('HOMOTOPY mode 3: distance between initial and final values should be at least %e for all variables', tol) 0080 %end 0081 iplus = find(targetvalues > last_values); 0082 iminus = find(targetvalues < last_values); 0083 0084 curvalues = last_values; 0085 inc = (targetvalues-last_values)/2; 0086 kplus = []; 0087 kminus = []; 0088 last_values = []; 0089 0090 disp('HOMOTOPY mode 3: launching solver at initial point...') 0091 0092 iter = 1; 0093 while iter <= step_nbr 0094 0095 M.params(values(ip,2)) = curvalues(ip); 0096 oo.exo_steady_state(values(ix,2)) = curvalues(ix); 0097 oo.exo_det_steady_state(values(ixd,2)) = curvalues(ixd); 0098 0099 old_ss = oo.steady_state; 0100 0101 [steady_state,params,info] = steady_(M,options,oo); 0102 if info(1) == 0 0103 oo.steady_state = steady_state; 0104 M.params = params; 0105 if length([kplus; kminus]) == nv 0106 return 0107 end 0108 if iter == 1 0109 disp('HOMOTOPY mode 3: successful step, now jumping to final point...') 0110 else 0111 disp('HOMOTOPY mode 3: successful step, now multiplying increment by 2...') 0112 end 0113 last_values = curvalues; 0114 old_params = params; 0115 old_exo_steady_state = oo.exo_steady_state; 0116 old_exo_det_steady_state = oo.exo_det_steady_state; 0117 inc = 2*inc; 0118 elseif iter == 1 0119 error('HOMOTOPY mode 3: can''t solve the model at 1st iteration') 0120 else 0121 disp('HOMOTOPY mode 3: failed step, now dividing increment by 2...') 0122 inc = inc/2; 0123 oo.steady_state = old_ss; 0124 end 0125 0126 curvalues = last_values + inc; 0127 kplus = find(curvalues(iplus) >= targetvalues(iplus)); 0128 curvalues(iplus(kplus)) = targetvalues(iplus(kplus)); 0129 kminus = find(curvalues(iminus) <= targetvalues(iminus)); 0130 curvalues(iminus(kminus)) = targetvalues(iminus(kminus)); 0131 0132 if max(abs(inc)) < tol 0133 disp('HOMOTOPY mode 3: failed, increment has become too small') 0134 M.params = old_params; 0135 oo.exo_steady_state = old_exo_steady_state; 0136 oo.exo_det_steady_state = old_exo_det_steady_state; 0137 return 0138 end 0139 0140 iter = iter + 1; 0141 end 0142 disp('HOMOTOPY mode 3: failed, maximum iterations reached') 0143 M.params = old_params; 0144 oo.exo_steady_state = old_exo_steady_state; 0145 oo.exo_det_steady_state = old_exo_det_steady_state;