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