diff --git a/matlab/gsa/map_ident_.m b/matlab/gsa/map_ident_.m
index 2a730faa3053ba1ab91c0a2e5faac80d9674ae54..a8cbaae3cc767b3f219ae96d1def8c75f82d19d5 100644
--- a/matlab/gsa/map_ident_.m
+++ b/matlab/gsa/map_ident_.m
@@ -60,7 +60,7 @@ if opt_gsa.load_ident_files==0
mss = teff(mss(:,istable),Nsam,istable);
yys = teff(yys(dr.order_var,istable),Nsam,istable);
if exist('T','var')
- [vdec, cc, ac] = mc_moments(T, lpmatx, dr);
+ [vdec, cc, ac] = mc_moments(T, lpmatx, dr, M_.exo_nbr, options_);
else
return
end
@@ -77,14 +77,14 @@ if opt_gsa.load_ident_files==0
ifig=0;
for j=1:M_.exo_nbr
if mod(j,6)==1
- hh_fig=dyn_figure(options_.nodisplay,'name',['Variance decomposition shocks']);
+ hh_fig=dyn_figure(options_.nodisplay,'name','Variance decomposition shocks');
ifig=ifig+1;
iplo=0;
end
iplo=iplo+1;
subplot(2,3,iplo)
myboxplot(squeeze(vdec(:,j,:))',[],'.',[],10)
- set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:size(options_.varobs,1)])
+ set(gca,'xticklabel',' ','fontsize',10,'xtick',1:size(options_.varobs,1))
set(gca,'xlim',[0.5 size(options_.varobs,1)+0.5])
set(gca,'ylim',[-2 102])
for ip=1:size(options_.varobs,1)
@@ -95,7 +95,7 @@ if opt_gsa.load_ident_files==0
title(M_.exo_names{j},'interpreter','none')
if mod(j,6)==0 || j==M_.exo_nbr
dyn_saveas(hh_fig,[OutputDirectoryName,'/',fname_,'_vdec_exo_',int2str(ifig)],options_.nodisplay,options_.graph_format);
- create_TeX_loader(options_,[OutputDirectoryName,'/',fname_,'_vdec_exo_',int2str(ifig)],ifig,['Variance decomposition shocks'],'vdec_exo',options_.figures.textwidth*min(iplo/3,1))
+ create_TeX_loader(options_,[OutputDirectoryName,'/',fname_,'_vdec_exo_',int2str(ifig)],ifig,'Variance decomposition shocks','vdec_exo',options_.figures.textwidth*min(iplo/3,1))
end
end
end
@@ -113,8 +113,8 @@ if opt_gsa.load_ident_files==0
iv = (1:endo_nbr)';
ic = [ nstatic+(1:nspred) endo_nbr+(1:size(dr.ghx,2)-nspred) ]';
- dr.ghx = T(:, [1:(nc1-M_.exo_nbr)],1);
- dr.ghu = T(:, [(nc1-M_.exo_nbr+1):end], 1);
+ dr.ghx = T(:, 1:(nc1-M_.exo_nbr),1);
+ dr.ghu = T(:, (nc1-M_.exo_nbr+1):end, 1);
[Aa,Bb] = kalman_transition_matrix(dr,iv,ic);
A = zeros(size(Aa,1),size(Aa,2)+size(Aa,1),length(istable));
if ~isempty(lpmatx)
@@ -122,8 +122,8 @@ if opt_gsa.load_ident_files==0
end
A(:,:,1)=[Aa, triu(Bb*M_.Sigma_e*Bb')];
for j=2:length(istable)
- dr.ghx = T(:, [1:(nc1-M_.exo_nbr)],j);
- dr.ghu = T(:, [(nc1-M_.exo_nbr+1):end], j);
+ dr.ghx = T(:, 1:(nc1-M_.exo_nbr),j);
+ dr.ghu = T(:, (nc1-M_.exo_nbr+1):end, j);
[Aa,Bb] = kalman_transition_matrix(dr, iv, ic);
if ~isempty(lpmatx)
set_shocks_param(lpmatx(j,:));
@@ -155,12 +155,12 @@ if opt_gsa.morris==1
SAvdec = squeeze(SAMorris(:,1,:))';
save([OutputDirectoryName,'/',fname_,'_morris_IDE.mat'],'SAvdec','vdec','ir_vdec','ic_vdec')
else
- load([OutputDirectoryName,'/',fname_,'_morris_IDE.mat'],'SAvdec','vdec','ir_vdec','ic_vdec')
+ load([OutputDirectoryName,'/',fname_,'_morris_IDE.mat'],'SAvdec')
end
hh_fig = dyn_figure(options_.nodisplay,'name','Screening identification: variance decomposition');
myboxplot(SAvdec,[],'.',[],10)
- set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:npT])
+ set(gca,'xticklabel',' ','fontsize',10,'xtick',1:npT)
set(gca,'xlim',[0.5 npT+0.5])
ydum = get(gca,'ylim');
set(gca,'ylim',[0 ydum(2)])
@@ -193,7 +193,7 @@ if opt_gsa.morris==1
hh_fig=dyn_figure(options_.nodisplay,'name','Screening identification: theoretical moments');
myboxplot(SAcc,[],'.',[],10)
- set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:npT])
+ set(gca,'xticklabel',' ','fontsize',10,'xtick',1:npT)
set(gca,'xlim',[0.5 npT+0.5])
set(gca,'ylim',[0 1])
set(gca,'position',[0.13 0.2 0.775 0.7])
@@ -231,11 +231,11 @@ if opt_gsa.morris==1
end
save([OutputDirectoryName,'/',fname_,'_morris_IDE.mat'],'SAnorm','SAmunorm','SAsignorm','-append')
else
- load([OutputDirectoryName,'/',fname_,'_morris_IDE'],'SAnorm','SAmunorm','SAsignorm')
+ load([OutputDirectoryName,'/',fname_,'_morris_IDE'],'SAnorm')
end
hh_fig=dyn_figure(options_.nodisplay,'name','Screening identification: model');
myboxplot(SAnorm',[],'.',[],10)
- set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:npT])
+ set(gca,'xticklabel',' ','fontsize',10,'xtick',1:npT)
set(gca,'xlim',[0.5 npT+0.5])
set(gca,'ylim',[0 1])
set(gca,'position',[0.13 0.2 0.775 0.7])
@@ -261,7 +261,7 @@ else % main effects analysis
fsuffix = '_rank';
end
- imap=[1:npT];
+ imap=1:npT;
if isempty(lpmat0)
x0=lpmat(istable,:);
@@ -309,12 +309,12 @@ else % main effects analysis
save([OutputDirectoryName,'/map_cc',fsuffix,'.mat'],'gsa_')
save([OutputDirectoryName,'/',fname_,'_main_eff.mat'],'imap_cc','SAcc','ccac','-append')
else
- load([OutputDirectoryName,'/',fname_,'_main_eff'],'imap_cc','SAcc','ccac')
+ load([OutputDirectoryName,'/',fname_,'_main_eff'],'SAcc')
end
hh_fig=dyn_figure(options_.nodisplay,'Name',['Identifiability indices in the ',fsuffix,' moments.']);
bar(sum(SAcc))
- set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:npT])
+ set(gca,'xticklabel',' ','fontsize',10,'xtick',1:npT)
set(gca,'xlim',[0.5 npT+0.5])
ydum = get(gca,'ylim');
set(gca,'ylim',[0 ydum(2)])
diff --git a/matlab/gsa/mc_moments.m b/matlab/gsa/mc_moments.m
index e14b39d48130787597f7f8126be3fd99b3b62264..5a04e34a32e38a79c0a5230d7eccafbf96e46fc1 100644
--- a/matlab/gsa/mc_moments.m
+++ b/matlab/gsa/mc_moments.m
@@ -1,6 +1,18 @@
-function [vdec, cc, ac] = mc_moments(mm, ss, dr)
+function [vdec, cc, ac] = mc_moments(mm, ss, dr, exo_nbr, options_)
+% [vdec, cc, ac] = mc_moments(mm, ss, dr, exo_nbr, options_)
+% Conduct Monte Carlo simulation of second moments for GSA
+% Inputs:
+% - dr [structure] decision rules
+% - exo_nbr [double] number of exogenous shocks
+% - options_ [structure] Matlab's structure describing the current options
+%
+% Outputs:
+% - vdec [double] variance decomposition matrix
+% - cc [double] vector of unique elements of cross correlation matrix
+% - ac [cell] autocorrelation matrix
+
-% Copyright © 2012-2018 Dynare Team
+% Copyright © 2012-2023 Dynare Team
%
% This file is part of Dynare.
%
@@ -17,27 +29,25 @@ function [vdec, cc, ac] = mc_moments(mm, ss, dr)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
-global options_ M_ estim_params_ oo_
-
-[nr1, nc1, nsam] = size(mm);
+[~, nc1, nsam] = size(mm);
nobs=length(options_.varobs);
-disp('Computing theoretical moments ...')
+disp('mc_moments: Computing theoretical moments ...')
h = dyn_waitbar(0,'Theoretical moments ...');
-vdec = zeros(nobs,M_.exo_nbr,nsam);
+vdec = zeros(nobs,exo_nbr,nsam);
cc = zeros(nobs,nobs,nsam);
ac = zeros(nobs,nobs*options_.ar,nsam);
for j=1:nsam
- oo_.dr.ghx = mm(:, [1:(nc1-M_.exo_nbr)],j);
- oo_.dr.ghu = mm(:, [(nc1-M_.exo_nbr+1):end], j);
+ dr.ghx = mm(:, 1:(nc1-exo_nbr),j);
+ dr.ghu = mm(:, (nc1-exo_nbr+1):end, j);
if ~isempty(ss)
set_shocks_param(ss(j,:));
end
- [vdec(:,:,j), corr, autocorr, z, zz] = th_moments(oo_.dr,options_.varobs);
+ [vdec(:,:,j), corr, autocorr] = th_moments(dr,options_,M_);
cc(:,:,j)=triu(corr);
- dum=[];
+ dum=NaN(nobs,nobs*options_.ar);
for i=1:options_.ar
- dum=[dum, autocorr{i}];
+ dum(:,(i-1)*nobs+1:i*nobs)=autocorr{i};
end
ac(:,:,j)=dum;
dyn_waitbar(j/nsam,h)
diff --git a/matlab/gsa/th_moments.m b/matlab/gsa/th_moments.m
index 98afaad7e087ab80e48361c310be0f2ef8416fed..fc5cb3de28a4d668243ed467e44e0cebdabfb4b7 100644
--- a/matlab/gsa/th_moments.m
+++ b/matlab/gsa/th_moments.m
@@ -1,7 +1,22 @@
-function [vdec, corr, autocorr, z, zz] = th_moments(dr,var_list)
-% [vdec, corr, autocorr, z, zz] = th_moments(dr,var_list)
+function [vdec, corr, autocorr, z, zz] = th_moments(dr,options_,M_)
+% [vdec, corr, autocorr, z, zz] = th_moments(dr,options_,M_)
+% Computes theoretical moments for GSA
+%
+% INPUTS
+% - dr [structure] model information structure
+% - options_ [structure] Matlab's structure describing the current options
+% - M_ [structure] Matlab's structure describing the model
+%
+% OUTPUTS
+% - vdec [double] variance decomposition matrix
+% - corr [double] correlation matrix
+% - autocorr [cell] contains autocorrelation or
+% auto- and cross-covariance matrices
+% - z [double] matrix containing mean, standard
+% deviation, and variance vector
+% - zz [double] autocorrelation matrix
-% Copyright © 2012-2018 Dynare Team
+% Copyright © 2012-2023 Dynare Team
%
% This file is part of Dynare.
%
@@ -18,18 +33,16 @@ function [vdec, corr, autocorr, z, zz] = th_moments(dr,var_list)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
-global M_ oo_ options_
-
-nvar = length(var_list);
+nvar = length(options_.var_list);
if nvar == 0
nvar = length(dr.order_var);
ivar = [1:nvar]';
else
ivar=zeros(nvar,1);
for i=1:nvar
- i_tmp = strmatch(var_list{i}, M_.endo_names, 'exact');
+ i_tmp = strmatch(options_.var_list{i}, M_.endo_names, 'exact');
if isempty(i_tmp)
- error(['One of the variables specified does not exist']) ;
+ error('th_moments: One of the variables specified does not exist');
else
ivar(i) = i_tmp;
end
@@ -39,21 +52,11 @@ end
[gamma_y,stationary_vars] = th_autocovariances(dr,ivar,M_, options_);
m = dr.ys(ivar(stationary_vars));
-
-% i1 = find(abs(diag(gamma_y{1})) > 1e-12);
-i1 = [1:length(ivar)];
+i1 = 1:length(ivar);
s2 = diag(gamma_y{1});
sd = sqrt(s2);
-
z = [ m sd s2 ];
-mean = m;
-var = gamma_y{1};
-
-
-%'THEORETICAL MOMENTS';
-%'MEAN','STD. DEV.','VARIANCE');
-z;
%'VARIANCE DECOMPOSITION (in percent)';
if M_.exo_nbr>1
@@ -69,6 +72,7 @@ else
corr = gamma_y{1}(i1,i1);
end
if options_.ar > 0
+ zz=NaN(length(ivar),options_.ar);
%'COEFFICIENTS OF AUTOCORRELATION';
for i=1:options_.ar
if options_.opt_gsa.useautocorr