Verified Commit fdf1d3fb authored by Willi Mutschler's avatar Willi Mutschler

Merge remote-tracking branch 'JohannesPfeifer/Andreasen' into moment_estimation

parents 52297f4c a7cbc637
Pipeline #3647 failed with stages
in 68 minutes and 36 seconds
function [dr] = Dynare_Unfold_Matrices(M_,options_,dr_)
% function [oo_] = Dynare_Unfold_Matrices(M_,options_,oo_)
% This function unfolds the matrices reported by Dynare up to a third order
% approximation in order to be linear in innovations
% Our notation is as follows
% z_t = f([x_t-1;u_t])
% where 'x' the state variables
% 'u' the innovations
% and v_t = [x_t-1;u_t]
% Based on code by Martin Andreasen
% INPUTS
% o M_ Matlab's structure describing the Model (initialized by dynare, see @ref{M_}).
% o options_ Matlab's structure describing the options (initialized by dynare, see @ref{options_}).
% o dr_ Matlab's decision rule structure
%
% OUTPUTS
% o dr Decision rule structure for the linear
% innovation state space
%
% SPECIAL REQUIREMENTS
% None.
% Copyright (C) 2013 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
% Dimensions
orderApp=options_.order;
nv = M_.nspred + M_.exo_nbr;
nz = M_.endo_nbr; % Number of control variables + state variables
% The cholesky factorisation of the covariance matrix
SS=M_.Sigma_e+1e-14*eye(M_.exo_nbr);
[sigma,error_mes] = chol(SS,'lower');
if error_mes ~= 0
sigma = diag(sqrt(diag(M_.Sigma_e)));
fprintf('\nThe covariance matrix of the structural shocks is not positive definite.\n');
fprintf('I continue the simulations using only the diagonal of the covariance matrix.\n');
end
% set singular covariance entries to 0
singular_entries=find(diag(M_.Sigma_e)==0);
if ~isempty(singular_entries)
sigma(singular_entries,:)=0;
sigma(:,singular_entries)=0;
end
% Steady states
f_s0 = dr_.ys(dr_.order_var); %bring in DR order
% Decision rule for the linear approximation
fv = [dr_.ghx dr_.ghu];
% Decision rule for the 2nd order.
% Given that g_2 are multiplied by 1/2, we multiply fvv by 2. See below.
if orderApp > 1
if ~isfield(dr_,'g_2')
error('GMM requires the use of k_order_solver')
end
f_2 = dr_.g_2;
fss = dr_.ghs2;
% we unfold the matrix f_22 and multiply by 2
fvv = zeros(nz,nv,nv);
index = 0;
for alfa1=1:nv
for alfa2=alfa1:nv
index = index + 1;
fvv(:,alfa1,alfa2) = f_2(:,index);
if alfa1 ~= alfa2
fvv(:,alfa2,alfa1) = f_2(:,index);
end
end
end
fvv = 2*fvv;
else
fss = zeros(nz,1);
fvv = zeros(nz,nv,nv);
end
% The third order terms are given by f_33 and by the difference between
% f_31 and f_11. We multiply by 6 because the terms in Dynare++ are already
% multiplied by 1/6 (see below)
if orderApp > 2
f_3 = dr_.g_3;
fssv = [dr_.ghxss dr_.ghuss];
% We unfold the matrix f_33 and multiply by 6
fvvv = zeros(nz,nv,nv,nv);
index = 0;
for alfa1=1:nv
for alfa2=alfa1:nv
for alfa3=alfa2:nv
index = index + 1;
fvvv(:,alfa1,alfa2,alfa3) = f_3(:,index);
% Using symmetry for alfa1 and alfa2
if alfa1 == alfa2 && alfa2 ~= alfa3 %alfa1==alfa2~=alfa3
fvvv(:,alfa1,alfa3,alfa1) = fvvv(:,alfa1,alfa2,alfa3);
fvvv(:,alfa3,alfa1,alfa1) = fvvv(:,alfa1,alfa2,alfa3);
end
% Using symmetry for alfa2 and alfa3
if alfa1 ~= alfa2 && alfa2 == alfa3 %alfa1~=alfa2==alfa3
fvvv(:,alfa2,alfa1,alfa2) = fvvv(:,alfa1,alfa2,alfa3);
fvvv(:,alfa2,alfa2,alfa1) = fvvv(:,alfa1,alfa2,alfa3);
end
% Using symmetry for alfa1,alfa2, and alfa3
if alfa1 ~= alfa2 && alfa1 ~= alfa3 && alfa2 ~= alfa3 %alfa1~=alfa2~=alfa3
fvvv(:,alfa1,alfa3,alfa2) = fvvv(:,alfa1,alfa2,alfa3);
fvvv(:,alfa3,alfa1,alfa2) = fvvv(:,alfa1,alfa2,alfa3);
fvvv(:,alfa3,alfa2,alfa1) = fvvv(:,alfa1,alfa2,alfa3);
fvvv(:,alfa2,alfa3,alfa1) = fvvv(:,alfa1,alfa2,alfa3);
fvvv(:,alfa2,alfa1,alfa3) = fvvv(:,alfa1,alfa2,alfa3);
end
end
end
end
fvvv = 6*fvvv;
fsss = zeros(nz,1); %Dynare++ solves only for Gaussian shocks
else
fsss = zeros(nz,1);
fssv = zeros(nz,nv);
fvvv = zeros(nz,nv,nv,nv);
end
dr.sig=1;
dr.f_s0=f_s0;
dr.f_v=fv;
dr.f_vv=fvv;
dr.f_ss=fss;
dr.f_vvv=fvvv;
dr.f_ssv=fssv;
dr.f_sss=fsss;
dr.sigma=sigma;
dr.order_var=dr_.order_var;
dr.inv_order_var=dr_.inv_order_var;
end
function [fval,info,exit_flag,moments_difference,modelMoments,junk1,junk2,Model,DynareOptions,GMMinfo,DynareResults]...
= GMM_Objective_Function(xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,GMMinfo,BoundsInfo,DynareResults)
% [fval,info,exit_flag,moments_difference,modelMoments,junk1,junk2,Model,DynareOptions,GMMinfo,DynareResults]...
% = GMM_Objective_Function(xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,GMMinfo,BoundsInfo,DynareResults)
% This function evaluates the objective function for GMM estimation
%
% INPUTS
% o xparam1: initial value of estimated parameters as returned by set_prior()
% o DynareDataset: data after required transformation
% o DynareOptions: Matlab's structure describing the options (initialized by dynare, see @ref{options_}).
% o Model Matlab's structure describing the Model (initialized by dynare, see @ref{M_}).
% o EstimatedParameters: Matlab's structure describing the estimated_parameters (initialized by dynare, see @ref{estim_params_}).
% o GMMInfo Matlab's structure describing the GMM settings (initialized by dynare, see @ref{bayesopt_}).
% o BoundsInfo Matlab's structure containing prior bounds
% o DynareResults Matlab's structure gathering the results (initialized by dynare, see @ref{oo_}).
%
% OUTPUTS
% o fval: value of the quadratic form of the moment difference
% o info: vector storing error code and penalty
% o exit_flag: 0 if no error, 1 of error
% o moments_difference: [numMom x 1] vector with difference of empirical and model moments
% o modelMoments: [numMom x 1] vector with model moments
% o Model: Matlab's structure describing the Model (initialized by dynare, see @ref{M_}).
% o DynareOptions: Matlab's structure describing the options (initialized by dynare, see @ref{options_}).
% o GMMinfo: Matlab's structure describing the GMM parameter options (initialized by dynare, see @ref{GMMinfo_}).
% o DynareResults: Matlab's structure gathering the results (initialized by dynare, see @ref{oo_}).
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 2013-17 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
% Initialization of the returned variables and others...
fval = NaN;
exit_flag = 1;
info = 0;
junk2 = [];
junk1 = [];
moments_difference=NaN(GMMinfo.numMom,1);
modelMoments=NaN(GMMinfo.numMom,1);
%------------------------------------------------------------------------------
% 1. Get the structural parameters & define penalties
%------------------------------------------------------------------------------
% Return, with endogenous penalty, if some parameters are smaller than the lower bound of the parameters.
if ~isequal(DynareOptions.mode_compute,1) && any(xparam1<BoundsInfo.lb)
k = find(xparam1<BoundsInfo.lb);
fval = Inf;
exit_flag = 0;
info(1) = 41;
info(4)= sum((BoundsInfo.lb(k)-xparam1(k)).^2);
return
end
% Return, with endogenous penalty, if some parameters are greater than the upper bound of the parameters.
if ~isequal(DynareOptions.mode_compute,1) && any(xparam1>BoundsInfo.ub)
k = find(xparam1>BoundsInfo.ub);
fval = Inf;
exit_flag = 0;
info(1) = 42;
info(4)= sum((xparam1(k)-BoundsInfo.ub(k)).^2);
return
end
% Set all parameters
Model = set_all_parameters(xparam1,EstimatedParameters,Model);
% Test if Q is positive definite.
if ~issquare(Model.Sigma_e) || EstimatedParameters.ncx || isfield(EstimatedParameters,'calibrated_covariances')
[Q_is_positive_definite, penalty] = ispd(Model.Sigma_e(EstimatedParameters.Sigma_e_entries_to_check_for_positive_definiteness,EstimatedParameters.Sigma_e_entries_to_check_for_positive_definiteness));
if ~Q_is_positive_definite
fval = Inf;
exit_flag = 0;
info(1) = 43;
info(4) = penalty;
return
end
if isfield(EstimatedParameters,'calibrated_covariances')
correct_flag=check_consistency_covariances(Model.Sigma_e);
if ~correct_flag
penalty = sum(Model.Sigma_e(EstimatedParameters.calibrated_covariances.position).^2);
fval = Inf;
exit_flag = 0;
info(1) = 71;
info(4) = penalty;
return
end
end
end
%------------------------------------------------------------------------------
% 2. call resol to compute steady state and model solution
%------------------------------------------------------------------------------
% Linearize the model around the deterministic steady state and extract the matrices of the state equation (T and R).
[dr_dynare_state_space,info,Model,DynareOptions,DynareResults] = resol(0,Model,DynareOptions,DynareResults);
% Return, with endogenous penalty when possible, if dynare_resolve issues an error code (defined in resol).
if info(1)
if info(1) == 3 || info(1) == 4 || info(1) == 5 || info(1)==6 ||info(1) == 19 ||...
info(1) == 20 || info(1) == 21 || info(1) == 23 || info(1) == 26 || ...
info(1) == 81 || info(1) == 84 || info(1) == 85 || info(1) == 86
%meaningful second entry of output that can be used
fval = Inf;
info(4) = info(2);
exit_flag = 0;
return
else
fval = Inf;
info(4) = 0.1;
exit_flag = 0;
return
end
end
% % check endogenous prior restrictions
% info=endogenous_prior_restrictions(T,R,Model,DynareOptions,DynareResults);
% if info(1)
% fval = Inf;
% info(4)=info(2);
% exit_flag = 0;
% return
% end
%------------------------------------------------------------------------------
% 3. Set up state-space with linear innovations
%------------------------------------------------------------------------------
% Transformation of the approximated solution
DynareResults.gmm.dr = Dynare_Unfold_Matrices(Model,DynareOptions,dr_dynare_state_space);
% We set up the alternative state space representation and use only selected endogenous variables
DynareResults.gmm.dr = State_Space_LinearInov(Model,DynareResults.gmm.dr,GMMinfo.control_indices,GMMinfo.state_indices);
%------------------------------------------------------------------------------
% 4. Compute Moments of the model solution for normal innovations
%------------------------------------------------------------------------------
DynareResults.gmm.unconditionalmoments = Get_Pruned_Unconditional_Moments(DynareResults.gmm.dr,DynareOptions,GMMinfo,DynareOptions.gmm.autolag);
% Get the moments implied by the model solution that are matched
if DynareOptions.gmm.centeredmoments
modelMoments = collect_Moments(DynareResults.gmm.unconditionalmoments.E_y,DynareResults.gmm.unconditionalmoments.Var_y,DynareResults.gmm.unconditionalmoments.autoCov_y,DynareOptions);
else
modelMoments = collect_Moments(DynareResults.gmm.unconditionalmoments.E_y,DynareResults.gmm.unconditionalmoments.E_yy,DynareResults.gmm.unconditionalmoments.autoE_yy,DynareOptions);
end
%------------------------------------------------------------------------------
% 4. Compute quadratic target function using weighting matrix W
%------------------------------------------------------------------------------
moments_difference = DynareResults.gmm.datamoments.momentstomatch-modelMoments;
fval = moments_difference'*DynareResults.gmm.W*moments_difference;
if DynareOptions.gmm.penalized_estimator
fval=fval+(xparam1-GMMinfo.p1)'/diag(GMMinfo.p2)*(xparam1-GMMinfo.p1);
end
end
function [GIRF] = Generalized_IRF(dr,Model,DynareOptions,shocksize,control_indices,vectorMom3,xf,xs)
% [GIRF] = Generalized_IRF(dr,Model,DynareOptions,shocksize,control_indices,vectorMom3,xf,xs)
% This function sets up the state space linear in innovations and calls the function computing
% the closed-form expressions for the impulse response functions using the pruning method
% when using the following definition of the impulse response function.
%
% INPUTS
% o dr_ Matlab's decision rule structure
% o Model Matlab's structure describing the Model (initialized by dynare, see @ref{M_}).
% o DynareOptions Matlab's structure describing the options (initialized by dynare, see @ref{options_}).
% o shocksize [double] (p*1) vector of shocks (measured in standard deviations).
% o control_indices [double] (ny*1) vector of indices of requested
% endogenous variables in decision rules (i.e. M_.endo_names(dr.inv_order_var,:))
% o vectorMom3 [double] (p*1) vector of third moments
% o xf [double] (p*1) vector of first order state at which
% to compute IRF
% o xs [double] (p*1) vector of second order state at which
% to compute IRF
%
% OUTPUT:
% o GIRF [structure] strucutre containing the GIRFs, including
% GIRF.y impulse response function for y = g(x,sig) (sum of all parts)
%
% GIRF.parts.yf impulse responses at first order
% GIRF.parts.ys impulse responses at second order
% GIRF.parts.yrd impulse responses at thrid order
%
% GIRF.x impulse response function for x = h(x,sig) (sum of all parts)
%
% GIRF.parts.xf impulse responses at first order
% GIRF.parts.xs impulse responses at second order
% GIRF.parts.xrd impulse responses at third order
%
% Copyright (C) 2013-17 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
state_indices=[Model.nstatic+1:Model.nstatic+Model.nspred];
if nargin<5
control_indices=dr.inv_order_var(1:Model.orig_endo_nbr); % variables in matrices are in order_var ordering and need to be mapped to declaration order using inv_order_var
end
if nargin<6
vectorMom3=zeros(1,length(shocksize));
end
if nargin<7
xf=[];
end
if nargin<8
xs=[];
end
if DynareOptions.order>1 && ~DynareOptions.pruning
error('GIRFs require the option pruning')
end
%------------------------------------------------------------------------------
% 1. Make sure state-space in linear innovations is set up, including
% updating when control variables have changed
%------------------------------------------------------------------------------
% Transformation of the approximated solution
dr = Dynare_Unfold_Matrices(Model,DynareOptions,dr);
% We set up the alternative state space representation and use only selected endogenous variables
dr = State_Space_LinearInov(Model,dr,control_indices,state_indices);
%------------------------------------------------------------------------------
% 2. Compute Generalized IRFs for state-space at parameters
%------------------------------------------------------------------------------
[GIRF] = Generalized_IRF_core(dr,DynareOptions,shocksize,vectorMom3,xf,xs);
\ No newline at end of file
function [GIRF] = Generalized_IRF_core(dr,options_,shockSize,vectorMom3,xf,xs)
% [GIRF] = Generalized_IRF_core(dr,options_,shockSize,vectorMom3,xf,xs)
% This function computes closed-form expressions for the impulse response
% functions using the pruning method when using the following definition
% of the impulse response function.
%
% impulse(x_t+l) = E_t[w_t|eps_t+1=eps_t+1+shockSize, eps_t+2, eps_t+3, ..., eps_t+l]
% -E_t[w_t|eps_t+1=eps_t+1 , eps_t+2, eps_t+3, ..., eps_t+l]
%
% Note that this definition is the one behind the
% so-called Generalized Impulse Reponse Function (GIRF). The definition is
% based on Andreasen (2012), European Economic Review
%
% INPUTS
% o dr_ Matlab's decision rule structure
% o options_ Matlab's structure describing the options (initialized by dynare, see @ref{options_}).
% o shocksize [double] (p*1) vector of shocks (measured in standard deviations).
% o vectorMom3 [double] (p*1) vector of third moments
% o xf [double] (p*1) vector of first order state at which
% to compute IRF
% o xs [double] (p*1) vector of second order state at which
% to compute IRF
%
% OUTPUT:
% o GIRF [structure] strucutre containing the GIRFs, including
% GIRF.y impulse response function for y = g(x,sig) (sum of all parts)
%
% GIRF.parts.yf impulse responses at first order
% GIRF.parts.ys impulse responses at second order
% GIRF.parts.yrd impulse responses at thrid order
%
% GIRF.x impulse response function for x = h(x,sig) (sum of all parts)
%
% GIRF.parts.xf impulse responses at first order
% GIRF.parts.xs impulse responses at second order
% GIRF.parts.xrd impulse responses at third order
%
% NOTE: if xf = [] and xs = [], then we use the default setting and apply
% the mean values of xf and xs.
%
% Based on code by Martin M. Andreasen, 2012
%
% Copyright (C) 2013-17 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
orderApp=options_.order;
IRFlength=options_.irf;
% read out decision rule coefficients
eta = dr.eta;
sig = dr.sig;
if ~isfield(dr,'h_v')
error('State Space with linear innovations has not been set up')
end
hx = dr.h_v;
gx = dr.g_v;
% The dimension of the model
[ny,nx] = size(gx);
ne = size(eta,2);
GIRFy = zeros(ny,IRFlength,3);
GIRFx = zeros(nx,IRFlength,3);
% read out decision rule coefficients
if orderApp>1
hxx = dr.h_vv;
gxx = dr.g_vv;
hss = dr.h_ss;
% Defining matrices
HHxxtil = 1/2*reshape(hxx,nx,nx^2);
GGxxtil = 1/2*reshape(gxx,ny,nx^2);
end
if orderApp>2
hxxx = dr.h_vvv;
gxxx = dr.g_vvv;
hssx = dr.h_ssv;
gssx = dr.g_ssv;
% Defining matrices
HHxxxtil = 1/6*reshape(hxxx,nx,nx^3);
GGxxxtil = 1/6*reshape(gxxx,ny,nx^3);
end
% Checking on shockSize
[n1,n2] = size(shockSize);
% ensuring that we have a column vector
if n2 > n1
shockSize = shockSize';
end
if size(shockSize,1) ~= ne
error('lenght of shockSize must be equal to ne');
end
%% Computing the mean values for xf and xs
sigeta = sig*eta;
eta_eta = eta*eta';
AA = (eye(nx)-hx)\eye(nx); %inv(eye(nx)-hx);
vec_var_xf = (eye(nx^2)-kron(hx,hx))\reshape(sig^2*eta_eta,nx^2,1);
var_xf = reshape(vec_var_xf,nx,nx);
mean_xfxf = var_xf;
if orderApp>1
mean_xs = AA*(0.5*reshape(hxx,nx,nx^2,1)*reshape(mean_xfxf,nx^2,1)+0.5*hss*sig^2);
end
if isempty(xf)
xf = zeros(nx,1);
end
if isempty(xs) & orderApp>1
xs = mean_xs;
end
%% Impulse response functions at first order
GIRFy1st = zeros(ny,IRFlength);
GIRFxf = zeros(nx,IRFlength);
for k=1:IRFlength
GIRFxf(:,k) = hx^(k-1)*sigeta*shockSize;
GIRFy1st(:,k) = gx*GIRFxf(:,k);
end
% Saving output
GIRFy(:,:,1) = GIRFy1st;
GIRFx(:,:,1) = GIRFxf;
GIRF.parts.xf=GIRFxf;
GIRF.y=GIRFy;
GIRF.x=GIRFx;
GIRF.parts.yf=GIRFy1st;
if orderApp == 1
GIRF.y=squeeze(GIRFy(:,:,1));
GIRF.x=squeeze(GIRFx(:,:,1));
return
end
%% Impulse response functions at second order
GIRFy2nd = zeros(ny,IRFlength);
GIRFxs = zeros(nx,IRFlength);
GIRFxfxf = zeros(nx^2,IRFlength);
consQ = kron(sigeta*shockSize,sigeta*shockSize);
% Computing GIRFxfxf
h= dyn_waitbar(0,'Computing GIRF for x1st^2');
for k=1:IRFlength
GIRFxfxf(:,k) = kron(hx^k*xf,hx^(k-1)*sigeta*shockSize) + kron(hx^(k-1)*sigeta*shockSize,hx^k*xf) + ...
kron(hx^(k-1),hx^(k-1))*consQ;
dyn_waitbar(k/IRFlength,h);
end
dyn_waitbar_close(h);
% Computing GIRFxs and GIRFy2nd
h = dyn_waitbar(0/IRFlength,'Computing GIRF for x2nd');
for k=1:IRFlength
if k > 1
GIRFxs(:,k) = hx*GIRFxs(:,k-1) + HHxxtil*GIRFxfxf(:,k-1);
end
GIRFy2nd(:,k) = gx*(GIRFxf(:,k)+GIRFxs(:,k)) + GGxxtil*GIRFxfxf(:,k);
dyn_waitbar(k/IRFlength,h);
end
dyn_waitbar_close(h);
% Saving output
GIRFy(:,:,2) = GIRFy2nd;
GIRFx(:,:,2) = GIRFxf+GIRFxs;
GIRF.parts.xs=GIRFxs;
GIRF.parts.ys=GIRFy2nd;
if orderApp == 2
GIRF.y=squeeze(GIRFy(:,:,2));
GIRF.x=squeeze(GIRFx(:,:,2));
return
end
%% Impulse response functions a third order
GIRFy3rd = zeros(ny,IRFlength);
GIRFxrd = zeros(nx,IRFlength);
GIRFxfxs = zeros(nx^2,IRFlength);
GIRFxfxfxf = zeros(nx^3,IRFlength);
xf_xf = kron(xf,xf);
etaEps_etaEps = zeros(nx^2,1);
idx = 0;
for gama2=1:nx
for gama1=1:nx
idx = idx + 1;
etaEps_etaEps(idx,1) = eta(gama2,:)*eta(gama1,:)'*sig^2;
end
end
% Computing GIRFxfxfxf
h = dyn_waitbar(0/IRFlength,'Computing GIRF for x1st^3');
for k=1:IRFlength
term1 = zeros(nx^3,1);
term2 = zeros(nx^3,1);
term3 = zeros(nx^3,1);
hxlm1 = sig*hx^(k-1);
for j=1:k
omegaj = zeros(nx^3,1);
hxlmj = sig*hx^(k-j);
idx = 0;
for gama3=1:nx
for gama2=1:nx
for gama1=1:nx
idx = idx + 1;
%for phi1=1:ne
% omegaj(idx,1) = omegaj(idx,1) + hxlmj(gama3,:)*sigeta(:,phi1)...
% *hxlm1(gama2,:)*sigeta*shockSize...
% *hxlmj(gama1,:)*sigeta(:,phi1);
%end
omegaj(idx,1) = sig^2*hxlmj(gama3,:)*eta_eta*hxlmj(gama1,:)'...
*hxlm1(gama2,:)*sigeta*shockSize;
end
end
end
term1 = term1 + omegaj(:,1);
term2 = term2 + kron(kron(hx^(k-j),hx^(k-j))*etaEps_etaEps,hx^(k-1)*sigeta*shockSize);
term3 = term3 + kron(hx^(k-1)*sigeta*shockSize,kron(hx^(k-j),hx^(k-j))*etaEps_etaEps);
end
GIRFxfxfxf(:,k) = kron(hx^k*xf,kron((hx^k*xf+hx^(k-1)*sigeta*shockSize),hx^(k-1)*sigeta*shockSize)) ...
+ kron(hx^k*xf,kron(hx^(k-1)*sigeta*shockSize,hx^k*xf)) ...