diff --git a/matlab/conditional_variance_decomposition.m b/matlab/conditional_variance_decomposition.m
index 4c79249efa18fe078770e888a4652a4ac799a310..088a62ddc8c63eb90063afd2efad09fa52b87221 100644
--- a/matlab/conditional_variance_decomposition.m
+++ b/matlab/conditional_variance_decomposition.m
@@ -1,11 +1,13 @@
-function [ConditionalVarianceDecomposition, ConditionalVarianceDecomposition_ME]= conditional_variance_decomposition(StateSpaceModel, Steps, SubsetOfVariables,sigma_e_is_diagonal)
+function [ConditionalVarianceDecomposition, ConditionalVarianceDecomposition_ME]= conditional_variance_decomposition(M_,options_,dr, Steps, SubsetOfVariables)
% This function computes the conditional variance decomposition of a given state space model
% for a subset of endogenous variables.
%
% INPUTS
-% StateSpaceModel [structure] Specification of the state space model.
-% Steps [integer] 1*h vector of dates.
-% SubsetOfVariables [integer] 1*q vector of indices (declaration order).
+% M_ [struct] Matlab's structure describing the Model (initialized by dynare, see @ref{M_}).
+% options_ [struct] Matlab's structure describing the options (initialized by dynare, see @ref{options_}).
+% dr : [struct] Dynare decision rules structure
+% Steps [integer] 1*h vector of dates.
+% SubsetOfVariables [integer] 1*q vector of indices (declaration order).
%
% OUTPUTS
% ConditionalVarianceDecomposition [double] [n h p] array, where
@@ -17,7 +19,7 @@ function [ConditionalVarianceDecomposition, ConditionalVarianceDecomposition_ME]
% h is the number of Steps
% p is the number of state innovations and
-% Copyright © 2010-2021 Dynare Team
+% Copyright © 2010-2023 Dynare Team
%
% This file is part of Dynare.
%
@@ -39,22 +41,19 @@ if any(Steps <= 0)
'positive'])
end
-number_of_state_innovations = ...
- StateSpaceModel.number_of_state_innovations;
-transition_matrix = StateSpaceModel.transition_matrix;
-number_of_state_equations = ...
- StateSpaceModel.number_of_state_equations;
-order_var = StateSpaceModel.order_var;
+[transition_matrix, impulse_matrix] = kalman_transition_matrix(dr,(1:M_.endo_nbr)',M_.nstatic+(1:M_.nspred)',M_.exo_nbr);
+
+number_of_state_innovations = M_.exo_nbr;
+number_of_state_equations = M_.endo_nbr;
+order_var = dr.order_var;
nSteps = length(Steps);
ConditionalVariance = zeros(number_of_state_equations,nSteps,number_of_state_innovations);
-if StateSpaceModel.sigma_e_is_diagonal
- B = StateSpaceModel.impulse_matrix.* ...
- repmat(sqrt(diag(StateSpaceModel.state_innovations_covariance_matrix)'),...
- number_of_state_equations,1);
+if M_.sigma_e_is_diagonal
+ B = impulse_matrix.* repmat(sqrt(diag(M_.Sigma_e)'),number_of_state_equations,1);
else
- B = StateSpaceModel.impulse_matrix*chol(StateSpaceModel.state_innovations_covariance_matrix)';
+ B = impulse_matrix*chol(M_.Sigma_e)';
end
for i=1:number_of_state_innovations
@@ -88,9 +87,9 @@ end
% get intersection of requested variables and observed variables with
% Measurement error
-if ~all(diag(StateSpaceModel.measurement_error)==0)
- [observable_pos,index_subset,index_observables]=intersect(SubsetOfVariables,StateSpaceModel.observable_pos,'stable');
- ME_Variance=diag(StateSpaceModel.measurement_error);
+if ~all(diag(M_.H)==0)
+ [observable_pos,index_subset,index_observables]=intersect(SubsetOfVariables,options_.varobs_id,'stable');
+ ME_Variance=diag(M_.H);
ConditionalVarianceDecomposition_ME = zeros(length(observable_pos),length(Steps),number_of_state_innovations+1);
for i=1:number_of_state_innovations
diff --git a/matlab/dsge_simulated_theoretical_conditional_variance_decomposition.m b/matlab/dsge_simulated_theoretical_conditional_variance_decomposition.m
index 9602eafad3d073305b38516c49480e85fd348f89..5918b75e475f62a37f8b7490a3a6993b8e8194e8 100644
--- a/matlab/dsge_simulated_theoretical_conditional_variance_decomposition.m
+++ b/matlab/dsge_simulated_theoretical_conditional_variance_decomposition.m
@@ -19,7 +19,7 @@ function [nvar,vartan,NumberOfConditionalDecompFiles] = ...
% vartan [char] array of characters (with nvar rows).
% NumberOfConditionalDecompFiles [integer] scalar, number of prior or posterior data files (for covariance).
-% Copyright © 2009-2021 Dynare Team
+% Copyright © 2009-2023 Dynare Team
%
% This file is part of Dynare.
%
@@ -46,8 +46,7 @@ elseif strcmpi(type,'prior')
CheckPath('prior/moments',M_.dname);
posterior = 0;
else
- disp('dsge_simulated_theoretical_conditional_variance_decomposition:: Unknown type!')
- error()
+ error('dsge_simulated_theoretical_conditional_variance_decomposition:: Unknown type requested!')
end
%delete old stale files before creating new ones
@@ -82,15 +81,12 @@ options_.ar = 0;
NumberOfSavedElementsPerSimulation = nvar*M_.exo_nbr*length(Steps);
MaXNumberOfConditionalDecompLines = ceil(options_.MaxNumberOfBytes/NumberOfSavedElementsPerSimulation/8);
-ME_present=0;
-if ~all(diag(M_.H)==0)
- [observable_pos_requested_vars,index_subset,index_observables]=intersect(ivar,options_.varobs_id,'stable');
- if ~isempty(observable_pos_requested_vars)
- ME_present=1;
- nobs_ME=length(observable_pos_requested_vars);
- NumberOfSavedElementsPerSimulation_ME = nobs_ME*(M_.exo_nbr+1)*length(Steps);
- MaXNumberOfConditionalDecompLines_ME = ceil(options_.MaxNumberOfBytes/NumberOfSavedElementsPerSimulation_ME/8);
- end
+[ME_present,observable_pos_requested_vars,index_subset,index_observables]=check_measurement_error_requested_vars(M_,options_,ivar);
+
+if ME_present && ~isempty(observable_pos_requested_vars)
+ nobs_ME=length(observable_pos_requested_vars);
+ NumberOfSavedElementsPerSimulation_ME = nobs_ME*(M_.exo_nbr+1)*length(Steps);
+ MaXNumberOfConditionalDecompLines_ME = ceil(options_.MaxNumberOfBytes/NumberOfSavedElementsPerSimulation_ME/8);
end
if SampleSize<=MaXNumberOfConditionalDecompLines
@@ -118,12 +114,6 @@ end
NumberOfConditionalDecompLines = size(Conditional_decomposition_array,4);
ConditionalDecompFileNumber = 0;
-
-StateSpaceModel.number_of_state_equations = M_.endo_nbr;
-StateSpaceModel.number_of_state_innovations = M_.exo_nbr;
-
-first_call = 1;
-
linea = 0;
linea_ME = 0;
for file = 1:NumberOfDrawsFiles
@@ -142,28 +132,12 @@ for file = 1:NumberOfDrawsFiles
dr = temp.pdraws{linee,2};
else
M_=set_parameters_locally(M_,temp.pdraws{linee,1});
- [dr,info,M_,oo_] = compute_decision_rules(M_,options_,oo_)
- end
- if first_call
- endo_nbr = M_.endo_nbr;
- nstatic = M_.nstatic;
- nspred = M_.nspred;
- iv = (1:endo_nbr)';
- ic = [ nstatic+(1:nspred) endo_nbr+(1:size(dr.ghx,2)-nspred) ]';
- StateSpaceModel.number_of_state_equations = M_.endo_nbr;
- StateSpaceModel.number_of_state_innovations = M_.exo_nbr;
- StateSpaceModel.sigma_e_is_diagonal = M_.sigma_e_is_diagonal;
- StateSpaceModel.order_var = dr.order_var;
- StateSpaceModel.observable_pos=options_.varobs_id;
- first_call = 0;
- clear('endo_nbr','nstatic','nspred','k');
+ [dr,info,M_,oo_] = compute_decision_rules(M_,options_,oo_);
end
- [StateSpaceModel.transition_matrix,StateSpaceModel.impulse_matrix] = kalman_transition_matrix(dr,iv,ic,M_.exo_nbr);
- StateSpaceModel.state_innovations_covariance_matrix = M_.Sigma_e;
+
M_ = set_measurement_errors(temp.pdraws{linee,1},temp.estim_params_,M_);
- StateSpaceModel.measurement_error=M_.H;
- clear('dr');
- [ConditionalVarianceDecomposition, ConditionalVarianceDecomposition_ME]=conditional_variance_decomposition(StateSpaceModel, Steps, ivar);
+ [ConditionalVarianceDecomposition, ConditionalVarianceDecomposition_ME] = ...
+ conditional_variance_decomposition(M_,options_,dr, Steps, ivar);
Conditional_decomposition_array(:,:,:,linea) =ConditionalVarianceDecomposition;
if ME_present
Conditional_decomposition_array_ME(:,:,:,linea) =ConditionalVarianceDecomposition_ME;
diff --git a/matlab/dsge_simulated_theoretical_variance_decomposition.m b/matlab/dsge_simulated_theoretical_variance_decomposition.m
index e3de660b5183699c546334a8d48f46dd74dbaad9..305dc76d227d3772a6bb01a809c729c1a68ad134 100644
--- a/matlab/dsge_simulated_theoretical_variance_decomposition.m
+++ b/matlab/dsge_simulated_theoretical_variance_decomposition.m
@@ -56,7 +56,8 @@ if posterior
delete_stale_file([M_.dname '/metropolis/' M_.fname '_PosteriorVarianceDecompME*']);
else
delete_stale_file([M_.dname '/prior/moments/' M_.fname '_PriorVarianceDecomposition*']);
- delete_stale_file([M_.dname '/prior/moments/' M_.fname '_PriorVarianceDecompME*']);end
+ delete_stale_file([M_.dname '/prior/moments/' M_.fname '_PriorVarianceDecompME*']);
+end
% Set varlist (vartan)
if ~posterior
@@ -78,22 +79,17 @@ nvar = length(ivar);
nar = options_.ar;
options_.ar = 0;
-
-
nexo = M_.exo_nbr;
NumberOfSavedElementsPerSimulation = nvar*(nexo+1);
MaXNumberOfDecompLines = ceil(options_.MaxNumberOfBytes/NumberOfSavedElementsPerSimulation/8);
-ME_present=0;
-if ~all(diag(M_.H)==0)
- [observable_pos_requested_vars,index_subset,index_observables]=intersect(ivar,options_.varobs_id,'stable');
- if ~isempty(observable_pos_requested_vars)
- ME_present=1;
- nobs_ME=length(observable_pos_requested_vars);
- NumberOfSavedElementsPerSimulation_ME = nobs_ME*(nexo+1);
- MaXNumberOfDecompLines_ME = ceil(options_.MaxNumberOfBytes/NumberOfSavedElementsPerSimulation_ME/8);
- end
+[ME_present,observable_pos_requested_vars,index_subset,index_observables]=check_measurement_error_requested_vars(M_,options_,ivar);
+
+if ME_present && ~isempty(observable_pos_requested_vars)
+ nobs_ME=length(observable_pos_requested_vars);
+ NumberOfSavedElementsPerSimulation_ME = nobs_ME*(nexo+1);
+ MaXNumberOfDecompLines_ME = ceil(options_.MaxNumberOfBytes/NumberOfSavedElementsPerSimulation_ME/8);
end
if SampleSize<=MaXNumberOfDecompLines
diff --git a/matlab/dyn_first_order_solver.m b/matlab/dyn_first_order_solver.m
index 41d49cbdd543f4ea0ba539dedd37e2915d127229..01226aaeaec4682a5866eaf7dbc49ec1a286981d 100644
--- a/matlab/dyn_first_order_solver.m
+++ b/matlab/dyn_first_order_solver.m
@@ -65,7 +65,6 @@ if isempty(reorder_jacobian_columns)
nstatic = DynareModel.nstatic;
ndynamic = DynareModel.ndynamic;
nsfwrd = DynareModel.nsfwrd;
- n = DynareModel.endo_nbr;
k1 = 1:(npred+nboth);
k2 = 1:(nfwrd+nboth);
@@ -76,16 +75,10 @@ if isempty(reorder_jacobian_columns)
nz = nnz(lead_lag_incidence);
lead_id = find(lead_lag_incidence(maximum_lag+2,:));
- lead_idx = lead_lag_incidence(maximum_lag+2,lead_id);
if maximum_lag
lag_id = find(lead_lag_incidence(1,:));
- lag_idx = lead_lag_incidence(1,lag_id);
- static_id = find((lead_lag_incidence(1,:) == 0) & ...
- (lead_lag_incidence(3,:) == 0));
else
lag_id = [];
- lag_idx = [];
- static_id = find(lead_lag_incidence(2,:)==0);
end
both_id = intersect(lead_id,lag_id);
@@ -277,7 +270,6 @@ if nstatic > 0
end
end
ghx = [temp; ghx];
- temp = [];
end
A_ = real([B_static C*gx+B_pred B_fyd]); % The state_variable of the block are located at [B_pred B_both]
diff --git a/matlab/dynare_config.m b/matlab/dynare_config.m
index c6617a02758fe5e7b832f43d01c2d629c4758c78..c272e937e5cc60591b07831ec9f559eafdaea111 100644
--- a/matlab/dynare_config.m
+++ b/matlab/dynare_config.m
@@ -59,6 +59,7 @@ p = {'/../contrib/jsonlab/' ; ...
'/modules/dseries/src/' ; ...
'/modules/reporting/macros/'; ...
'/modules/reporting/src/' ; ...
+ '/moments/'; ...
'/ms-sbvar/' ; ...
'/ms-sbvar/identification/' ; ...
'/nonlinear-filters/' ; ...
diff --git a/matlab/moments/check_measurement_error_requested_vars.m b/matlab/moments/check_measurement_error_requested_vars.m
new file mode 100644
index 0000000000000000000000000000000000000000..76bee009d90b1b135418bd15a58f90bea7320729
--- /dev/null
+++ b/matlab/moments/check_measurement_error_requested_vars.m
@@ -0,0 +1,45 @@
+function [ME_present,observable_pos_requested_vars,index_subset,index_observables]=check_measurement_error_requested_vars(M_,options_,ivar)
+% [ME_present,observable_pos_requested_vars,index_subset,index_observables]=check_measurement_error_requested_vars(M_,options_,ivar)
+% This function checks for the presence of measurement error and outputs
+% the indices of the affected variables
+%
+% INPUTS
+% M_ [struct] Matlab's structure describing the Model
+% options_ [struct] Matlab's structure describing the options
+% i_var : [double] Index of requested variables in declaration order
+%
+% OUTPUTS
+% ME_present [boolean] indicator whether measurement error is present for requested variables
+% observable_pos_requested_vars [integer] index of observables in list of endogenous variables
+% index_subset [integer] index of observables in ivar
+% index_observables [integer] index of requested i_var in observables
+
+% Copyright © 2023 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 <https://www.gnu.org/licenses/>.
+
+
+index_subset=[];
+index_observables=[];
+observable_pos_requested_vars=[];
+
+ME_present=false;
+if ~all(diag(M_.H)==0)
+ [observable_pos_requested_vars,index_subset,index_observables]=intersect(ivar,options_.varobs_id,'stable');
+ if ~isempty(observable_pos_requested_vars)
+ ME_present=true;
+ end
+end
diff --git a/matlab/moments/compute_variance_decomposition.m b/matlab/moments/compute_variance_decomposition.m
new file mode 100644
index 0000000000000000000000000000000000000000..97b7afde20570ce3ecbc783b54c4e27dee1aab70
--- /dev/null
+++ b/matlab/moments/compute_variance_decomposition.m
@@ -0,0 +1,64 @@
+function var_decomp=compute_variance_decomposition(M_,options_,var_stationary,A,A_stationary,ghu,ghu_states_only,stationary_vars,index_stationary_vars,nvars_requested)
+% function var_decomp=compute_variance_decomposition(M_,options_,var_stationary,A,A_stationary,ghu,ghu_states_only,stationary_vars,index_stationary_vars,nvars_requested)
+% Computes the theoretical auto-covariances, Gamma_y, for an AR(p) process
+% with coefficients dr.ghx and dr.ghu and shock variances Sigma_e
+% for a subset of variables ivar.
+% Theoretical HP-filtering and band-pass filtering is available as an option
+%
+% INPUTS
+% M_ [structure] Global dynare's structure, description of the DSGE model.
+% options_ [structure] Global dynare's structure.
+% var_stationary [double] unconditional variance of stationary
+% variables
+% A [double] State transition matrix
+% A_stationary [double] Reaction of stationary variables to states
+% ghu [double] Reaction of variables to shocks
+% ghu_states_only [double] Reaction of stationary variables to shocks
+% stationary_vars [integer] index of stationary vars in requested output
+% index_stationary_vars [integer] index of stationary vars in decision rules
+% nvars_requested [integer] number of originally requested variables
+%
+% OUTPUTS
+% stationary_vars [double] [#stationary vars by shocks] Matrix containing the variance decomposition
+%
+% Copyright © 2001-2023 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 <https://www.gnu.org/licenses/>.
+
+if M_.exo_nbr == 1
+ var_decomp = ones(size(ghu,1),1);
+else
+ var_decomp = NaN(nvars_requested,M_.exo_nbr);
+ cs = get_lower_cholesky_covariance(M_.Sigma_e,options_.add_tiny_number_to_cholesky);
+ b1 = ghu_states_only*cs;
+ b2 = ghu(index_stationary_vars,:)*cs;
+ variance_sum_loop = 0;
+ for i=1:M_.exo_nbr
+ if i==1 % first time, do Schur decomposition
+ variance_states = lyapunov_symm(A,b1(:,i)*b1(:,i)',options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold,1,options_.debug);
+ else
+ variance_states = lyapunov_symm(A,b1(:,i)*b1(:,i)',options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold,2,options_.debug);
+ end
+ vx2 = diag(A_stationary*variance_states*A_stationary'+b2(:,i)*b2(:,i)');
+ var_decomp(stationary_vars,i) = vx2;
+ variance_sum_loop = variance_sum_loop +vx2; %track overall variance over shocks
+ end
+ if ~options_.pruning && max(abs(variance_sum_loop-var_stationary)./var_stationary) > 1e-4
+ warning(['Aggregate variance and sum of variances by shocks ' ...
+ 'differ by more than 0.01 %'])
+ end
+ var_decomp(stationary_vars,:) = var_decomp(stationary_vars,:)./repmat(variance_sum_loop,1,M_.exo_nbr);
+end
\ No newline at end of file
diff --git a/matlab/disp_th_moments.m b/matlab/moments/disp_th_moments.m
similarity index 58%
rename from matlab/disp_th_moments.m
rename to matlab/moments/disp_th_moments.m
index 46a05223f9c7bd8670118ed03af3036006c7c4dd..01c59e0ddb89f9f130458394a5fe223c99c2dc41 100644
--- a/matlab/disp_th_moments.m
+++ b/matlab/moments/disp_th_moments.m
@@ -1,8 +1,29 @@
function oo_ = disp_th_moments(dr, var_list, M_, options_, oo_)
-
+%oo_ = disp_th_moments(dr, var_list, M_, options_, oo_)
% Display theoretical moments of variables
+% INPUTS:
+% dr : [struct] Dynare decision rules structure
+% var_list [cell] list of variables considered
+% M_ [struct] structure describing the Model
+% options_ [struct] structure describing the options
+% oo_ [struct] structure describing the Model
+%
+% OUTPUTS:
+% oo_ [struct] structure describing the Model, containing
+% gamma_y [cell] Matlab cell of nar+1 arrays, where nar is the order of the autocorrelation function.
+% gamma_y{1} [double] Covariance matrix.
+% gamma_y{i+1} [double] Autocorrelation function (for i=1,...,options_.ar).
+% mean [vector] Unconditional mean
+% var [matrix] Unconditional covariance matrix
+% autocorr [cell] Cell storing the theoretical autocorrelation
+% contemporaneous_correlation [matrix] matrix of contemporaneous correlations
+% autocorr [cell] Cell storing the theoretical autocorrelation
+% variance_decomposition [matrix] Unconditional variance decomposition matrix
+% variance_decomposition_ME [matrix] Unconditional variance decomposition matrix with measurement error
+% conditional_variance_decomposition [array] Conditional variance decomposition array
+% conditional_variance_decomposition_ME [array] Conditional variance decomposition array with measurement error
-% Copyright © 2001-2021 Dynare Team
+% Copyright © 2001-2023 Dynare Team
%
% This file is part of Dynare.
%
@@ -53,24 +74,9 @@ z = [ m sd s2 ];
oo_.mean = m;
oo_.var = oo_.gamma_y{1};
-ME_present=0;
-if ~all(diag(M_.H)==0)
- [observable_pos_requested_vars,index_subset,index_observables]=intersect(ivar,options_.varobs_id,'stable');
- if ~isempty(observable_pos_requested_vars)
- ME_present=1;
- end
-end
if size(stationary_vars, 1) > 0
- if ~nodecomposition
- oo_.variance_decomposition=100*oo_.gamma_y{options_.ar+2};
- if ME_present
- ME_Variance=diag(M_.H);
- oo_.variance_decomposition_ME=oo_.variance_decomposition(index_subset,:).*repmat(diag(oo_.var(index_subset,index_subset))./(diag(oo_.var(index_subset,index_subset))+ME_Variance(index_observables)),1,M_.exo_nbr);
- oo_.variance_decomposition_ME(:,end+1)=100-sum(oo_.variance_decomposition_ME,2);
- end
- end
- if ~options_.noprint %options_.nomoments == 0
+ if ~options_.noprint
if options_.order == 2
title = 'APPROXIMATED THEORETICAL MOMENTS';
else
@@ -86,61 +92,34 @@ if size(stationary_vars, 1) > 0
lh = cellofchararraymaxlength(labels)+2;
dyn_latex_table(M_, options_, title, 'th_moments', headers, labels, z, lh, 11, 4);
end
+ end
+ [ME_present,observable_pos_requested_vars,index_subset,index_observables]=check_measurement_error_requested_vars(M_,options_,ivar);
+ %store unconditional variance decomposition
+ if ~nodecomposition
+ oo_.variance_decomposition=100*oo_.gamma_y{options_.ar+2};
+ if ME_present
+ ME_Variance=diag(M_.H);
+ oo_.variance_decomposition_ME=oo_.variance_decomposition(index_subset,:).*repmat(diag(oo_.var(index_subset,index_subset))./(diag(oo_.var(index_subset,index_subset))+ME_Variance(index_observables)),1,M_.exo_nbr);
+ oo_.variance_decomposition_ME(:,end+1)=100-sum(oo_.variance_decomposition_ME,2);
+ end
+ end
+ if ~options_.noprint %options_.nomoments == 0
if M_.exo_nbr > 1 && ~nodecomposition
- skipline()
- if options_.order == 2
- title = 'APPROXIMATED VARIANCE DECOMPOSITION (in percent)';
- else
- title = 'VARIANCE DECOMPOSITION (in percent)';
- end
- title = add_filter_subtitle(title, options_);
- headers = M_.exo_names;
- headers = vertcat(' ', headers);
- labels=get_labels_transformed_vars(M_.endo_names,ivar(stationary_vars),options_,false);
- lh = cellofchararraymaxlength(labels)+2;
- dyntable(options_, title, headers, labels, 100*oo_.gamma_y{options_.ar+2}(stationary_vars,:), lh, 8, 2);
- if ME_present
- [stationary_observables, pos_index_subset] = intersect(index_subset, stationary_vars, 'stable');
- headers_ME = vertcat(headers, 'ME');
- labels=get_labels_transformed_vars(M_.endo_names,ivar(stationary_observables),options_,false);
- dyntable(options_, [title,' WITH MEASUREMENT ERROR'], headers_ME, labels, ...
- oo_.variance_decomposition_ME(pos_index_subset,:), lh, 8, 2);
- end
- if options_.TeX
- headers = M_.exo_names_tex;
- headers = vertcat(' ', headers);
- labels=get_labels_transformed_vars(M_.endo_names_tex,ivar(stationary_vars),options_,true);
- lh = cellofchararraymaxlength(labels)+2;
- dyn_latex_table(M_, options_, title, 'th_var_decomp_uncond', headers, labels, 100*oo_.gamma_y{options_.ar+2}(stationary_vars,:), lh, 8, 2);
- if ME_present
- headers_ME = vertcat(headers, 'ME');
- labels=get_labels_transformed_vars(M_.endo_names_tex,ivar(stationary_observables),options_,true);
- dyn_latex_table(M_, options_, [title,' WITH MEASUREMENT ERROR'], ...
- 'th_var_decomp_uncond_ME', headers_ME, labels, oo_.variance_decomposition_ME(pos_index_subset,:), lh, 8, 2);
- end
- end
+ display_unconditional_variance_decomposition(M_,options_,oo_,ivar,stationary_vars,index_subset,ME_present)
end
end
- conditional_variance_steps = options_.conditional_variance_decomposition;
- if ~isempty(conditional_variance_steps)
- StateSpaceModel.number_of_state_equations = M_.endo_nbr;
- StateSpaceModel.number_of_state_innovations = M_.exo_nbr;
- StateSpaceModel.sigma_e_is_diagonal = M_.sigma_e_is_diagonal;
- [StateSpaceModel.transition_matrix, StateSpaceModel.impulse_matrix] = ...
- kalman_transition_matrix(dr,(1:M_.endo_nbr)',M_.nstatic+(1:M_.nspred)',M_.exo_nbr);
- StateSpaceModel.state_innovations_covariance_matrix = M_.Sigma_e;
- StateSpaceModel.order_var = dr.order_var;
- StateSpaceModel.measurement_error = M_.H;
- StateSpaceModel.observable_pos = options_.varobs_id;
- [oo_.conditional_variance_decomposition, oo_.conditional_variance_decomposition_ME] = ...
- conditional_variance_decomposition(StateSpaceModel, conditional_variance_steps, ivar);
- if ~options_.noprint
- display_conditional_variance_decomposition(oo_.conditional_variance_decomposition, conditional_variance_steps, ivar, M_, options_);
- if ME_present
- display_conditional_variance_decomposition(oo_.conditional_variance_decomposition_ME, conditional_variance_steps, ...
- observable_pos_requested_vars, M_, options_);
- end
+end
+%% Conditional variance decomposition
+conditional_variance_steps = options_.conditional_variance_decomposition;
+if ~isempty(conditional_variance_steps)
+ [oo_.conditional_variance_decomposition, oo_.conditional_variance_decomposition_ME] = ...
+ conditional_variance_decomposition(M_,options_,dr, conditional_variance_steps, ivar);
+ if ~options_.noprint
+ display_conditional_variance_decomposition(oo_.conditional_variance_decomposition, conditional_variance_steps, ivar, M_, options_);
+ if ME_present
+ display_conditional_variance_decomposition(oo_.conditional_variance_decomposition_ME, conditional_variance_steps, ...
+ observable_pos_requested_vars, M_, options_);
end
end
end
@@ -195,13 +174,13 @@ if options_.ar > 0 && size(stationary_vars, 1) > 0
title = 'COEFFICIENTS OF AUTOCORRELATION';
end
title = add_filter_subtitle(title, options_);
- labels=get_labels_transformed_vars(M_.endo_names,ivar(i1),options_,false);
- headers = vertcat('Order ', cellstr(int2str([1:options_.ar]')));
+ labels=get_labels_transformed_vars(M_.endo_names,ivar(i1),options_,false);
+ headers = vertcat('Order ', cellstr(int2str((1:options_.ar)')));
lh = cellofchararraymaxlength(labels)+2;
dyntable(options_, title, headers, labels, z, lh, 8, 4);
if options_.TeX
labels=get_labels_transformed_vars(M_.endo_names_tex,ivar(i1),options_,true);
- headers = vertcat('Order ', cellstr(int2str([1:options_.ar]')));
+ headers = vertcat('Order ', cellstr(int2str((1:options_.ar)')));
lh = cellofchararraymaxlength(labels)+2;
dyn_latex_table(M_, options_, title, 'th_autocorr_matrix', headers, labels, z, lh, 8, 4);
end
diff --git a/matlab/disp_th_moments_pruned_state_space.m b/matlab/moments/disp_th_moments_pruned_state_space.m
similarity index 51%
rename from matlab/disp_th_moments_pruned_state_space.m
rename to matlab/moments/disp_th_moments_pruned_state_space.m
index 5675fabe3f8bbf8da35d94970ca1c1efa22b964c..8db2dad713d9b03f2aae2e0aec0f4b047d8b4757 100644
--- a/matlab/disp_th_moments_pruned_state_space.m
+++ b/matlab/moments/disp_th_moments_pruned_state_space.m
@@ -4,23 +4,27 @@ function oo_=disp_th_moments_pruned_state_space(dr,M_,options_,i_var,oo_)
% pruned state-space
%
% INPUTS:
-% dr : Dynare decision rules structure
-% M_ : Matlab's structure describing the Model (initialized by dynare, see @ref{M_}).
-% options_ : Matlab's structure describing the options (initialized by dynare, see @ref{options_}).
-% i_var : Index of requested variables in declaration order
-% oo_ : Matlab's structure describing the Model (initialized by dynare, see @ref{M_}), containing
+% dr : [struct] Dynare decision rules structure
+% M_ [struct] structure describing the Model
+% options_ [struct] structure describing the options
+% i_var [double] Index of requested variables in declaration order
+% oo_ [struct] structure describing the Model
%
% OUTPUTS:
-% oo_ : Matlab's structure describing the Model (initialized by dynare, see @ref{M_}), containing
-% gamma_y [cell] Matlab cell of nar+1 arrays, where nar is the order of the autocorrelation function.
-% gamma_y{1} [double] Covariance matrix.
-% gamma_y{i+1} [double] Autocorrelation function (for i=1,...,options_.ar).
-% mean [vector] Unconditional mean
-% var [matrix] Unconditional covariance matrix
-% autocorr [cell] Cell storing the theoretical autocorrelation
-% contemporaneous_correlation [matrix] matrix of contemporaneous correlations
-%
-% Copyright © 2020 Dynare Team
+% gamma_y [cell] Matlab cell of nar+1 arrays, where nar is the order of the autocorrelation function.
+% gamma_y{1} [double] Covariance matrix.
+% gamma_y{i+1} [double] Autocorrelation function (for i=1,...,options_.ar).
+% mean [vector] Unconditional mean
+% var [matrix] Unconditional covariance matrix
+% autocorr [cell] Cell storing the theoretical autocorrelation
+% contemporaneous_correlation [matrix] matrix of contemporaneous correlations
+% autocorr [cell] Cell storing the theoretical autocorrelation
+% variance_decomposition [matrix] Unconditional variance decomposition matrix
+% variance_decomposition_ME [matrix] Unconditional variance decomposition matrix with measurement error
+% conditional_variance_decomposition [array] Conditional variance decomposition array
+% conditional_variance_decomposition_ME [array] Conditional variance decomposition array with measurement error
+
+% Copyright © 2020-2023 Dynare Team
%
% This file is part of Dynare.
%
@@ -123,4 +127,38 @@ if options_.ar > 0 %&& size(stationary_vars, 1) > 0
dyn_latex_table(M_,options_,title,'th_autocorr_matrix',headers,labels,z,lh,8,4);
end
end
+end
+
+if options_.order==2 && ~options_.nodecomposition && M_.exo_nbr > 1% do variance decomposition
+ index_stationary_vars = dr.inv_order_var(i_var); %no nonstationary variables in pruning
+ ghu_states_only = zeros(M_.nspred,M_.exo_nbr);
+ ghu_states_only(1:M_.nspred,:) = dr.ghu(M_.nstatic+(1:M_.nspred),:); %get shock impact on states only
+ [A] = kalman_transition_matrix(dr,M_.nstatic+(1:M_.nspred)',1:M_.nspred);
+ oo_.gamma_y{options_.ar+2}=compute_variance_decomposition(M_,options_,s2,A,dr.ghx(index_stationary_vars,:),dr.ghu,ghu_states_only,1:length(i_var),index_stationary_vars,nvars);
+
+ [ME_present,observable_pos_requested_vars,index_subset,index_observables]=check_measurement_error_requested_vars(M_,options_,i_var);
+ %store unconditional variance decomposition
+ oo_.variance_decomposition=100*oo_.gamma_y{options_.ar+2};
+ if ME_present
+ ME_Variance=diag(M_.H);
+ oo_.variance_decomposition_ME=oo_.variance_decomposition(index_subset,:).*repmat(diag(oo_.var(index_subset,index_subset))./(diag(oo_.var(index_subset,index_subset))+ME_Variance(index_observables)),1,M_.exo_nbr);
+ oo_.variance_decomposition_ME(:,end+1)=100-sum(oo_.variance_decomposition_ME,2);
+ end
+ if ~options_.noprint %options_.nomoments == 0
+ display_unconditional_variance_decomposition(M_,options_,oo_,i_var,1:length(i_var),index_subset,ME_present)
+ end
+
+ %% Conditional variance decomposition
+ conditional_variance_steps = options_.conditional_variance_decomposition;
+ if ~isempty(conditional_variance_steps)
+ [oo_.conditional_variance_decomposition, oo_.conditional_variance_decomposition_ME] = ...
+ conditional_variance_decomposition(M_,options_,dr, conditional_variance_steps, i_var);
+ if ~options_.noprint
+ display_conditional_variance_decomposition(oo_.conditional_variance_decomposition, conditional_variance_steps, i_var, M_, options_);
+ if ME_present
+ display_conditional_variance_decomposition(oo_.conditional_variance_decomposition_ME, conditional_variance_steps, ...
+ observable_pos_requested_vars, M_, options_);
+ end
+ end
+ end
end
\ No newline at end of file
diff --git a/matlab/display_conditional_variance_decomposition.m b/matlab/moments/display_conditional_variance_decomposition.m
similarity index 100%
rename from matlab/display_conditional_variance_decomposition.m
rename to matlab/moments/display_conditional_variance_decomposition.m
diff --git a/matlab/moments/display_unconditional_variance_decomposition.m b/matlab/moments/display_unconditional_variance_decomposition.m
new file mode 100644
index 0000000000000000000000000000000000000000..df547d02cf12d4083605c385c3905d20f8ca921a
--- /dev/null
+++ b/matlab/moments/display_unconditional_variance_decomposition.m
@@ -0,0 +1,67 @@
+function display_unconditional_variance_decomposition(M_,options_,oo_,ivar,stationary_vars,index_subset,ME_present)
+% display_unconditional_variance_decomposition(M_,options_,oo_,ivar,stationary_vars,index_subset,ME_present)
+% This function displays the unconditional variance decomposition
+%
+% INPUTS
+% M_ [struct] Matlab's structure describing the Model (initialized by dynare, see @ref{M_}).
+% options_ [struct] Matlab's structure describing the options (initialized by dynare, see @ref{options_}).
+% oo_ [struct] structure describing the Model
+% i_var [double] Index of requested variables in declaration order
+% stationary_vars [double] index of stationary vars in requested output
+% index_subset [integer] index of observables in requested variables
+% ME_present [boolean] indicator whether measurement error is present for requested variables
+%
+% OUTPUTS
+% None
+
+% Copyright © 2023 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 <https://www.gnu.org/licenses/>.
+
+if M_.exo_nbr > 1
+ skipline()
+ if options_.order == 2
+ title = 'APPROXIMATED VARIANCE DECOMPOSITION (in percent)';
+ else
+ title = 'VARIANCE DECOMPOSITION (in percent)';
+ end
+ title = add_filter_subtitle(title, options_);
+ headers = M_.exo_names;
+ headers = vertcat(' ', headers);
+ labels=get_labels_transformed_vars(M_.endo_names,ivar(stationary_vars),options_,false);
+ lh = cellofchararraymaxlength(labels)+2;
+ dyntable(options_, title, headers, labels, 100*oo_.gamma_y{options_.ar+2}(stationary_vars,:), lh, 8, 2);
+ if ME_present
+ [stationary_observables, pos_index_subset] = intersect(index_subset, stationary_vars, 'stable');
+ headers_ME = vertcat(headers, 'ME');
+ labels=get_labels_transformed_vars(M_.endo_names,ivar(stationary_observables),options_,false);
+ dyntable(options_, [title,' WITH MEASUREMENT ERROR'], headers_ME, labels, ...
+ oo_.variance_decomposition_ME(pos_index_subset,:), lh, 8, 2);
+ end
+ if options_.TeX
+ headers = M_.exo_names_tex;
+ headers = vertcat(' ', headers);
+ labels=get_labels_transformed_vars(M_.endo_names_tex,ivar(stationary_vars),options_,true);
+ lh = cellofchararraymaxlength(labels)+2;
+ dyn_latex_table(M_, options_, title, 'th_var_decomp_uncond', headers, labels, 100*oo_.gamma_y{options_.ar+2}(stationary_vars,:), lh, 8, 2);
+ if ME_present
+ headers_ME = vertcat(headers, 'ME');
+ labels=get_labels_transformed_vars(M_.endo_names_tex,ivar(stationary_observables),options_,true);
+ dyn_latex_table(M_, options_, [title,' WITH MEASUREMENT ERROR'], ...
+ 'th_var_decomp_uncond_ME', headers_ME, labels, oo_.variance_decomposition_ME(pos_index_subset,:), lh, 8, 2);
+ end
+ end
+end
\ No newline at end of file
diff --git a/matlab/th_autocovariances.m b/matlab/moments/th_autocovariances.m
similarity index 85%
rename from matlab/th_autocovariances.m
rename to matlab/moments/th_autocovariances.m
index dbcbcd56e42ad50ee2e8ab92310528a1b486f6b6..a3255e237a27288b18443185e36f40ff8bce0ff6 100644
--- a/matlab/th_autocovariances.m
+++ b/matlab/moments/th_autocovariances.m
@@ -100,16 +100,16 @@ ghu_states_only(1:M_.nspred,:) = ghu(index_states,:); %get shock impact on state
% Compute stationary variables for unfiltered moments (filtering will remove unit roots)
if options_.hp_filter ~= 0 || options_.bandpass.indicator
% By construction, all variables are stationary when filtered
- iky = inv_order_var(ivar);
+ index_stationary_vars = inv_order_var(ivar);
stationary_vars = (1:length(ivar))';
else
[variance_states, unit_root_Schur_vector] = lyapunov_symm(A,B*M_.Sigma_e*B',options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold,[],options_.debug);
- iky = inv_order_var(ivar);
+ index_stationary_vars = inv_order_var(ivar);
stationary_vars = (1:length(ivar))';
if ~isempty(unit_root_Schur_vector)
x = abs(ghx*unit_root_Schur_vector);
- iky = iky(all(x(iky,:) < options_.schur_vec_tol,2));
stationary_vars = find(all(x(inv_order_var(ivar),:) < options_.schur_vec_tol,2));
+ index_stationary_vars = inv_order_var(ivar(stationary_vars));
end
end
@@ -119,17 +119,17 @@ if local_order == 2 % mean correction for 2nd order with no filters; oth
Ex = (dr.ghs2(index_states)+dr.ghxx(index_states,:)*variance_states(:)+dr.ghuu(index_states,:)*M_.Sigma_e(:))/2;
Ex = (eye(length(M_.nspred))-ghx(index_states,:))\Ex;
Gamma_y{nar+3} = NaN*ones(nvar, 1);
- Gamma_y{nar+3}(stationary_vars) = ghx(iky,:)*Ex+(dr.ghs2(iky)+dr.ghxx(iky,:)*variance_states(:)+...
- dr.ghuu(iky,:)*M_.Sigma_e(:))/2;
+ Gamma_y{nar+3}(stationary_vars) = ghx(index_stationary_vars,:)*Ex+(dr.ghs2(index_stationary_vars)+dr.ghxx(index_stationary_vars,:)*variance_states(:)+...
+ dr.ghuu(index_stationary_vars,:)*M_.Sigma_e(:))/2;
else %no static and no predetermined
Gamma_y{nar+3} = NaN*ones(nvar, 1);
- Gamma_y{nar+3}(stationary_vars) = (dr.ghs2(iky)+ dr.ghuu(iky,:)*M_.Sigma_e(:))/2;
+ Gamma_y{nar+3}(stationary_vars) = (dr.ghs2(index_stationary_vars)+ dr.ghuu(index_stationary_vars,:)*M_.Sigma_e(:))/2;
end
end
if options_.hp_filter == 0 && ~options_.bandpass.indicator
- aa = ghx(iky,:);
- bb = ghu(iky,:);
+ aa = ghx(index_stationary_vars,:);
+ bb = ghu(index_stationary_vars,:);
% unconditional variance
v = NaN*ones(nvar,nvar);
v(stationary_vars,stationary_vars) = aa*variance_states*aa'+ bb*M_.Sigma_e*bb';
@@ -151,35 +151,11 @@ if options_.hp_filter == 0 && ~options_.bandpass.indicator
end
% variance decomposition
if ~nodecomposition && M_.exo_nbr > 0 && size(stationary_vars, 1) > 0
- if M_.exo_nbr == 1
- Gamma_y{nar+2} = ones(nvar,1);
- else
- Gamma_y{nar+2} = NaN(nvar,M_.exo_nbr);
- cs = get_lower_cholesky_covariance(M_.Sigma_e,options_.add_tiny_number_to_cholesky);
- b1 = ghu_states_only*cs;
- b2 = ghu(iky,:)*cs;
- variance_states = lyapunov_symm(A,b1*b1',options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold,1,options_.debug);
- vv = diag(aa*variance_states*aa'+b2*b2');
- variance_sum_loop = 0;
- for i=1:M_.exo_nbr
- variance_states = lyapunov_symm(A,b1(:,i)*b1(:,i)',options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold,2,options_.debug);
- vx2 = diag(aa*variance_states*aa'+b2(:,i)*b2(:,i)');
- Gamma_y{nar+2}(stationary_vars,i) = vx2;
- variance_sum_loop = variance_sum_loop +vx2; %track overall variance over shocks
- end
- if max(abs(variance_sum_loop-vv)./vv) > 1e-4
- warning(['Aggregate variance and sum of variances by shocks ' ...
- 'differ by more than 0.01 %'])
- end
- for i=1:M_.exo_nbr
- Gamma_y{nar+2}(stationary_vars,i) = Gamma_y{nar+ ...
- 2}(stationary_vars,i)./variance_sum_loop;
- end
- end
+ Gamma_y{nar+2}=compute_variance_decomposition(M_,options_,diag(v(stationary_vars,stationary_vars)),A,aa,ghu,ghu_states_only,stationary_vars,index_stationary_vars,nvar);
end
else% ==> Theoretical filters.
- aa = ghx(iky,:); %R in Uhlig (2001)
- bb = ghu(iky,:); %S in Uhlig (2001)
+ aa = ghx(index_stationary_vars,:); %R in Uhlig (2001)
+ bb = ghu(index_stationary_vars,:); %S in Uhlig (2001)
ngrid = options_.filtered_theoretical_moments_grid;
freqs = 0 : ((2*pi)/ngrid) : (2*pi*(1 - .5/ngrid)); %[0,2*pi)
@@ -230,7 +206,7 @@ else% ==> Theoretical filters.
cs = get_lower_cholesky_covariance(M_.Sigma_e); %make sure Covariance matrix is positive definite
SS = cs*cs';
b1 = ghu_states_only;
- b2 = ghu(iky,:);
+ b2 = ghu(index_stationary_vars,:);
mathp_col = NaN(ngrid,length(ivar)^2);
IA = eye(size(A,1));
IE = eye(M_.exo_nbr);