1function out = identification_numerical_objective(params, outputflag, estim_params, M, oo, options, indpmodel, indpstderr, indpcorr, indvar, useautocorr, nlags, grid_nbr)
2%function out = identification_numerical_objective(params, outputflag, estim_params, M, oo, options, indpmodel, indpstderr, indpcorr, indvar, useautocorr, nlags, grid_nbr)
3% -------------------------------------------------------------------------
4% Objective function to compute numerically the Jacobians used for identification analysis
5% Previously this function was called thet2tau.m
6% =========================================================================
7% INPUTS
8%   params:         [vector] parameter values at which to evaluate objective function
9%                   stderr parameters come first, corr parameters second, model parameters third
10%   outputflag:     [integer] flag which objective to compute (see below)
11%   estim_params:   [structure] storing the estimation information
12%   M:              [structure] storing the model information
13%   oo:             [structure] storing the reduced form solution results
14%   options:        [structure] storing the options
15%   indpmodel:      [vector] Index of model parameters
16%   indpstderr:     [vector] Index of stderr parameters
17%   indpcorr:       [matrix] Index of corr parameters
18%   indvar:         [vector] Index of selected or observed variables
19% -------------------------------------------------------------------------
20% OUTPUTS
21%   out:    dependent on outputflag
22%   *  0:   out = [Yss; vec(A); vec(B); dyn_vech(Sig_e)]; of indvar variables only, in DR order. This is needed to compute dTAU and Komunjer and Ng's D.
23%           Note that Jacobian of Om is computed in get_identification_Jacobians.m (previously getJJ.m) or get_first_order_solution_params_deriv.m (previously getH.m) from Jacobian of B and Sigma_e, because this is more efficient due to some testing with analytical derivatives from An and Schorfheide model
24%   *  1:   out = [vech(cov(Y_t,Y_t)); vec(cov(Y_t,Y_{t-1}); ...; vec(cov(Y_t,Y_{t-nlags})] of indvar variables, in DR order. This is needed to compute Iskrev's J.
25%   *  2:   out = vec(spectral density) with dimension [var_nbr^2*grid_nbr,1] Spectral density of indvar variables evaluated at (grid_nbr/2+1) discretized points in the interval [0;pi]. This is needed for Qu and Tkachenko's G.
26%   * -1:   out = g1(:); of all variables, in DR order. This is needed to compute dLRE.
27%   * -2:   out = [Yss; vec(A); dyn_vech(B*Sigma_e*B')]; of indvar variables only, in DR order. This is used to compute numerically second derivatives d2A, d2Om d2Yss in get_first_order_solution_params_deriv.m (previously getH.m) for kronflag=1
28% where Yss is steady in DR order, A and B solution matrices of Kalman
29% transition equation, Sig_e the covariance of exogenous shocks, g1 the
30% Jacobian of the dynamic model equations, and Y_t selected variables
31% -------------------------------------------------------------------------
32% This function is called by
33%   * get_identification_jacobians.m (previously getJJ.m)
34% -------------------------------------------------------------------------
35% This function calls
36%   * [M.fname,'.dynamic']
37%   * dyn_vech
38%   * resol
39%   * vec
40% =========================================================================
41% Copyright (C) 2011-2020 Dynare Team
42%
43% This file is part of Dynare.
44%
45% Dynare is free software: you can redistribute it and/or modify
46% it under the terms of the GNU General Public License as published by
47% the Free Software Foundation, either version 3 of the License, or
48% (at your option) any later version.
49%
50% Dynare is distributed in the hope that it will be useful,
51% but WITHOUT ANY WARRANTY; without even the implied warranty of
52% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
53% GNU General Public License for more details.
54%
55% You should have received a copy of the GNU General Public License
56% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
57% =========================================================================
58
59%% Update stderr, corr and model parameters
60%note that if no estimated_params_block is given, then all stderr and model parameters are selected but no corr parameters
61if length(params) > length(indpmodel)
62    if isempty(indpstderr)==0 && isempty(estim_params.var_exo) %if there are stderr parameters but no estimated_params_block
63        %provide temporary necessary information for stderr parameters
64        estim_params.nvx = length(indpstderr);
65        estim_params.var_exo = indpstderr';
66    end
67    if isempty(indpmodel)==0 && isempty(estim_params.param_vals) %if there are model parameters but no estimated_params_block
68        %provide temporary necessary information for model parameters
69        estim_params.np = length(indpmodel);
70        estim_params.param_vals = indpmodel';
71    end
72    M = set_all_parameters(params,estim_params,M); %this function can only be used if there is some information in estim_params
73else
74    %if there are only model parameters, we don't need to use set_all_parameters
75    M.params(indpmodel) = params;
76end
77
78%% compute Kalman transition matrices and steady state with updated parameters
79[~,info,M,options,oo] = resol(0,M,options,oo);
80options = rmfield(options,'options_ident');
81pruned = pruned_state_space_system(M, options, oo.dr, indvar, nlags, useautocorr, 0);
82
83%% out = [vech(cov(Y_t,Y_t)); vec(cov(Y_t,Y_{t-1}); ...; vec(cov(Y_t,Y_{t-nlags})] of indvar variables, in DR order. This is Iskrev (2010)'s J matrix.
84if outputflag == 1
85    out = dyn_vech(pruned.Var_y);
86    for i = 1:nlags
87        if useautocorr
88            out = [out;vec(pruned.Corr_yi(:,:,i))];
89        else
90            out = [out;vec(pruned.Var_yi(:,:,i))];
91        end
92    end
93end
94
95%% out = vec(g_omega). This is needed for Qu and Tkachenko (2012)'s G matrix.
96if outputflag == 2
97    % This computes the spectral density g_omega where the interval [-pi;\pi] is discretized by grid_nbr points
98    freqs = (0 : pi/(grid_nbr/2):pi);% we focus only on positive values including the 0 frequency
99    tpos  = exp( sqrt(-1)*freqs); %Fourier frequencies
100    IA = eye(size(pruned.A,1));
101    var_nbr = size(pruned.C,1);
102    out = zeros(var_nbr^2*length(freqs),1);
103    kk = 0;
104    for ig = 1:length(freqs)
105        Transferfct = pruned.D + pruned.C*((tpos(ig)*IA-pruned.A)\pruned.B);
106        g_omega = (1/(2*pi))*(Transferfct*pruned.Varinov*Transferfct'); % note that ' is the conjugate transpose
107        kk = kk+1;
108        out(1 + (kk-1)*var_nbr^2 : kk*var_nbr^2) = g_omega(:);
109    end
110end