Add the possibility to specify bounds for the BGP solver.

bgp-dev
Stéphane Adjemian (Ryûk) 2022-10-23 17:16:36 +02:00
parent 872de4e781
commit 572db78143
Signed by: stepan
GPG Key ID: 295C1FE89E17EB3C
6 changed files with 382 additions and 55 deletions

101
matlab/+bgp/+transform/in.m Normal file
View File

@ -0,0 +1,101 @@
function Variables = in(DynareModel, constraints, variables, id, lb, ub)
if ~iscell(constraints)
error('Second argument must be a one-dimensional cell array.')
end
if iskeyvalue(constraints)
if nargin<4
id = zeros(length(constraints)/2,1);
j = 1;
for i=1:2:length(constraints)
if isempty(regexp(constraints{i}, 'G\(\w+\)')) i
tmp = find(strcmp(constraints{i}, DynareModel.endo_names));
else
tmp = find(strcmp(constraints{i}(3:end-1), DynareModel.endo_names));
if ~isempty(tmp)
tmp = tmp+DynareModel.endo_nbr;
end
end
if isempty(tmp)
error('Unknown variable (%s).', constraints{i})
else
id(j) = tmp;
end
j = j+1;
end
end
Variables = NaN(2*DynareModel.endo_nbr, 1);
jj = 1:length(Variables);
ll = setdiff(jj, id);
Variables(ll) = variables;
values = NaN(length(constraints)/2, 1);
j = 1;
for i=2:2:length(constraints)
values(j) = constraints{i};
j = j+1;
end
Variables(id) = values;
return
end
if iskeyvalues(constraints)
halfline.right = @(x, a) a+exp(x);
halfline.left = @(x, a) a-exp(x);
interval = @(x, a, b) a+(b-a)/(1+exp(x));
if nargin<4
id = zeros(length(constraints)/3,1);
lb = [-Inf(DynareModel.endo_nbr, 1); zeros(DynareModel.endo_nbr, 1)];
ub = Inf(2*DynareModel.endo_nbr, 1);
j = 1;
for i=1:3:length(constraints)
if isempty(regexp(constraints{i}, 'G\(\w+\)'))
tmp = find(strcmp(constraints{i}, DynareModel.endo_names));
else
tmp = find(strcmp(constraints{i}(3:end-1), DynareModel.endo_names));
if ~isempty(tmp)
tmp = tmp + DynareModel.endo_nbr;
end
end
assert(~isempty(tmp), 'Unknown variable (%s).', constraints{i})
lb(tmp) = constraints{i+1};
ub(tmp) = constraints{i+2};
if abs(lb(tmp)-ub(tmp))<1e-6
id(j) = tmp;
end
j = j+1;
end
id = id(id>0);
end
Variables = NaN(2*DynareModel.endo_nbr, 1);
jj = 1:length(Variables);
ll = setdiff(jj, id);
Variables(ll) = variables;
values = NaN(length(id), 1);
j = 1;
for i=1:length(id)
if id(i)>DynareModel.endo_nbr
j = find(strcmp(sprintf('G(%s)', DynareModel.endo_names{id(i)-DynareModel.endo_nbr}), constraints));
else
j = find(strcmp(sprintf('%s', DynareModel.endo_names{id(i)}), constraints));
end
assert(~isempty(j), 'Unknown variable (indexing issue).')
values(i) = constraints{j+1};
end
Variables(id) = values;
for i=1:length(Variables)
if isfinite(lb(i)) && isinf(ub(i))
Variables(i) = halfline.right(Variables(i), lb(i));
elseif isinf(lb(i)) && isfinite(ub(i))
Variables(i) = halfline.left(Variables(i), ub(i));
elseif isfinite(lb(i)) && isfinite(ub(i)) && ub(i)-lb(i)>1e-6
Variables(i) = interval(Variables(i), lb(i), ub(i));
end
end
return
end
error('Second argument must be a cell array of key/value pairs or key/values triplets.')
end

View File

@ -0,0 +1,78 @@
function [variables, id, lb, ub] = out(DynareModel, constraints, variables)
if ~iscell(constraints)
error('Second argument must be a one-dimensional cell array.')
end
if iskeyvalue(constraints)
% Calibrated variables.
id = zeros(length(constraints)/2,1);
j = 1;
for i=1:2:length(constraints)
if isempty(regexp(constraints{i}, 'G\(\w+\)'))
tmp = find(strcmp(constraints{i}, DynareModel.endo_names));
else
tmp = DynareModel.endo_nbr+find(strcmp(constraints{i}(3:end-1), DynareModel.endo_names));
end
assert(~isempty(tmp), 'Unknown variable (%s).', constraints{i})
id(j) = tmp;
j = j+1;
end
variables(id) = [];
lb = []; ub = [];
return
end
if iskeyvalues(constraints)
% Boundaries on variables and calibrated variables
halfline.right = @(x, a) log(x-a);
halfline.x = @(left, a) log(a-x);
interval = @(x, a, b) log(x-a)-log(b-x);
id = zeros(length(constraints)/3,1);
lb = [-Inf(DynareModel.endo_nbr, 1); zeros(DynareModel.endo_nbr, 1)];
ub = Inf(2*DynareModel.endo_nbr, 1);
j = 1;
for i=1:3:length(constraints)
if isempty(regexp(constraints{i}, 'G\(\w+\)'))
tmp = find(strcmp(constraints{i}, DynareModel.endo_names));
else
tmp = find(strcmp(constraints{i}(3:end-1), DynareModel.endo_names));
if isempty(tmp)
tmp = tmp + DynareModel.endo_nbr;
end
end
assert(~isempty(tmp), 'Unknown variable (%s).', constraints{i})
lb(tmp) = constraints{i+1};
ub(tmp) = constraints{i+2};
if abs(lb(tmp)-ub(tmp))<1e-6
id(j) = tmp;
end
j = j+1;
end
id = id(id>0);
id
for i=1:2*DynareModel.endo_nbr
if ~ismember(i, id)
if i>DynareModel.endo_nbr
assert(lb(i)<variables(i) && ub(i)>variables(i), 'Initial condition not consistent with declared boundaries (G(%s)).', DynareModel.endo_names{i-DynareModel.endo_nbr})
else
assert(lb(i)<variables(i) && ub(i)>variables(i), 'Initial condition not consistent with declared boundaries (%s)', DynareModel.endo_names{i})
end
if isfinite(lb(i)) && isinf(ub(i))
variables(i) = halfline.right(variables(i), lb(i));
elseif isinf(lb(i)) && isfinite(ub(i))
variables(i) = halfline.left(variables(i), ub(i));
elseif isfinite(lb(i)) && isfinite(ub(i)) && lb(i)<ub(i)
variables(i) = interval(variables(i), lb(i), ub(i));
end
end
end
length(variables)
variables(id) = [];
length(variables)
return
end
error('Second argument must be a cell array of key/value pairs or key/values triplets.')
end

View File

@ -15,7 +15,7 @@ function write(DynareModel, varargin)
% REMARKS
% - The trends are assumed to be multiplicative.
% Copyright © 2019 Dynare Team
% Copyright © 2019-2022 Dynare Team
%
% This file is part of Dynare.
%
@ -51,6 +51,42 @@ function write(DynareModel, varargin)
error('Wrong input arguments (%s is not an endogenous variable).', varargin{1})
end
end
elseif iskeyvalues(varargin)
idy = NaN(length(varargin)/3, 1);
vly = NaN(length(varargin)/3, 1);
idg = NaN(length(varargin)/3, 1);
vlg = NaN(length(varargin)/3, 1);
lb = [-Inf(DynareModel.endo_nbr, 1); zeros(DynareModel.endo_nbr, 1)];
ub = Inf(2*DynareModel.endo_nbr, 1);
prn = 0;
for i=1:3:length(varargin)
prn = prn+1;
if ismember(varargin{i}, DynareModel.endo_names)
j = find(strcmp(varargin{i}, DynareModel.endo_names));
lb(j) = varargin{i+1};
ub(j) = varargin{i+2};
if ub(j)-lb(j)<0
error('An upper bound cannot be smaller than a lower bound (%s).', varargin{i})
end
if ub(j)-lb(j)<1e-6
idy(prn) = find(strcmp(varargin{i}, DynareModel.endo_names));
vly(prn) = varargin{i+1};
end
elseif ~isempty(regexp(varargin{i}, 'G\(\w+\)')) && ismember(varargin{i}(3:end-1), DynareModel.endo_names)
j = DynareModel.endo_nbr+find(strcmp(varargin{i}(3:end-1), DynareModel.endo_names));
lb(j) = varargin{i+1};
ub(j) = varargin{i+2};
if ub(j)-lb(j)<0
error('An upper bound cannot be smaller than a lower bound (%s).', varargin{i})
end
if ub(j)-lb(j)<1e-6
idg(prn) = find(strcmp(varargin{i}, DynareModel.endo_names));
vlg(prn) = varargin{i+1};
end
else
error('Wrong input arguments (%s is not an endogenous variable).', varargin{1})
end
end
else
error('Wrong input arguments.')
end
@ -354,56 +390,101 @@ function write(DynareModel, varargin)
fprintf(fid,'end\n');
fclose(fid);
if nargin>1
if iskeyvalue(varargin)
% Create function in mod namespace.
fid = fopen(sprintf('+%s/bgpfunhardcon.m', DynareModel.fname), 'w');
% Write header.
fprintf(fid, 'function [F, JAC] = bgpfunhardcon(z)\n\n');
fprintf(fid, '%% This file has been generated by dynare (%s).\n\n', datestr(now));
j1 = setdiff(transpose(1:DynareModel.endo_nbr), idy(~isnan(idy)));
j2 = setdiff(transpose(1:DynareModel.endo_nbr), idg(~isnan(idg)));
tmp = sprintf(' %u', [j1; DynareModel.endo_nbr+j2]);
fprintf(fid, 'id = [%s];\n\n', strtrim(tmp));
fprintf(fid, 'Z = NaN(%u, 1);\n\n', 2*DynareModel.endo_nbr);
fprintf(fid, 'Z(id) = z;\n');
for i=1:length(idy)
if ~isnan(idy(i))
fprintf(fid, 'Z(%u) = %s;\n', idy(i), num2str(vly(i)));
end
if ~isnan(idg(i))
fprintf(fid, 'Z(%u) = %s;\n', length(idy)+idg(i), num2str(vlg(i)));
end
end
skipline(1, fid)
fprintf(fid, '[F, JAC] = %s.bgpfun(Z);\n\n', DynareModel.fname);
for i=1:length(idy)
if ~isnan(idy(i))
fprintf(fid, 'JAC(:,%u) = [];\n', idy(i));
end
if ~isnan(idg(i))
fprintf(fid, 'JAC(:,%u) = [];\n', length(idy)+idg(i));
end
end
fclose(fid);
elseif iskeyvalues(varargin)
% Create function in mod namespace.
fid = fopen(sprintf('+%s/bgpfuncon.m', DynareModel.fname), 'w');
% Write header.
fprintf(fid, 'function [F, JAC] = bgpfuncon(z)\n\n');
fprintf(fid, '%% This file has been generated by dynare (%s).\n\n', datestr(now));
%fprintf(fid, 'keyboard\n\n');
j1 = setdiff(transpose(1:DynareModel.endo_nbr), idy(~isnan(idy)));
j2 = setdiff(transpose(1:DynareModel.endo_nbr), idg(~isnan(idg)));
jj = [j1; DynareModel.endo_nbr+j2];
fprintf(fid, 'Z = NaN(%u, 1);\n\n', 2*DynareModel.endo_nbr);
for i=1:length(jj)
j = jj(i);
if isinf(lb(j)) && isinf(ub(j))
fprintf(fid, 'Z(%u) = z(%u);\n', j, i);
elseif isfinite(lb(j)) && isinf(ub(j))
fprintf(fid, 'Z(%u) = %s+exp(z(%u));\n', j, num2str(lb(j), 9), i);
elseif isinf(lb(j)) && isfinite(ub(j))
fprintf(fid, 'Z(%u) = %s-exp(z(%u));\n', j, num2str(ub(j), 9), i);
elseif isfinite(lb(j)) && isfinite(ub(j))
fprintf(fid, 'Z(%u) = %s+%s/(1+exp(-z(%u)));\n', j, num2str(lb(j), 9), num2str(ub(j)-lb(j), 9), i);
end
end
skipline(1, fid)
for i=1:length(idy)
if ~isnan(idy(i))
fprintf(fid, 'Z(%u) = %s;\n', idy(i), num2str(vly(i)));
end
if ~isnan(idg(i))
fprintf(fid, 'Z(%u) = %s;\n', length(idy)+idg(i), num2str(vlg(i)));
end
end
skipline(1, fid)
fprintf(fid, '[F, JAC] = %s.bgpfun(Z);\n\n', DynareModel.fname);
%
for i=1:length(jj)
j = jj(i);
if isinf(lb(j)) && isinf(ub(j))
% Do nothing (no transformation of the input)
elseif isfinite(lb(j)) && isinf(ub(j))
fprintf(fid, 'JAC(:,%u) = JAC(:,%u)*exp(z(%u));\n', j, j, i);
elseif isinf(lb(j)) && isfinite(ub(j))
fprintf(fid, 'JAC(:,%u) = -JAC(:,%u)*exp(z(%u));\n', j, j, i);
elseif isfinite(lb(j)) && isfinite(ub(j))
fprintf(fid, 'tmp = 1/(1+exp(-z(%u)));\n', i);
fprintf(fid, 'JAC(:,%u) = JAC(:,%u)*tmp*(1-tmp)*%s;\n', j, j, num2str(ub(j)-lb(j), 9));
end
end
%
skipline(1, fid)
for i=1:length(idy)
if ~isnan(idy(i))
fprintf(fid, 'JAC(:,%u) = [];\n', idy(i));
end
if ~isnan(idg(i))
fprintf(fid, 'JAC(:,%u) = [];\n', length(idy)+idg(i));
end
end
fclose(fid);
end
if nargin>1 && iskeyvalue(varargin)
% Create function in mod namespace.
fid = fopen(sprintf('+%s/bgpfunhardcon.m', DynareModel.fname), 'w');
% Write header.
fprintf(fid, 'function [F, JAC] = bgpfunhardcon(z)\n\n');
fprintf(fid, '%% This file has been generated by dynare (%s).\n\n', datestr(now));
j1 = setdiff(transpose(1:DynareModel.endo_nbr), idy(~isnan(idy)));
j2 = setdiff(transpose(1:DynareModel.endo_nbr), idg(~isnan(idg)));
tmp = sprintf(' %u', [j1; DynareModel.endo_nbr+j2]);
fprintf(fid, 'id = [%s];\n\n', strtrim(tmp));
fprintf(fid, 'Z = NaN(%u, 1);\n\n', 2*DynareModel.endo_nbr);
fprintf(fid, 'Z(id) = z;\n');
for i=1:length(idy)
if ~isnan(idy(i))
fprintf(fid, 'Z(%u) = %s;\n', idy(i), num2str(vly(i)));
end
if ~isnan(idg(i))
fprintf(fid, 'Z(%u) = %s;\n', length(idy)+idg(i), num2str(vlg(i)));
end
end
skipline(1, fid)
fprintf(fid, '[F, JAC] = %s.bgpfun(Z);\n\n', DynareModel.fname);
for i=1:length(idy)
if ~isnan(idy(i))
fprintf(fid, 'JAC(:,%u) = [];\n', idy(i));
end
if ~isnan(idg(i))
fprintf(fid, 'JAC(:,%u) = [];\n', length(idy)+idg(i));
end
end
fclose(fid);
end
rehash()
end
function is = iskeyvalue(rowcellarray)
is = true;
if ~iscell(rowcellarray) || ~isrow(rowcellarray) || ~iseven(length(rowcellarray))
is = false;
return
end
if any(~cellfun(@ischar, rowcellarray(1:2:length(rowcellarray))))
is = false;
return
end
if any(~cellfun(@(x) isnumeric(x) & isscalar(x), rowcellarray(2:2:length(rowcellarray))))
is = false;
return
end
end

View File

@ -0,0 +1,15 @@
function is = iskeyvalue(rowcellarray)
is = true;
if ~iscell(rowcellarray) || ~isrow(rowcellarray) || ~iseven(length(rowcellarray))
is = false;
return
end
if any(~cellfun(@ischar, rowcellarray(1:2:length(rowcellarray))))
is = false;
return
end
if any(~cellfun(@(x) isnumeric(x) & isscalar(x), rowcellarray(2:2:length(rowcellarray))))
is = false;
return
end
end

View File

@ -0,0 +1,19 @@
function is = iskeyvalues(rowcellarray)
is = true;
if ~iscell(rowcellarray) || ~isrow(rowcellarray) || ~iseven(length(rowcellarray))
is = false;
return
end
if any(~cellfun(@ischar, rowcellarray(1:3:length(rowcellarray))))
is = false;
return
end
if any(~cellfun(@(x) isnumeric(x) & isscalar(x), rowcellarray(2:3:length(rowcellarray))))
is = false;
return
end
if any(~cellfun(@(x) isnumeric(x) & isscalar(x), rowcellarray(3:3:length(rowcellarray))))
is = false;
return
end
end

View File

@ -1,3 +1,4 @@
// --+ options: json=compute +--
/*
* This file is a modified version of 'fs2000.mod'.
*
@ -61,12 +62,38 @@ end;
verbatim;
HardCodedLevels = {'k', 10, 'M', 1};
Constraints = {'M', 1, Inf, ...
'A', 1, Inf, ...
'k', 10, 10, ...
'c', 0, Inf, ...
'y', 0, Inf, ...
'P', 10, 10, ...
'W', 0, Inf, ...
'l', 0, Inf, ...
'd', 0, Inf, ...
'n', 0, Inf, ...
'R', 0, Inf, ...
'gM', 1, Inf, ...
'gA', 1, Inf, ...
'G(M)', 1, 1.1, ...
'G(A)', 1, 1.1, ...
'G(k)', 1, 1.1, ...
'G(c)', 1, 1.1, ...
'G(y)', 1, 1.1, ...
'G(P)', 1, 1.1, ...
'G(W)', 1, 1.1, ...
xb 'G(l)', 1, 1.1, ...
'G(d)', 1, 1.1, ...
'G(n)', .9, 1.1, ...
'G(R)', .9, 1.1, ...
'G(gM)', .9, 1.1, ...
'G(gA)', .9, 1.1, ...
};
bgp.write(M_, HardCodedLevels{:});
bgp.write(M_, Constraints{:});
y = ones(M_.endo_nbr,1)+(rand(M_.endo_nbr,1))*.1;
G = ones(M_.endo_nbr,1);% 1+(rand(M_.endo_nbr,1)-.5)*.1;
y = ones(M_.endo_nbr,1) + rand(M_.endo_nbr,1)*.1;
G = ones(M_.endo_nbr,1) + rand(M_.endo_nbr,1)*.1;
y = [y;G];
if isoctave
@ -78,11 +105,17 @@ verbatim;
options = optimoptions('fsolve','Display','iter','Algorithm','levenberg-marquardt','MaxFunctionEvaluations',10000000,'MaxIterations',1000000,'SpecifyObjectiveGradient',true,'FunctionTolerance',1e-4,'StepTolerance',1e-6);
end
[y, id] = bgp.hdl.out(M_, HardCodedLevels, y);
[y, id, lb, ub] = bgp.transform.out(M_, Constraints, y);
[y, fval, exitflag] = fsolve(@fs2000.bgpfunhardcon, y, options);
[y, fval, exitflag] = fsolve(@fs2000.bgpfuncon, y, options);
y = bgp.hdl.in(M_, HardCodedLevels, y, id);
y(1:M_.orig_endo_nbr-2)
y(M_.endo_nbr-2+(1:M_.orig_endo_nbr))
y = bgp.transform.in(M_, Constraints, y, id, lb, ub);
y(1:M_.orig_endo_nbr)
y(M_.endo_nbr+(1:M_.orig_endo_nbr))
if exitflag<1
error('Solution not found')