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