1function [out,info] = get_perturbation_params_derivs_numerical_objective(params, outputflag, estim_params, M, oo, options)
2%function [out,info] = get_perturbation_params_derivs_numerical_objective(params, outputflag, estim_params, M, oo, options)
3% -------------------------------------------------------------------------
4% Objective function to compute numerically the Jacobians used for get_perturbation_params_derivs
5% =========================================================================
6% INPUTS
7%   params:         [vector] parameter values at which to evaluate objective function
8%                   stderr parameters come first, corr parameters second, model parameters third
9%   outputflag:     [string] flag which objective to compute (see below)
10%   estim_params:   [structure] storing the estimation information
11%   M:              [structure] storing the model information
12%   oo:             [structure] storing the solution results
13%   options:        [structure] storing the options
14% -------------------------------------------------------------------------
15%
16% OUTPUT
17%   out (dependent on outputflag and order of approximation):
18%     - 'perturbation_solution':  out = out1 = [vec(Sigma_e);vec(ghx);vec(ghu)]; (order==1)
19%                                 out = out2 = [out1;vec(ghxx);vec(ghxu);vec(ghuu);vec(ghs2)]; (order==2)
20%                                 out = out3 = [out1;out2;vec(ghxxx);vec(ghxxu);vec(ghxuu);vec(ghuuu);vec(ghxss);vec(ghuss)]; (order==3)
21%     - 'dynamic_model':          out = [Yss; vec(g1); vec(g2); vec(g3)]
22%     - 'Kalman_Transition':      out = [Yss; vec(KalmanA); dyn_vech(KalmanB*Sigma_e*KalmanB')];
23%     all in DR-order
24%   info            [integer] output from resol
25% -------------------------------------------------------------------------
26% This function is called by
27%   * get_perturbation_params_derivs.m (previously getH.m)
28% -------------------------------------------------------------------------
29% This function calls
30%   * [M.fname,'.dynamic']
31%   * resol
32%   * dyn_vech
33% =========================================================================
34% Copyright (C) 2019-2020 Dynare Team
35%
36% This file is part of Dynare.
37%
38% Dynare is free software: you can redistribute it and/or modify
39% it under the terms of the GNU General Public License as published by
40% the Free Software Foundation, either version 3 of the License, or
41% (at your option) any later version.
42%
43% Dynare is distributed in the hope that it will be useful,
44% but WITHOUT ANY WARRANTY; without even the implied warranty of
45% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
46% GNU General Public License for more details.
47%
48% You should have received a copy of the GNU General Public License
49% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
50% =========================================================================
51
52%% Update stderr, corr and model parameters and compute perturbation approximation and steady state with updated parameters
53M = set_all_parameters(params,estim_params,M);
54[~,info,M,options,oo] = resol(0,M,options,oo);
55Sigma_e = M.Sigma_e;
56
57if info(1) > 0
58    % there are errors in the solution algorithm
59    out = [];
60    return
61else
62    ys = oo.dr.ys; %steady state of model variables in declaration order
63    ghx = oo.dr.ghx; ghu = oo.dr.ghu;
64    if options.order > 1
65        ghxx = oo.dr.ghxx; ghxu = oo.dr.ghxu; ghuu = oo.dr.ghuu; ghs2 = oo.dr.ghs2;
66    end
67    if options.order > 2
68        ghxxx = oo.dr.ghxxx; ghxxu = oo.dr.ghxxu; ghxuu = oo.dr.ghxuu; ghxss = oo.dr.ghxss; ghuuu = oo.dr.ghuuu; ghuss = oo.dr.ghuss;
69    end
70end
71Yss = ys(oo.dr.order_var); %steady state of model variables in DR order
72
73%% out = [vec(Sigma_e);vec(ghx);vec(ghu);vec(ghxx);vec(ghxu);vec(ghuu);vec(ghs2);vec(ghxxx);vec(ghxxu);vec(ghxuu);vec(ghuuu);vec(ghxss);vec(ghuss)]
74if strcmp(outputflag,'perturbation_solution')
75    out = [Sigma_e(:); ghx(:); ghu(:)];
76    if options.order > 1
77        out = [out; ghxx(:); ghxu(:); ghuu(:); ghs2(:);];
78    end
79    if options.order > 2
80        out = [out; ghxxx(:); ghxxu(:); ghxuu(:); ghuuu(:); ghxss(:); ghuss(:)];
81    end
82end
83
84%% out = [Yss; vec(g1); vec(g2); vec(g3)]; of all endogenous variables, in DR order
85if strcmp(outputflag,'dynamic_model')
86    [I,~] = find(M.lead_lag_incidence'); %I is used to evaluate dynamic model files
87    if options.order == 1
88        [~, g1] = feval([M.fname,'.dynamic'], ys(I), oo.exo_steady_state', M.params, ys, 1);
89        out = [Yss; g1(:)];
90    elseif options.order == 2
91        [~, g1, g2] = feval([M.fname,'.dynamic'], ys(I), oo.exo_steady_state', M.params, ys, 1);
92        out = [Yss; g1(:); g2(:)];
93    elseif options.order == 3
94        [~, g1, g2, g3] = feval([M.fname,'.dynamic'], ys(I), oo.exo_steady_state', M.params, ys, 1);
95        g3 = unfold_g3(g3, length(ys(I))+M.exo_nbr);
96        out = [Yss; g1(:); g2(:); g3(:)];
97    end
98end
99
100%% out = [Yss; vec(KalmanA); dyn_vech(KalmanB*Sigma_e*KalmanB')]; in DR order, where A and B are Kalman transition matrices
101if strcmp(outputflag,'Kalman_Transition')
102    if options.order == 1
103        KalmanA = zeros(M.endo_nbr,M.endo_nbr);
104        KalmanA(:,M.nstatic+(1:M.nspred)) = ghx;
105        Om = ghu*Sigma_e*transpose(ghu);
106        out = [Yss; KalmanA(:); dyn_vech(Om)];
107    else
108        error('''get_perturbation_params_derivs_numerical_objective.m'': Kalman_Transition works only at order=1');
109    end
110end
111