From 53b1910201dda28aefe2a9f10a88f8e0bd9e65f5 Mon Sep 17 00:00:00 2001 From: Michel Juillard Date: Thu, 5 Oct 2017 14:19:24 +0200 Subject: [PATCH] adding missing occbin functions --- matlab/occbin/evaluate_model.m | 19 ++++ matlab/occbin/get_coef.m | 26 ++++++ .../occbin/get_complementarity_conditions.m | 72 +++++++++++++++ .../get_occbin_complementarity_conditions.m | 75 ++++++++++++++++ matlab/occbin/get_occbin_constraints.m | 88 +++++++++++++++++++ matlab/occbin/get_residuals.m | 10 +++ matlab/occbin/occbin.m | 19 ++++ matlab/occbin/test_constraint.m | 16 ++++ 8 files changed, 325 insertions(+) create mode 100644 matlab/occbin/evaluate_model.m create mode 100644 matlab/occbin/get_coef.m create mode 100644 matlab/occbin/get_complementarity_conditions.m create mode 100644 matlab/occbin/get_occbin_complementarity_conditions.m create mode 100644 matlab/occbin/get_occbin_constraints.m create mode 100644 matlab/occbin/get_residuals.m create mode 100644 matlab/occbin/occbin.m create mode 100644 matlab/occbin/test_constraint.m diff --git a/matlab/occbin/evaluate_model.m b/matlab/occbin/evaluate_model.m new file mode 100644 index 000000000..f01d58a08 --- /dev/null +++ b/matlab/occbin/evaluate_model.m @@ -0,0 +1,19 @@ +function [r,g1,g2,g3] = evaluate_model(z,x,M,ss) + +ll = M.lead_lag_incidence'; +y = z(find(ll(:))); + +switch nargout + case 1 + r = feval([M.fname '_dynamic'],y,x, ... + M.params, ss, 1); + case 2 + [r,g1] = feval([M.fname '_dynamic'],y,x, ... + M.params, ss, 1); + case 3 + [r,g1,g2] = feval([M.fname '_dynamic'],y,x, ... + M.params, ss, 1); + case 4 + [r,g1,g2,g3] = feval([M.fname '_dynamic'],y,x, ... + M.params, ss, 1); +end \ No newline at end of file diff --git a/matlab/occbin/get_coef.m b/matlab/occbin/get_coef.m new file mode 100644 index 000000000..f2dbde5a6 --- /dev/null +++ b/matlab/occbin/get_coef.m @@ -0,0 +1,26 @@ +function [coef_y,coef_u] = get_coef(jacobian,M) + +ll = M.lead_lag_incidence; +endo_nbr = M.endo_nbr; +coef_y = zeros(endo_nbr,3*endo_nbr); +coef_u = zeros(endo_nbr,M.exo_nbr); + +if M.maximum_lag > 0 + [junk,c1,c2] = find(ll(1,:)); + coef_y(:,c1) = jacobian(:,c2); + [junk,c1,c2] = find(ll(2,:)); + coef_y(:,c1+endo_nbr) = jacobian(:,c2); + if M.maximum_lead > 0 + [junk,c1,c2] = find(ll(3,:)); + coef_y(:,c1+2*endo_nbr) = jacobian(:,c2); + end +else + [junk,c1,c2] = find(ll(1,:)); + coef_y(:,c1+endo_nbr) = jacobian(:,c2); + if M.maximum_lead > 0 + [junk,c1,c2] = find(ll(2,:)); + coef_y(:,c1+2*endo_nbr) = jacobian(:,c2); + end +end + +coef_u = jacobian(:,max(ll(end,:))+1:end); \ No newline at end of file diff --git a/matlab/occbin/get_complementarity_conditions.m b/matlab/occbin/get_complementarity_conditions.m new file mode 100644 index 000000000..cf5a40866 --- /dev/null +++ b/matlab/occbin/get_complementarity_conditions.m @@ -0,0 +1,72 @@ +function [lb,ub,eq_index] = get_complementarity_conditions(M,ramsey_policy) + +% Copyright (C) 2014 Dynare Team +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see . + +ub = inf(M.endo_nbr,1); +lb = -ub; +eq_index = (1:M.endo_nbr)'; +if ramsey_policy + if isfield(M,'ramsey_model_constraints') + rc = M.ramsey_model_constraints; + for i = 1:length(rc) + switch rc{i}{2} + case {'>','>='} + lb(rc{i}{1}) = eval(rc{i}{3}); + case {'<','<='} + ub(rc{i}{1}) = eval(rc{i}{3}); + otherwise + error('Wrong operator in get_complementarity_conditions') + end + eq_index(i) = 1; + end + end +end + +etags = M.equations_tags; +for i=1:size(etags,1) + if strcmp(etags{i,2},'mcp') + str = etags{i,3}; + kop = strfind(etags{i,3},'<'); + if ~isempty(kop) + k = find(strcmp(strtrim(str(1:kop-1)),cellstr(M.endo_names))); + if isempty(k) + error(sprintf(['Complementarity condition %s: variable %s is ' ... + 'not recognized',etags{i,3},b{1}])) + end + ub(k) = str2num(str(kop+1:end)); + eq_index(etags{i,1}) = k; + eq_index(k) = etags{i,1}; + else + kop = strfind(etags{i,3},'>'); + if ~isempty(kop) + k = find(strcmp(strtrim(str(1:kop-1)),cellstr(M.endo_names))); + if isempty(k) + error(sprintf(['Complementarity condition %s: variable %s is ' ... + 'not recognized',etags{i},b{1}])) + end + lb(k) = str2num(str(kop+1:end)); + eq_index(etags{i,1}) = k; + eq_index(k) = etags{i,1}; + else + error(sprintf(['Complementarity condition %s can''t be ' ... + 'parsed'],etags{i,3})) + end + end + end +end + diff --git a/matlab/occbin/get_occbin_complementarity_conditions.m b/matlab/occbin/get_occbin_complementarity_conditions.m new file mode 100644 index 000000000..29c865587 --- /dev/null +++ b/matlab/occbin/get_occbin_complementarity_conditions.m @@ -0,0 +1,75 @@ +function [ivar,ieq,lb,ub] = get_occbin_complementarity_conditions(M,ramsey_policy) + +% Copyright (C) 2015 Dynare Team +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see . + +nrow = 1; +if ramsey_policy + if isfield(M,'ramsey_model_constraints') + rc = M.ramsey_model_constraints; + for i = 1:length(rc) + switch rc{i}{2} + case {'>','>='} + ivar(nrow) = rc{i}{1}; + ieq(nrow) = rc{i}{1}; + lb(nrow) = eval(rc{i}{3}); + case {'<','<='} + ivar(nrow) = rc{i}{1}; + ieq(nrow) = rc{i}{1}; + ub(nrow) = eval(rc{i}{3}); + otherwise + error('Wrong operator in get_complementarity_conditions') + end + nrow = nrow + 1; + end + end +end + +etags = M.equations_tags; +for i=1:size(etags,1) + if strcmp(etags{i,2},'mcp') + str = etags{i,3}; + kop = strfind(etags{i,3},'<'); + if ~isempty(kop) + k = find(strcmp(strtrim(str(1:kop-1)),cellstr(M.endo_names))); + if isempty(k) + error(sprintf(['Complementarity condition %s: variable %s is ' ... + 'not recognized',etags{i,3},b{1}])) + end + ivar(nrow) = k; + ieq(nrow) = etags{i,1}; + ub(nrow) = eval(str(kop+1:end)); + else + kop = strfind(etags{i,3},'>'); + if ~isempty(kop) + k = find(strcmp(strtrim(str(1:kop-1)),cellstr(M.endo_names))); + if isempty(k) + error(sprintf(['Complementarity condition %s: variable %s is ' ... + 'not recognized',etags{i},b{1}])) + end + ivar(nrow) = k; + ieq(nrow) = etags{i,1}; + lb(k) = eval(str(kop+1:end)); + else + error(sprintf(['Complementarity condition %s can''t be ' ... + 'parsed'],etags{i,3})) + end + end + nrow = nrow + 1; + end +end + diff --git a/matlab/occbin/get_occbin_constraints.m b/matlab/occbin/get_occbin_constraints.m new file mode 100644 index 000000000..6723f4542 --- /dev/null +++ b/matlab/occbin/get_occbin_constraints.m @@ -0,0 +1,88 @@ +function [i_base,i_alt,c_base,c_alt] = get_occbin_constraints(M,steady_state,ramsey_policy) + +% Copyright (C) 2015 Dynare Team +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see . + +nrow = 1; +if ramsey_policy + if isfield(M,'ramsey_model_constraints') + rc = M.ramsey_model_constraints; + for i = 1:length(rc) + switch rc{i}{2} + case {'>','>='} + ivar(nrow) = rc{i}{1}; + ieq(nrow) = rc{i}{1}; + lb(nrow) = eval(rc{i}{3}); + case {'<','<='} + ivar(nrow) = rc{i}{1}; + ieq(nrow) = rc{i}{1}; + ub(nrow) = eval(rc{i}{3}); + otherwise + error('Wrong operator in get_complementarity_conditions') + end + nrow = nrow + 1; + end + end +end + +i_base = {}; +i_alt = {}; +etags = M.equations_tags; +m = 1; +base = true; +for i=1:size(etags,1) + [iv,boundary,operator] = parse_constraint(etags{i,3},M.endo_names,M.params,M.param_names); + if strcmp(etags{i,2},'OCCBIN') + if base + i_alt{m} = 1:M.eq_nbr; + i_alt{m}(etags{i,1}) = []; + c_base{m,1} = etags{i,1}; + c_base{m,2} = iv; + c_base{m,3} = boundary - steady_state(iv); + c_base{m,4} = operator; + base = false; + else + i_base{m} = 1:M.eq_nbr; + i_base{m}(etags{i,1}) = []; + c_alt{m,1} = etags{i,1}; + c_alt{m,2} = iv; + c_alt{m,3} = boundary - steady_state(iv); + c_alt{m,4} = operator; + base = true; + m = m + 1; + end + end +end +if ~base + error('OCCBIN: constraints must come by pair') +end + +function [iv,boundary,operator] = parse_constraint(str,endo_names,params,param_names) +delim = {'<=','>=','<','>'}; +[c,operator] = strsplit(str,delim); +operator = operator{1}; +iv = strmatch(strtrim(c{1}),endo_names); +% try for a number +boundary = str2num(strtrim(c{2})); +% if not a number try for a parameter name +if isempty(boundary) + k = strmatch(strtrim(c{2}),param_names); + if isempty(k) + error(['OCCBIN: illegal constraint ' str]); + end + boundary = params(k); +end diff --git a/matlab/occbin/get_residuals.m b/matlab/occbin/get_residuals.m new file mode 100644 index 000000000..f86d05cfb --- /dev/null +++ b/matlab/occbin/get_residuals.m @@ -0,0 +1,10 @@ +function r = get_residuals(ivar,lb,ub,M,oo) + +ss = oo.steady_state; +for i = 1:length(ivar) + % only one is different from zero + ss(ivar(i)) = lb(i) + ub(i); +end +oo.steady_state = ss; + +r = evaluate_model(M,oo); \ No newline at end of file diff --git a/matlab/occbin/occbin.m b/matlab/occbin/occbin.m new file mode 100644 index 000000000..ea0780d21 --- /dev/null +++ b/matlab/occbin/occbin.m @@ -0,0 +1,19 @@ +function [endo_simul,endo_simul_no_constraint,status] = occbin(M,oo,options) +% function oo=occbin(M,oo,options) solves linear models with occasionally +% binding constraints using OCCBIN by L. Guerrieri + +status = 1; +constraint_nbr = sum(strcmp(upper(M.equations_tags(:,2)),'OCCBIN'))/2; + +switch(constraint_nbr) + case 1 + [zdatalinear_ zdatapiecewise_ zdatass_ oobase_ ] = ... + solve_one_constraint(M,oo,options); + case 2 + [zdatalinear_ zdatapiecewise_ zdatass_ oobase_ ] = ... + solve_two_constraints(M,oo,options); + otherwise + error('OCCBIN can only handle two constraints in a model') +end +endo_simul = zdatapiecewise_'; +endo_simul_no_constraint = zdatalinear_'; \ No newline at end of file diff --git a/matlab/occbin/test_constraint.m b/matlab/occbin/test_constraint.m new file mode 100644 index 000000000..96e41839a --- /dev/null +++ b/matlab/occbin/test_constraint.m @@ -0,0 +1,16 @@ +function b=test_constraint(x,constr) +b = zeros(size(x,1),size(constr,1)); +for i=1:size(constr,1) + switch constr{i,4} + case '<' + b(:,i) = ~(x(:,constr{i,2}) < constr{i,3}); + case '<=' + b(:,i) = ~(x(:,constr{i,2}) <= constr{i,3}); + case '>' + b(:,i) = ~(x(:,constr{i,2}) > constr{i,3}); + case '>=' + b(:,i) = ~(x(:,constr{i,2}) >= constr{i,3}); + otherwise + error('OCCBIN: wrong inequality sign') + end +end \ No newline at end of file