From 152991864dfb8dcd24343d63a31a0ded5745d8cc Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Tue, 28 Nov 2023 16:53:43 +0100 Subject: [PATCH] GSA: cleanup and removal of globals in filt_mc_.m --- matlab/dynare_sensitivity.m | 2 +- matlab/gsa/filt_mc_.m | 153 +++++++++++++++++------------------- 2 files changed, 73 insertions(+), 82 deletions(-) diff --git a/matlab/dynare_sensitivity.m b/matlab/dynare_sensitivity.m index ed9af2e1d..e8297e973 100644 --- a/matlab/dynare_sensitivity.m +++ b/matlab/dynare_sensitivity.m @@ -434,7 +434,7 @@ if options_gsa.rmse end clear a; % filt_mc_(OutputDirectoryName,data_info); - filt_mc_(OutputDirectoryName,options_gsa,dataset_,dataset_info); + filt_mc_(OutputDirectoryName,options_gsa,dataset_,dataset_info,M_,oo_,options_,bayestopt_,estim_params_); end options_.opt_gsa = options_gsa; options_.prior_trunc=original_prior_trunc; diff --git a/matlab/gsa/filt_mc_.m b/matlab/gsa/filt_mc_.m index ebe9ca9ee..26ceb689d 100644 --- a/matlab/gsa/filt_mc_.m +++ b/matlab/gsa/filt_mc_.m @@ -1,5 +1,21 @@ -function [rmse_MC, ixx] = filt_mc_(OutDir,options_gsa_,dataset_,dataset_info) -% function [rmse_MC, ixx] = filt_mc_(OutDir) +function [rmse_MC, ixx] = filt_mc_(OutDir,options_gsa_,dataset_,dataset_info,M_,oo_,options_,bayestopt_,estim_params_) +% function [rmse_MC, ixx] = filt_mc_(OutDir,options_gsa_,dataset_,dataset_info,M_,oo_,options_,bayestopt_,estim_params_ +% Inputs: +% - OutputDirectoryName [string] name of the output directory +% - options_gsa_ [structure] GSA options +% - dataset_ [dseries] object storing the dataset +% - dataset_info [structure] storing informations about the sample. +% - M_ [structure] Matlab's structure describing the model +% - oo_ [structure] storing the results +% - options_ [structure] Matlab's structure describing the current options +% - bayestopt_ [structure] describing the priors +% - estim_params_ [structure] characterizing parameters to be estimated +% +% Outputs: +% - rmse_MC [double] RMSE by nvar matrix of the RMSEs +% - ixx [double] RMSE by nvar matrix of sorting +% indices (descending order of RMSEs) + % inputs (from opt_gsa structure) % vvarvecm = options_gsa_.var_rmse; % loadSA = options_gsa_.load_rmse; @@ -7,7 +23,6 @@ function [rmse_MC, ixx] = filt_mc_(OutDir,options_gsa_,dataset_,dataset_info) % alpha = options_gsa_.alpha_rmse; % alpha2 = options_gsa_.alpha2_rmse; % istart = options_gsa_.istart_rmse; -% alphaPC = 0.5; % % Written by Marco Ratto % Joint Research Centre, The European Commission, @@ -31,9 +46,6 @@ function [rmse_MC, ixx] = filt_mc_(OutDir,options_gsa_,dataset_,dataset_info) % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . -global bayestopt_ estim_params_ M_ options_ oo_ - -% options_gsa_=options_.opt_gsa; vvarvecm = options_gsa_.var_rmse; if options_.TeX vvarvecm_tex = options_gsa_.var_rmse_tex; @@ -46,11 +58,8 @@ alpha = options_gsa_.alpha_rmse; alpha2 = 0; pvalue = options_gsa_.alpha2_rmse; istart = max(2,options_gsa_.istart_rmse); -alphaPC = 0.5; fname_ = M_.fname; -lgy_ = M_.endo_names; -dr_ = oo_.dr; skipline(2) disp('Starting sensitivity analysis') @@ -61,12 +70,12 @@ if ~options_.nograph a=dir([OutDir,filesep,'*.*']); tmp1='0'; if options_.opt_gsa.ppost - tmp=['_rmse_post']; + tmp='_rmse_post'; else if options_.opt_gsa.pprior - tmp=['_rmse_prior']; + tmp='_rmse_prior'; else - tmp=['_rmse_mc']; + tmp='_rmse_mc'; end if options_gsa_.lik_only tmp1 = [tmp,'_post_SA']; @@ -95,7 +104,7 @@ if options_.opt_gsa.ppost c=load([M_.dname filesep 'Output' filesep fname_,'_mean.mat'],'xparam1'); xparam1_mean=c.xparam1; clear c -elseif ~isempty(options_.mode_file) && exist([M_.dname filesep 'Output' filesep fname_,'_mean.mat'])==2 +elseif ~isempty(options_.mode_file) && exist([M_.dname filesep 'Output' filesep fname_,'_mean.mat'],'file')==2 c=load([M_.dname filesep 'Output' filesep fname_,'_mean.mat'],'xparam1'); xparam1_mean=c.xparam1; clear c @@ -124,31 +133,11 @@ if loadSA end end if ~loadSA - if exist('xparam1','var') - M_ = set_all_parameters(xparam1,estim_params_,M_); - ys_mode=evaluate_steady_state(oo_.steady_state,[oo_.exo_steady_state; oo_.exo_det_steady_state],M_,options_,~options_.steadystate.nocheck); - end - if exist('xparam1_mean','var') - M_ = set_all_parameters(xparam1_mean,estim_params_,M_); - ys_mean=evaluate_steady_state(oo_.steady_state,[oo_.exo_steady_state; oo_.exo_det_steady_state],M_,options_,~options_.steadystate.nocheck); - end Y = transpose(dataset_.data); gend = dataset_.nobs; data_index = dataset_info.missing.aindex; missing_value = dataset_info.missing.state; - for jx=1:gend - data_indx(jx,data_index{jx})=true; - end - load([DirectoryName filesep M_.fname '_data.mat']); filfilt = dir([DirectoryName filesep M_.fname '_filter_step_ahead*.mat']); - temp_smooth_file_list = dir([DirectoryName filesep M_.fname '_smooth*.mat']); - jfile=0; - for j=1:length(temp_smooth_file_list) - if isempty(strfind(temp_smooth_file_list(j).name,'smoothed')), - jfile=jfile+1; - filsmooth(jfile)=temp_smooth_file_list(j); - end - end filupdate = dir([DirectoryName filesep M_.fname '_update*.mat']); filparam = dir([DirectoryName filesep M_.fname '_param*.mat']); x=[]; @@ -156,11 +145,11 @@ if ~loadSA sto_ys=[]; for j=1:length(filparam) if isempty(strmatch([M_.fname '_param_irf'],filparam(j).name)) - load([DirectoryName filesep filparam(j).name]); - x=[x; stock]; - logpo2=[logpo2; stock_logpo]; - sto_ys=[sto_ys; stock_ys]; - clear stock stock_logpo stock_ys; + temp=load([DirectoryName filesep filparam(j).name]); + x=[x; temp.stock]; + logpo2=[logpo2; temp.stock_logpo]; + sto_ys=[sto_ys; temp.stock_ys]; + clear temp; end end nruns=size(x,1); @@ -168,38 +157,41 @@ if ~loadSA if options_.opt_gsa.ppost || (options_.opt_gsa.ppost==0 && options_.opt_gsa.lik_only==0) skipline() disp('Computing RMSE''s...') + jxj=NaN(length(vvarvecm),1); + js=NaN(length(vvarvecm),1); + yss=NaN(length(vvarvecm),gend,size(sto_ys,1)); for i = 1:length(vvarvecm) vj = vvarvecm{i}; - jxj(i) = strmatch(vj, lgy_(dr_.order_var), 'exact'); - js(i) = strmatch(vj, lgy_, 'exact'); + jxj(i) = strmatch(vj, M_.endo_names(oo_.dr.order_var), 'exact'); + js(i) = strmatch(vj, M_.endo_names, 'exact'); yss(i,:,:)=repmat(sto_ys(:,js(i))',[gend,1]); end if exist('xparam1','var') - [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK] = DsgeSmoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_); - y0 = reshape( squeeze(aK(1,jxj,1:gend)),[gend length(jxj)]);% + kron(ys_mode(js),ones(1,gend))); - yobs = transpose( ahat(jxj,:));% + kron(ys_mode(js),ones(1,gend))); + [~,~,~,ahat,~,~,aK] = DsgeSmoother(xparam1,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_); + y0 = reshape( squeeze(aK(1,jxj,1:gend)),[gend length(jxj)]); + yobs = transpose( ahat(jxj,:)); rmse_mode = sqrt(mean((yobs(istart:end,:)-y0(istart:end,:)).^2)); r2_mode = 1-sum((yobs(istart:end,:)-y0(istart:end,:)).^2)./sum(yobs(istart:end,:).^2); end + y0=-yss; nbb=0; for j=1:length(filfilt) - load([DirectoryName filesep M_.fname '_filter_step_ahead',num2str(j),'.mat']); - nb = size(stock,4); - y0(:,:,nbb+1:nbb+nb)=y0(:,:,nbb+1:nbb+nb)+reshape(stock(1,js,1:gend,:),[length(js) gend nb]); + temp=load([DirectoryName filesep M_.fname '_filter_step_ahead',num2str(j),'.mat']); + nb = size(temp.stock,4); + y0(:,:,nbb+1:nbb+nb)=y0(:,:,nbb+1:nbb+nb)+reshape(temp.stock(1,js,1:gend,:),[length(js) gend nb]); nbb=nbb+nb; - clear stock; + clear temp; end yobs=-yss; nbb=0; for j=1:length(filupdate) - load([DirectoryName filesep M_.fname '_update',num2str(j),'.mat']); - nb = size(stock,3); - yobs(:,:,nbb+1:nbb+nb)=yobs(:,:,nbb+1:nbb+nb)+reshape(stock(js,1:gend,:),[length(js) gend nb]); + temp=load([DirectoryName filesep M_.fname '_update',num2str(j),'.mat']); + nb = size(temp.stock,3); + yobs(:,:,nbb+1:nbb+nb)=yobs(:,:,nbb+1:nbb+nb)+reshape(temp.stock(js,1:gend,:),[length(js) gend nb]); nbb=nbb+nb; - clear stock; + clear temp; end - y0M=mean(y0,2); rmse_MC=zeros(nruns,length(js)); r2_MC=zeros(nruns,length(js)); for j=1:nruns @@ -207,14 +199,15 @@ if ~loadSA r2_MC(j,:) = 1-mean((yobs(:,istart:end,j)'-y0(:,istart:end,j)').^2)./mean((yobs(:,istart:end,j)').^2); end if exist('xparam1_mean','var') - [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK] = DsgeSmoother(xparam1_mean,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_); - y0 = reshape( squeeze(aK(1,jxj,1:gend)),[gend length(jxj)]);% + kron(ys_mean(js),ones(1,gend))); - yobs = transpose( ahat(jxj,:));% + kron(ys_mean(js),ones(1,gend))); + [~,~,~,ahat,~,~,aK] = DsgeSmoother(xparam1_mean,gend,Y,data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_); + y0 = reshape( squeeze(aK(1,jxj,1:gend)),[gend length(jxj)]); + yobs = transpose( ahat(jxj,:)); rmse_pmean = sqrt(mean((yobs(istart:end,:)-y0(istart:end,:)).^2)); r2_pmean = 1-mean((yobs(istart:end,:)-y0(istart:end,:)).^2)./mean(yobs(istart:end,:).^2); end clear stock_filter; end + lnprior=NaN(nruns,1); for j=1:nruns lnprior(j,1) = priordens(x(j,:)',bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7,bayestopt_.p3,bayestopt_.p4); end @@ -242,7 +235,7 @@ if ~loadSA end end end -else +else % loadSA if options_.opt_gsa.lik_only && options_.opt_gsa.ppost==0 load([OutDir,filesep,fnamtmp, '.mat'],'x','logpo2','likelihood'); else @@ -252,12 +245,12 @@ else nruns=size(x,1); nfilt=floor(pfilt*nruns); end -% smirnov tests +% Smirnov tests nfilt0 = nfilt*ones(length(vvarvecm), 1); logpo2=logpo2(:); if ~options_.opt_gsa.ppost - [dum, ipost]=sort(-logpo2); - [dum, ilik]=sort(-likelihood); + [~, ipost]=sort(-logpo2); + [~, ilik]=sort(-likelihood); end % visual scatter analysis! @@ -333,8 +326,8 @@ else rmse_txt=rmse_pmean; r2_txt=r2_pmean; else - if options_.opt_gsa.pprior || ~exist('rmse_pmean') - if exist('rmse_mode') + if options_.opt_gsa.pprior || ~exist('rmse_pmean','var') + if exist('rmse_mode','var') rmse_txt=rmse_mode; r2_txt=r2_mode; else @@ -346,18 +339,19 @@ else r2_txt=r2_pmean; end end + ixx=NaN(size(rmse_MC,1),length(vvarvecm)); for i = 1:length(vvarvecm) - [dum, ixx(:,i)] = sort(rmse_MC(:,i)); + [~, ixx(:,i)] = sort(rmse_MC(:,i)); end PP = ones(npar+nshock, length(vvarvecm)); PPV = ones(length(vvarvecm), length(vvarvecm), npar+nshock); SS = zeros(npar+nshock, length(vvarvecm)); for j = 1:npar+nshock for i = 1:length(vvarvecm) - [H, P, KSSTAT] = smirnov(x(ixx(nfilt0(i)+1:end,i),j),x(ixx(1:nfilt0(i),i),j), alpha); - [H1, P1, KSSTAT1] = smirnov(x(ixx(nfilt0(i)+1:end,i),j),x(ixx(1:nfilt0(i),i),j),alpha,1); - [H2, P2, KSSTAT2] = smirnov(x(ixx(nfilt0(i)+1:end,i),j),x(ixx(1:nfilt0(i),i),j),alpha,-1); - if H1 & H2==0 + [~, P] = smirnov(x(ixx(nfilt0(i)+1:end,i),j),x(ixx(1:nfilt0(i),i),j), alpha); + [H1] = smirnov(x(ixx(nfilt0(i)+1:end,i),j),x(ixx(1:nfilt0(i),i),j),alpha,1); + [H2] = smirnov(x(ixx(nfilt0(i)+1:end,i),j),x(ixx(1:nfilt0(i),i),j),alpha,-1); + if H1==0 && H2==0 SS(j,i)=1; elseif H1==0 SS(j,i)=-1; @@ -369,7 +363,7 @@ else for i = 1:length(vvarvecm) for l = 1:length(vvarvecm) if l~=i && PP(j,i)1)); - snam=param_names(find(nsp>0)); + snam0=param_names(nsp==0); + snam1=param_names(nsp==1); + snam2=param_names(nsp>1); nsnam=(find(nsp>1)); skipline(2) disp('These parameters do not affect significantly the fit of ANY observed series:') @@ -767,7 +758,7 @@ else set(h0,'color',a00(i,:),'linewidth',2) end ydum=get(gca,'ylim'); - if exist('xparam1') + if exist('xparam1','var') xdum=xparam1(ipar(j)); h1=plot([xdum xdum],ydum); set(h1,'color',[0.85 0.85 0.85],'linewidth',2) @@ -824,7 +815,7 @@ else set(h0,'color',a00(i,:),'linewidth',2) end ydum=get(gca,'ylim'); - if exist('xparam1') + if exist('xparam1','var') xdum=xparam1(nsnam(j)); h1=plot([xdum xdum],ydum); set(h1,'color',[0.85 0.85 0.85],'linewidth',2)