stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

v4 matlab+preprocessor: the three homotopy modes should now work 

git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@1769 ac1d8469-bf42-47a9-8791-bf33cf982152
time-shift
sebastien 2008-03-31 16:19:16 +00:00
parent 06b6429b04
commit 367948e4db
13 changed files with 279 additions and 125 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

BIN
matlab/dynare_m.exe
View File

Binary file not shown.

100
matlab/homotopy1.m
View File

@ -1,66 +1,62 @@
% homotopy1 implements a homotopy method with a fixed number of intervals
function homotopy1(values, step_nbr)
% function homotopy1(values, step_nbr)
%
% Implements homotopy (mode 1) for steady-state computation.
% The multi-dimensional vector going from the set of initial values
% to the set of final values is divided in as many sub-vectors as
% there are steps, and the problem is solved as many times.
%
% 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.
% step_nbr: number of steps for homotopy
%
% OUTPUTS
% none
%
% SPECIAL REQUIREMENTS
% none
%
% part of DYNARE, copyright Dynare Team (2008)
% Gnu Public License.
function homotopy1(params,exo,exodet, step_nbr)
global M_ oo_ options_
options_.jacobian_flag = 1;
np = size(params,1);
ip = zeros(np,1);
vp = zeros(np,step_nbr+1);
nv = size(values, 1);
ip = find(values(:,1) == 4); % Parameters
ix = find(values(:,1) == 1); % Exogenous
ixd = find(values(:,1) == 2); % Exogenous deterministic
nx = size(exo,1);
ix = zeros(nx,1);
vx = zeros(nx,step_nbr+1);
nxd = size(exodet,1);
ixd = zeros(nxd,1);
vxd = zeros(nxd,step_nbr+1);
for i = 1:np
temp1 = strmatch(params{i,1},M_.param_names,'exact');
if isempty(temp1)
error(['HOMOTOPY: unknown parameter name: ' params{i,1}])
end
ip(i) = temp1;
vp(i,:) = params{i,2}:(params{i,3}-params{i,2})/step_nbr:params{i,3};
if length([ip, ix, ixd]) ~= nv
error('HOMOTOPY: incorrect variable types specified')
end
if any(values(:,3) == values(:,4))
error('HOMOTOPY: initial and final values should be different')
end
for i = 1:nx
temp1 = strmatch(exo{i,1},M_.exo_names,'exact');
if isempty(temp1)
error(['HOMOTOPY: unknown exogenous variable name: ' exo{i,1}])
end
ix(i) = temp1;
vp(i,:) = exo{i,2}:(exo{i,3}-exo{i,2})/step_nbr:exo{i,3};
points = zeros(nv, step_nbr+1);
for i = 1:nv
points(i,:) = values(i,3):(values(i,4)-values(i,3))/step_nbr:values(i,4);
end
for i = 1:nxd
temp1 = strmatch(exodet{i,1},M_.exodet_names,'exact');
if isempty(temp1)
error(['HOMOTOPY: unknown deterministic exogenous name: ' exodet{i,1}])
end
ixd(i) = temp1;
vxd(i,:) = exodet{i,2}:(exodet{i,3}-exodet{i,2})/step_nbr:exodet{i,3};
end
for i=1:step_nbr+1
M_.params(ip) = vp(:,i);
for j=1:np
assignin('base',params{j,1},vp(j,i));
end
oo_.exo_steady_state(ix) = vx(:,i);
oo_.exo_det_steady_state(ixd) = vxd(:,i);
for i=1:step_nbr+1
M_.params(values(ip,2)) = points(ip,i);
oo_.exo_steady_state(values(ix,2)) = points(ix,i);
oo_.exo_det_steady_state(values(ixd,2)) = points(ixd,i);
[oo_.steady_state,check] = dynare_solve([M_.fname '_static'],...
oo_.steady_state,...
options_.jacobian_flag, ...
[oo_.exo_steady_state; ...
oo_.exo_det_steady_state]);
oo_.steady_state,...
options_.jacobian_flag, ...
[oo_.exo_steady_state; ...
oo_.exo_det_steady_state]);
if check
error('HOMOTOPY didn''t succeed')
end
steady_;
end
end

74
matlab/homotopy2.m Executable file
View File

@ -0,0 +1,74 @@
function homotopy2(values, step_nbr)
% function homotopy2(values, step_nbr)
%
% Implements homotopy (mode 2) for steady-state computation.
% Only one parameter/exogenous is changed at a time.
% Computation jumps to next variable only when current variable has been
% brought to its final value.
% Variables are processed in the order in which they appear in "values".
% The problem is solved var_nbr*step_nbr times.
%
% 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.
% step_nbr: number of steps for homotopy
%
% OUTPUTS
% none
%
% SPECIAL REQUIREMENTS
% none
%
% part of DYNARE, copyright Dynare Team (2008)
% Gnu Public License.
global M_ oo_ options_
nv = size(values, 1);
% Initialize all variables with initial value
for i = 1:nv
switch values(i,1)
case 1
oo_.exo_steady_state(values(i,2)) = values(i,3);
case 2
oo_.exo_det_steady_state(values(i,2)) = values(i,3);
case 4
M_.params(values(i,2)) = values(i,3);
otherwise
error('HOMOTOPY: incorrect variable types specified')
end
end
if any(values(:,3) == values(:,4))
error('HOMOTOPY: initial and final values should be different')
end
% Actually do the homotopy
for i = 1:nv
for v = values(i,3):(values(i,4)-values(i,3))/step_nbr:values(i,4)
switch values(i,1)
case 1
oo_.exo_steady_state(values(i,2)) = v;
case 2
oo_.exo_det_steady_state(values(i,2)) = v;
case 4
M_.params(values(i,2)) = v;
end
[oo_.steady_state,check] = dynare_solve([M_.fname '_static'],...
oo_.steady_state,...
options_.jacobian_flag, ...
[oo_.exo_steady_state; ...
oo_.exo_det_steady_state]);
if check
error('HOMOTOPY didn''t succeed')
end
end
end

120
matlab/homotopy3.m
View File

@ -1,65 +1,91 @@
% homotopy3 implements a homotopy method that reduces the step as much as necessary
function homotopy3(params,exo,exodet, step_nbr)
function homotopy3(values, step_nbr)
% 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.
% step_nbr: maximum number of steps to try before aborting
%
% OUTPUTS
% none
%
% SPECIAL REQUIREMENTS
% none
%
% part of DYNARE, copyright Dynare Team (2008)
% Gnu Public License.
global M_ oo_ options_
options_.jacobian_flag = 1;
np = length(param_names);
ip = zeros(np,1);
oldvalues = zeros(np,1);
iplus = [];
iminus = [];
for i = 1:np
temp1 = strmatch(param_names{i},M_.param_names,'exact');
if isempty(temp1)
error(['HOMOTOPY: unknown parameter name: ' param_names{i}])
end
ip(i) = temp1;
oldvalues(i) = param_values{i}(1);
targetvalues(i) = param_values{i}(2);
if targetvalues(i) > oldvalues(i)
iplus = [iplus i];
else
iminus = [iminus i];
end
end
tol = 1e-8;
iter = 1
maxiter = 500;
values = oldvalues;
nv = size(values,1);
ip = find(values(:,1) == 4); % Parameters
ix = find(values(:,1) == 1); % Exogenous
ixd = find(values(:,1) == 2); % Exogenous deterministic
if length([ip, ix, ixd]) ~= nv
error('HOMOTOPY: incorrect variable types specified')
end
oldvalues = values(:,3);
targetvalues = values(:,4);
if min(abs(targetvalues - oldvalues)) < tol
error('HOMOTOPY: distance between initial and final values should be at least %e for all variables', tol)
end
iplus = find(targetvalues > oldvalues);
iminus = find(targetvalues < oldvalues);
curvalues = oldvalues;
inc = (targetvalues-oldvalues)/2;
k = [];
kplus = [];
kminus = [];
old_ss = oo_.steady_state;
while iter < maxiter
for j=1:np
M_.params(ip(j)) = values(j,i);
assignin('base',param_names{1},values(j,1));
end
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);
[oo_.steady_state,check] = dynare_solve([M_.fname '_static'],...
oo_.steady_state,...
options_.jacobian_flag, ...
[oo_.exo_steady_state; ...
oo_.exo_det_steady_state]);
oo_.steady_state,...
options_.jacobian_flag, ...
[oo_.exo_steady_state; ...
oo_.exo_det_steady_state]);
if check
inc = inc/2;
oo_.steady_state = old_ss;
else
if length(k) == np
return
if length([kplus, kminus]) == nv
return
end
oldvalues = values;
oldvalues = curvalues;
inc = 2*inc;
end
values = oldvalues + inc
k = find(values(iplus) > targetvalues(iplus));
values(k) = targetvalues(k);
k = find(values(iminus) < targetvalues(iminus));
values(k) = targetvalues(k);
values
if max(abs(inc)) < 1e-8
curvalues = oldvalues + inc;
kplus = find(curvalues(iplus) >= targetvalues(iplus));
curvalues(kplus) = targetvalues(kplus);
kminus = find(curvalues(iminus) <= targetvalues(iminus));
curvalues(kminus) = targetvalues(kminus);
if max(abs(inc)) < tol
error('HOMOTOPY didn''t succeed')
end
iter = iter + 1;
end
error('HOMOTOPY didn''t succeed')
error('HOMOTOPY didn''t succeed')

15
matlab/steady.m
View File

@ -23,7 +23,20 @@ function steady()
if exist([M_.fname '_steadystate.m'])
options_.steadystate_flag = 1;
end
if options_.steadystate_flag && options_.homotopy_mode
error('STEADY: Can''t use homotopy when providing a steady state external file');
end
switch options_.homotopy_mode
case 1
homotopy1(options_.homotopy_values, options_.homotopy_steps);
case 2
homotopy2(options_.homotopy_values, options_.homotopy_steps);
case 3
homotopy3(options_.homotopy_values, options_.homotopy_steps);
end
steady_;
disp(' ')

30
preprocessor/NumericalInitialization.cc
View File

@ -154,7 +154,7 @@ HistValStatement::writeOutput(ostream &output, const string &basename) const
}
HomotopyStatement::HomotopyStatement(const homotopy_values_type &homotopy_values_arg,
const SymbolTable &symbol_table_arg) :
const SymbolTable &symbol_table_arg) :
homotopy_values(homotopy_values_arg),
symbol_table(symbol_table_arg)
{
@ -163,11 +163,10 @@ HomotopyStatement::HomotopyStatement(const homotopy_values_type &homotopy_values
void
HomotopyStatement::writeOutput(ostream &output, const string &basename) const
{
output << interfaces::comment() << "\n" << interfaces::comment() << "HOMOTOPY_SETUP instructions\n"
<< interfaces::comment() << "\n";
output << "homotopy_param = {};" << endl;
output << "homotopy_exo = {};" << endl;
output << "homotopy_exodet = {};" << endl;
output << interfaces::comment() << endl
<< interfaces::comment() << "HOMOTOPY_SETUP instructions" << endl
<< interfaces::comment() << endl
<< "options_.homotopy_values = [];" << endl;
for(homotopy_values_type::const_iterator it = homotopy_values.begin();
it != homotopy_values.end(); it++)
@ -176,24 +175,13 @@ HomotopyStatement::writeOutput(ostream &output, const string &basename) const
const NodeID expression1 = it->second.first;
const NodeID expression2 = it->second.second;
Type type = symbol_table.getType(name);
const Type type = symbol_table.getType(name);
const int id = symbol_table.getID(name) + 1;
if (type == eParameter)
output << "homotopy_param = vertcat(homotopy_param,{ '" << name << "', ";
else if (type == eExogenous)
output << "homotopy_exo = vertcat(homotopy_exo,{ '" << name << "', ";
else if (type != eExogenousDet)
output << "homotopy_exodet = vertcat(homotopy_exodet,{ '" << name << "', ";
output << "options_.homotopy_values = vertcat(options_.homotopy_values, [ " << type << ", " << id << ", ";
expression1->writeOutput(output);
output << ", ";
expression2->writeOutput(output);
output << "});" << endl;
output << "]);" << endl;
}
output << "options_.homotopy_param = homotopy_param;" << endl;
output << "options_.homotopy_exo = homotopy_exo;" << endl;
output << "options_.homotopy_exodet = homotopy_exodet;" << endl;
}

3
preprocessor/ParsingDriver.cc
View File

@ -350,8 +350,7 @@ ParsingDriver::homotopy_val(string *name, NodeID val1, NodeID val2)
if (homotopy_values.find(*name) != homotopy_values.end())
error("homotopy_val: " + *name +" declared twice");
pair<NodeID, NodeID> expressions(val1, val2);
homotopy_values[*name] = expressions;
homotopy_values[*name] = make_pair(val1, val2);
delete name;
}

2
preprocessor/include/CodeInterpreter.hh
View File

@ -58,7 +58,7 @@ const int EVALUATE_FOREWARD_R=8;
const int EVALUATE_BACKWARD_R=9;
//! Enumeration of possible symbol types
/*! Warning: do not to change the order of the enumeration, it matters for VariableTable (at least ensure that eEndogenous is the first one) */
/*! Warning: do not to change existing values: the order matters for VariableTable (at least ensure that eEndogenous is the first one), and the values matter for homotopy_setup command */
enum Type
{
eEndogenous = 0, //!< Endogenous

2
preprocessor/include/NumericalInitialization.hh
View File

@ -105,7 +105,7 @@ private:
const SymbolTable &symbol_table;
public:
HomotopyStatement(const homotopy_values_type &homotopy_values_arg,
const SymbolTable &symbol_table_arg);
const SymbolTable &symbol_table_arg);
virtual void writeOutput(ostream &output, const string &basename) const;
};
#endif

13
tests/homotopy/common.mod Executable file
View File

@ -0,0 +1,13 @@
var c k;
varexo x;
parameters alph gam delt bet aa;
alph=0.5;
gam=0.5;
delt=0.02;
aa=0.5;
model;
c + k - aa*x*k(-1)^alph - (1-delt)*k(-1);
c^(-gam) - (1+bet)^(-1)*(aa*alph*x(+1)*k^(alph-1) + 1 - delt)*c(+1)^(-gam);
end;

13
tests/homotopy/ramst_final.mod Executable file
View File

@ -0,0 +1,13 @@
// Computes the steady state which should be arrived at in ramst_homotopy.mod
@include "common.mod"
bet = 0.1;
initval;
x = 2;
k = ((delt+bet)/(1.0*aa*alph))^(1/(alph-1));
c = aa*k^alph-delt*k;
end;
steady;

19
tests/homotopy/ramst_homotopy.mod Executable file
View File

@ -0,0 +1,19 @@
// Contains the test for homotopy.
// The values in initval block are obtained from ramst_initial.mod
// Result of the computation should be the same than in ramst_final.mod
@include "common.mod"
initval;
k = 12.75;
c = 1.5;
end;
homotopy_setup;
bet, 0.05, 0.1;
x, 1, 2;
end;
steady(homotopy_mode = 1, homotopy_steps = 50);
//steady(homotopy_mode = 2, homotopy_steps = 50);
//steady(homotopy_mode = 3, homotopy_steps = 50);

13
tests/homotopy/ramst_initial.mod Executable file
View File

@ -0,0 +1,13 @@
// Computes the steady state used in initval block of ramst_homotopy.mod
@include "common.mod"
bet=0.05;
initval;
x = 1;
k = ((delt+bet)/(1.0*aa*alph))^(1/(alph-1));
c = aa*k^alph-delt*k;
end;
steady;