1function [steady_state,params,check] = dyn_ramsey_static(ys_init,M,options_,oo) 2 3% function [steady_state,params,check] = dyn_ramsey_static_(ys_init,M,options_,oo) 4% Computes the steady state for optimal policy 5% 6% INPUTS 7% ys_init: vector of endogenous variables or instruments 8% M: Dynare model structure 9% options: Dynare options structure 10% oo: Dynare results structure 11% 12% OUTPUTS 13% steady_state: steady state value 14% params: parameters at steady state, potentially updated by 15% steady_state file 16% check: error indicator, 0 if everything is OK 17% 18% SPECIAL REQUIREMENTS 19% none 20 21% Copyright (C) 2003-2018 Dynare Team 22% 23% This file is part of Dynare. 24% 25% Dynare is free software: you can redistribute it and/or modify 26% it under the terms of the GNU General Public License as published by 27% the Free Software Foundation, either version 3 of the License, or 28% (at your option) any later version. 29% 30% Dynare is distributed in the hope that it will be useful, 31% but WITHOUT ANY WARRANTY; without even the implied warranty of 32% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 33% GNU General Public License for more details. 34% 35% You should have received a copy of the GNU General Public License 36% along with Dynare. If not, see <http://www.gnu.org/licenses/>. 37 38 39params = M.params; 40check = 0; 41options_.steadystate.nocheck = 1; %locally disable checking because Lagrange multipliers are not accounted for in evaluate_steady_state_file 42 % dyn_ramsey_static_1 is a subfunction 43nl_func = @(x) dyn_ramsey_static_1(x,M,options_,oo); 44 45% check_static_model is a subfunction 46if check_static_model(ys_init,M,options_,oo) && ~options_.steadystate_flag 47 steady_state = ys_init; 48 return 49elseif options_.steadystate_flag 50 k_inst = []; 51 inst_nbr = size(options_.instruments,1); 52 for i = 1:inst_nbr 53 k_inst = [k_inst; strmatch(options_.instruments{i}, M.endo_names, 'exact')]; 54 end 55 if inst_nbr == 1 56 %solve for instrument, using univariate solver, starting at initial value for instrument 57 [inst_val, info1]= csolve(nl_func,ys_init(k_inst),'',options_.solve_tolf,options_.ramsey.maxit); 58 if info1==1 || info1==3 59 check=81; 60 end 61 if info1==4 62 check=87; 63 end 64 else 65 %solve for instrument, using multivariate solver, starting at 66 %initial value for instrument 67 opt = options_; 68 opt.jacobian_flag = false; 69 [inst_val,info1] = dynare_solve(nl_func,ys_init(k_inst), ... 70 opt); 71 if info1~=0 72 check=81; 73 end 74 end 75 ys_init(k_inst) = inst_val; 76 exo_ss = [oo.exo_steady_state oo.exo_det_steady_state]; 77 [xx,params] = evaluate_steady_state_file(ys_init,exo_ss,M,options_,~options_.steadystate.nocheck); %run steady state file again to update parameters 78 [~,~,steady_state] = nl_func(inst_val); %compute and return steady state 79else 80 n_var = M.orig_endo_nbr; 81 xx = oo.steady_state(1:n_var); 82 opt = options_; 83 opt.jacobian_flag = false; 84 [xx,info1] = dynare_solve(nl_func,xx,opt); 85 if info1~=0 86 check=81; 87 end 88 [~,~,steady_state] = nl_func(xx); 89end 90 91 92 93function [resids,rJ,steady_state] = dyn_ramsey_static_1(x,M,options_,oo) 94resids = []; 95rJ = []; 96mult = []; 97 98% recovering usefull fields 99params = M.params; 100endo_nbr = M.endo_nbr; 101endo_names = M.endo_names; 102orig_endo_nbr = M.orig_endo_nbr; 103aux_vars_type = [M.aux_vars.type]; 104orig_endo_aux_nbr = orig_endo_nbr + min(find(aux_vars_type == 6)) - 1; 105orig_eq_nbr = M.orig_eq_nbr; 106inst_nbr = orig_endo_aux_nbr - orig_eq_nbr; 107% indices of Lagrange multipliers 108fname = M.fname; 109 110 111if options_.steadystate_flag 112 k_inst = []; 113 instruments = options_.instruments; 114 for i = 1:size(instruments,1) 115 k_inst = [k_inst; strmatch(instruments{i}, endo_names, 'exact')]; 116 end 117 ys_init=zeros(size(oo.steady_state)); %create starting vector for steady state computation as only instrument value is handed over 118 ys_init(k_inst) = x; %set instrument, the only value required for steady state computation, to current value 119 [x,params,check] = evaluate_steady_state_file(ys_init,... %returned x now has size endo_nbr as opposed to input size of n_instruments 120 [oo.exo_steady_state; ... 121 oo.exo_det_steady_state], ... 122 M,options_,~options_.steadystate.nocheck); 123 if any(imag(x(1:M.orig_endo_nbr))) %return with penalty 124 resids=ones(inst_nbr,1)+sum(abs(imag(x(1:M.orig_endo_nbr)))); %return with penalty 125 steady_state=NaN(endo_nbr,1); 126 return 127 end 128 if check %return 129 resids=ones(inst_nbr,1)+sum(abs(x(1:M.orig_endo_nbr))); %return with penalty 130 steady_state=NaN(endo_nbr,1); 131 return 132 end 133 134end 135 136xx = zeros(endo_nbr,1); %initialize steady state vector 137xx(1:M.orig_endo_nbr) = x(1:M.orig_endo_nbr); %set values of original endogenous variables based on steady state file or initial value 138 139% setting steady state of auxiliary variables that depends on original endogenous variables 140if any([M.aux_vars.type] ~= 6) %auxiliary variables other than multipliers 141 needs_set_auxiliary_variables = 1; 142 if M.set_auxiliary_variables 143 fh = str2func([M.fname '.set_auxiliary_variables']); 144 s_a_v_func = @(z) fh(z,... 145 [oo.exo_steady_state,... 146 oo.exo_det_steady_state],... 147 params); 148 else 149 s_a_v_func = z; 150 end 151 xx = s_a_v_func(xx); 152else 153 needs_set_auxiliary_variables = 0; 154end 155 156% value and Jacobian of objective function 157ex = zeros(1,M.exo_nbr); 158[U,Uy,Uyy] = feval([fname '.objective.static'],x,ex, params); 159Uyy = reshape(Uyy,endo_nbr,endo_nbr); 160 161% set multipliers and auxiliary variables that 162% depends on multipliers to 0 to compute residuals 163if (options_.bytecode) 164 [chck, res, junk] = bytecode('static',xx,[oo.exo_steady_state oo.exo_det_steady_state], ... 165 params, 'evaluate'); 166 fJ = junk.g1; 167else 168 [res,fJ] = feval([fname '.static'],xx,[oo.exo_steady_state oo.exo_det_steady_state], ... 169 params); 170end 171% index of multipliers and corresponding equations 172% the auxiliary variables before the Lagrange multipliers are treated 173% as ordinary endogenous variables 174aux_eq = [1:orig_endo_aux_nbr, orig_endo_aux_nbr+orig_eq_nbr+1:size(fJ,1)]; 175A = fJ(1:orig_endo_aux_nbr,orig_endo_nbr+find(aux_vars_type==6)); 176y = res(1:orig_endo_aux_nbr); 177mult = -A\y; 178 179resids1 = y+A*mult; 180if inst_nbr == 1 181 r1 = sqrt(resids1'*resids1); 182else 183 [q,r,e] = qr([A y]'); 184 k = size(A,1)+(1-inst_nbr:0); 185 r1 = r(end,k)'; 186end 187if options_.steadystate_flag 188 resids = r1; 189else 190 resids = [res(orig_endo_nbr+(1:orig_endo_nbr-inst_nbr)); r1]; 191end 192rJ = []; 193if needs_set_auxiliary_variables 194 steady_state = s_a_v_func([xx(1:orig_endo_aux_nbr); mult]); 195else 196 steady_state = [xx(1:orig_endo_aux_nbr); mult]; 197end 198 199function result = check_static_model(ys,M,options_,oo) 200result = false; 201if (options_.bytecode) 202 [chck, res, ~] = bytecode('static',ys,[oo.exo_steady_state oo.exo_det_steady_state], ... 203 M.params, 'evaluate'); 204else 205 res = feval([M.fname '.static'],ys,[oo.exo_steady_state oo.exo_det_steady_state], ... 206 M.params); 207end 208if norm(res) < options_.solve_tolf 209 result = true; 210end 211