From 28da4a09c5e36d5f3792345e4ef17528a7a3e928 Mon Sep 17 00:00:00 2001 From: ratto Date: Wed, 4 Apr 2007 08:19:51 +0000 Subject: [PATCH] Updated sensitivity package !!! git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@1235 ac1d8469-bf42-47a9-8791-bf33cf982152 --- matlab/dat_fil_.m | 26 +++ matlab/dynare_MC.m | 249 +++++----------------------- matlab/dynare_sensitivity.m | 200 ++++++++++++++++++----- matlab/filt_mc_.m | 269 +++++++++++++++++------------- matlab/ghx2transition.m | 41 +++++ matlab/gsa_sdp_dyn.m | 74 +++++++++ matlab/prior_draw_gsa.m | 125 ++++++++++++++ matlab/read_data.m | 39 +++++ matlab/redform_map.m | 314 ++++++++++++++++++++++++++++++++++++ matlab/redform_screen.m | 166 +++++++++++++++++++ matlab/smirnov.m | 15 +- matlab/stab_map_.m | 273 ++++++++++++++++++++----------- matlab/stab_map_1.m | 20 ++- matlab/stab_map_2.m | 23 +-- matlab/teff.m | 34 ++++ matlab/trank.m | 25 +++ 16 files changed, 1426 insertions(+), 467 deletions(-) create mode 100644 matlab/dat_fil_.m create mode 100644 matlab/ghx2transition.m create mode 100644 matlab/gsa_sdp_dyn.m create mode 100644 matlab/prior_draw_gsa.m create mode 100644 matlab/read_data.m create mode 100644 matlab/redform_map.m create mode 100644 matlab/redform_screen.m create mode 100644 matlab/teff.m create mode 100644 matlab/trank.m diff --git a/matlab/dat_fil_.m b/matlab/dat_fil_.m new file mode 100644 index 000000000..04c3cbdcb --- /dev/null +++ b/matlab/dat_fil_.m @@ -0,0 +1,26 @@ +function c = dat_fil_(data_file); +% +% Part of the Sensitivity Analysis Toolbox for DYNARE +% +% Written by Marco Ratto, 2006 +% Joint Research Centre, The European Commission, +% (http://eemc.jrc.ec.europa.eu/), +% marco.ratto@jrc.it +% +% Disclaimer: This software is not subject to copyright protection and is in the public domain. +% It is an experimental system. The Joint Research Centre of European Commission +% assumes no responsibility whatsoever for its use by other parties +% and makes no guarantees, expressed or implied, about its quality, reliability, or any other +% characteristic. We would appreciate acknowledgement if the software is used. +% Reference: +% M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006. +% + +eval(data_file); +clear data_file; + +a=who; + +for j=1:length(a) + eval(['c.',a{j},'=',a{j},';']); +end \ No newline at end of file diff --git a/matlab/dynare_MC.m b/matlab/dynare_MC.m index ac59ae645..d947ab1f2 100644 --- a/matlab/dynare_MC.m +++ b/matlab/dynare_MC.m @@ -1,179 +1,44 @@ function dynare_MC(var_list_,OutDir) +% +% Adapted by M. Ratto from dynare_estimation.m and posteriorsmoother.m +% (dynare_estimation.m and posteriorsmoother.m are part of DYNARE, +% copyright M. Juillard) +% +% Part of the Sensitivity Analysis Toolbox for DYNARE +% +% Written by Marco Ratto, 2006 +% Joint Research Centre, The European Commission, +% (http://eemc.jrc.ec.europa.eu/), +% marco.ratto@jrc.it +% +% Disclaimer: This software is not subject to copyright protection and is in the public domain. +% It is an experimental system. The Joint Research Centre of European Commission +% assumes no responsibility whatsoever for its use by other parties +% and makes no guarantees, expressed or implied, about its quality, reliability, or any other +% characteristic. We would appreciate acknowledgement if the software is used. +% Reference: +% M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006. +% global M_ options_ oo_ estim_params_ global bayestopt_ -% temporary fix until M_.H is initialized by the parser -M_.H = []; - -options_.varlist = var_list_; -options_.lgyidx2varobs = zeros(size(M_.endo_names,1),1); -for i = 1:size(M_.endo_names,1) - tmp = strmatch(deblank(M_.endo_names(i,:)),options_.varobs,'exact'); - if ~isempty(tmp) - options_.lgyidx2varobs(i,1) = tmp; - end -end - -options_ = set_default_option(options_,'first_obs',1); -options_ = set_default_option(options_,'prefilter',0); -options_ = set_default_option(options_,'presample',0); -options_ = set_default_option(options_,'lik_algo',1); -options_ = set_default_option(options_,'lik_init',1); -options_ = set_default_option(options_,'nograph',0); -options_ = set_default_option(options_,'mh_conf_sig',0.90); -options_ = set_default_option(options_,'mh_replic',20000); -options_ = set_default_option(options_,'mh_drop',0.5); -options_ = set_default_option(options_,'mh_jscale',0.2); -options_ = set_default_option(options_,'mh_init_scale',2*options_.mh_jscale); -options_ = set_default_option(options_,'mode_file',''); -options_ = set_default_option(options_,'mode_compute',4); -options_ = set_default_option(options_,'mode_check',0); -options_ = set_default_option(options_,'prior_trunc',1e-10); -options_ = set_default_option(options_,'mh_mode',1); -options_ = set_default_option(options_,'mh_nblck',2); -options_ = set_default_option(options_,'load_mh_file',0); -options_ = set_default_option(options_,'nodiagnostic',0); -options_ = set_default_option(options_,'loglinear',0); -options_ = set_default_option(options_,'unit_root_vars',[]); -options_ = set_default_option(options_,'XTick',[]); -options_ = set_default_option(options_,'XTickLabel',[]); -options_ = set_default_option(options_,'bayesian_irf',0); -options_ = set_default_option(options_,'bayesian_th_moments',0); -options_ = set_default_option(options_,'TeX',0); -options_ = set_default_option(options_,'irf',40); -options_ = set_default_option(options_,'relative_irf',0); -options_ = set_default_option(options_,'order',1); -options_ = set_default_option(options_,'ar',5); -options_ = set_default_option(options_,'dr_algo',0); -options_ = set_default_option(options_,'linear',0); -options_ = set_default_option(options_,'drop',0); -options_ = set_default_option(options_,'replic',1); -options_ = set_default_option(options_,'hp_filter',0); -options_ = set_default_option(options_,'forecast',0); -options_ = set_default_option(options_,'smoother',0); -options_ = set_default_option(options_,'moments_varendo',0); -options_ = set_default_option(options_,'filtered_vars',0); -options_ = set_default_option(options_,'kalman_algo',1); -options_ = set_default_option(options_,'kalman_tol',10^(-12)); -options_ = set_default_option(options_,'posterior_mode_estimation',1); -options_ = set_default_option(options_,'MaxNumberOfBytes',1e6); -options_ = set_default_option(options_,'xls_sheet',''); -options_ = set_default_option(options_,'xls_range',''); -options_ = set_default_option(options_,'filter_step_ahead',0); -options_ = set_default_option(options_,'diffuse_d',[]); -options_ = set_default_option(options_,'Opt6Iter',3); -options_ = set_default_option(options_,'Opt6Numb',100000); -options_ = set_default_option(options_,'steadystate_flag',0); -options_ = set_default_option(options_,'logdata',0); -options_ = set_default_option(options_,'use_mh_covariance_matrix',0); - -if exist([M_.fname '_steadystate.m']) - options_.steadystate_flag = 1; -end - if options_.filtered_vars ~= 0 & options_.filter_step_ahead == 0 options_.filter_step_ahead = 1; end -if options_.bayesian_irf & options_.irf == 0 - options_.irf = 40; -end if options_.filter_step_ahead ~= 0 options_.nk = max(options_.filter_step_ahead); else options_.nk = 0; end -%% Add something to the parser ++> -M_.dname = M_.fname; % The user should be able to choose another name - % for the directory... - - -pnames = [' ';'beta ';'gamm ';'norm ';'invg ';'unif ';'invg2']; -n_varobs = size(options_.varobs,1); - -[xparam1,estim_params_,bayestopt_,lb,ub] = set_prior(estim_params_); - -bounds = prior_bounds(bayestopt_); -bounds(:,1)=max(bounds(:,1),lb); -bounds(:,2)=min(bounds(:,2),ub); - -if any(xparam1 < bounds(:,1)) | any(xparam1 > bounds(:,2)) - find(xparam1 < bounds(:,1)) - find(xparam1 > bounds(:,2)) - error('Initial parameter values are outside parameter bounds') -end -lb = bounds(:,1); -ub = bounds(:,2); -bayestopt_.lb = lb; -bayestopt_.ub = ub; - -if ~isfield(options_,'trend_coeffs') - bayestopt_.with_trend = 0; -else - bayestopt_.with_trend = 1; - bayestopt_.trend_coeff = {}; - trend_coeffs = options_.trend_coeffs; - nt = length(trend_coeffs); - for i=1:n_varobs - if i > length(trend_coeffs) - bayestopt_.trend_coeff{i} = '0'; - else - bayestopt_.trend_coeff{i} = trend_coeffs{i}; - end - end -end - -bayestopt_.penalty = 1e8; % penalty nvx = estim_params_.nvx; nvn = estim_params_.nvn; ncx = estim_params_.ncx; ncn = estim_params_.ncn; np = estim_params_.np ; -nx = nvx+nvn+ncx+ncn+np; - -dr = set_state_space([]); - -%% Initialization with unit-root variables -if ~isempty(options_.unit_root_vars) - n_ur = size(options_.unit_root_vars,1); - i_ur = zeros(n_ur,1); - for i=1:n_ur - i1 = strmatch(deblank(options_.unit_root_vars(i,:)),M_.endo_names(dr.order_var,:),'exact'); - if isempty(i1) - error('Undeclared variable in unit_root_vars statement') - end - i_ur(i) = i1; - end - if M_.maximum_lag > 1 - l1 = flipud([cumsum(M_.lead_lag_incidence(1:M_.maximum_lag-1,dr.order_var),1);ones(1,M_.endo_nbr)]); - n1 = nnz(l1); - bayestopt_.Pinf = zeros(n1,n1); - l2 = find(l1'); - l3 = zeros(M_.endo_nbr,M_.maximum_lag); - l3(i_ur,:) = l1(:,i_ur)'; - l3 = l3(:); - i_ur1 = find(l3(l2)); - i_stable = ones(M_.endo_nbr,1); - i_stable(i_ur) = zeros(n_ur,1); - i_stable = find(i_stable); - bayestopt_.Pinf(i_ur1,i_ur1) = diag(ones(1,length(i_ur1))); - bayestopt_.i_var_stable = i_stable; - l3 = zeros(M_.endo_nbr,M_.maximum_lag); - l3(i_stable,:) = l1(:,i_stable)'; - l3 = l3(:); - bayestopt_.i_T_var_stable = find(l3(l2)); - else - n1 = M_.endo_nbr; - bayestopt_.Pinf = zeros(n1,n1); - bayestopt_.Pinf(i_ur,i_ur) = diag(ones(1,length(i_ur))); - l1 = ones(M_.endo_nbr,1); - l1(i_ur,:) = zeros(length(i_ur),1); - bayestopt_.i_T_var_stable = find(l1); - end - options_.lik_init = 3; -end % if ~isempty(options_.unit_root_vars) +npar = nvx+nvn+ncx+ncn+np; if isempty(options_.datafile) error('ESTIMATION: datafile option is missing') @@ -183,36 +48,9 @@ if isempty(options_.varobs) error('ESTIMATION: VAROBS is missing') end - -%% If jscale isn't specified for an estimated parameter, use -%% global option options_.jscale, set to 0.2, by default -k = find(isnan(bayestopt_.jscale)); -bayestopt_.jscale(k) = options_.mh_jscale; - %% Read and demean data rawdata = read_variables(options_.datafile,options_.varobs,[],options_.xls_sheet,options_.xls_range); -k = []; -k1 = []; -for i=1:n_varobs - k = [k strmatch(deblank(options_.varobs(i,:)),M_.endo_names(dr.order_var,:),'exact')]; - k1 = [k1 strmatch(deblank(options_.varobs(i,:)),M_.endo_names, 'exact')]; -end -% union of observed and state variables -k2 = union(k',[dr.nstatic+1:dr.nstatic+dr.npred]'); -% including variables in t-2 and earlier, if any -k2 = [k2;[M_.endo_nbr+(1:dr.nspred-dr.npred)]']; - -% set restrict_state to postion of observed + state variables -% in expanded state vector -bayestopt_.restrict_state = k2; -% set mf1 to positions of observed variables in restricted state vector -% for likelihood computation -[junk,bayestopt_.mf1] = ismember(k,k2); -% set mf2 to positions of observed variables in expanded state vector -% for filtering and smoothing -bayestopt_.mf2 = k; -bayestopt_.mfys = k1; options_ = set_default_option(options_,'nobs',size(rawdata,1)-options_.first_obs+1); gend = options_.nobs; @@ -232,12 +70,6 @@ if ~isreal(rawdata) ' transformation']) end -nvx = estim_params_.nvx; -nvn = estim_params_.nvn; -ncx = estim_params_.ncx; -ncn = estim_params_.ncn; -np = estim_params_.np ; -npar = nvx+nvn+ncx+ncn+np; offset = npar-np; fname_=M_.fname; @@ -257,39 +89,40 @@ load(options_.mode_file) x=[lpmat0(istable,:) lpmat(istable,:)]; clear lpmat lpmat0 istable %iunstable egg yys T B = size(x,1); -if M_.maximum_lag > 1 - l1 = flipud([cumsum(M_.lead_lag_incidence(1:M_.maximum_lag-1,dr.order_var),1);ones(1,M_.endo_nbr)]); - n1 = nnz(l1); -else - n1 = M_.endo_nbr; -end +[atT,innov,measurement_error,filtered_state_vector,ys,trend_coeff, aK] = DsgeSmoother(x(1,:)',gend,data); +n1=size(atT,1); + nfil=B/40; -stock_smooth = zeros(gend,n1,40); -stock_filter = zeros(gend+1,n1,40); -stock_ys = zeros(M_.endo_nbr,40); +stock_smooth = zeros(M_.endo_nbr,gend,40); +stock_filter = zeros(M_.endo_nbr,gend+1,40); +stock_ys = zeros(40, M_.endo_nbr); logpo2=zeros(B,1); %% h = waitbar(0,'MC smoother ...'); +delete([OutDir,'\',namfile,'_*.mat']) ib=0; ifil=0; opt_gsa=options_.opt_gsa; for b=1:B ib=ib+1; - deep = x(b,:); + deep = x(b,:)'; + set_all_parameters(deep); + dr = resol(oo_.steady_state,0); %deep(1:offset) = xparam1(1:offset); - logpo2(b,1) = DsgeLikelihood(deep',gend,data); + logpo2(b,1) = DsgeLikelihood(deep,gend,data); if opt_gsa.lik_only==0, - [atT,innov,measurement_error,filtered_state_vector,ys,trend_coeff] = DsgeSmoother(deep,gend,data); - stock_smooth(:,:,ib)=atT'; - stock_filter(:,:,ib)=filtered_state_vector'; - stock_ys(:,ib)=ys; + [atT,innov,measurement_error,filtered_state_vector,ys,trend_coeff, aK] = DsgeSmoother(deep,gend,data); + stock_smooth(:,:,ib)=atT(1:M_.endo_nbr,:); + stock_filter(:,:,ib)=filtered_state_vector(1:M_.endo_nbr,:); + %stock_filter(:,:,ib)=aK(options_.filter_step_ahead,1:M_.endo_nbr,:); + stock_ys(ib,:)=ys'; if ib==40, ib=0; ifil=ifil+1; save([OutDir,'\',namfile,'_',num2str(ifil)],'stock_smooth','stock_filter','stock_ys') - stock_smooth = zeros(gend,n1,40); - stock_filter = zeros(gend+1,n1,40); - stock_ys = zeros(M_.endo_nbr,40); + stock_smooth = zeros(M_.endo_nbr,gend,40); + stock_filter = zeros(M_.endo_nbr,gend+1,40); + stock_ys = zeros(40, M_.endo_nbr); end end waitbar(b/B,h,['MC smoother ...',num2str(b),'/',num2str(B)]); @@ -300,7 +133,7 @@ if ib>0, ifil=ifil+1; stock_smooth = stock_smooth(:,:,1:ib); stock_filter = stock_filter(:,:,1:ib); - stock_ys = stock_ys(:,1:ib); + stock_ys = stock_ys(1:ib,:); save([OutDir,'\',namfile,'_',num2str(ifil)],'stock_smooth','stock_filter','stock_ys') end end diff --git a/matlab/dynare_sensitivity.m b/matlab/dynare_sensitivity.m index 40bf4fee6..df0c1e895 100644 --- a/matlab/dynare_sensitivity.m +++ b/matlab/dynare_sensitivity.m @@ -1,5 +1,19 @@ function x0=dynare_sensitivity() -% copyright Marco Ratto 2006 +% Part of the Sensitivity Analysis Toolbox for DYNARE +% +% Written by Marco Ratto, 2006 +% Joint Research Centre, The European Commission, +% (http://eemc.jrc.ec.europa.eu/), +% marco.ratto@jrc.it +% +% Disclaimer: This software is not subject to copyright protection and is in the public domain. +% It is an experimental system. The Joint Research Centre of European Commission +% assumes no responsibility whatsoever for its use by other parties +% and makes no guarantees, expressed or implied, about its quality, reliability, or any other +% characteristic. We would appreciate acknowledgement if the software is used. +% Reference: +% M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006. +% global M_ options_ oo_ bayestopt_ @@ -10,29 +24,100 @@ x0=[]; options_ = set_default_option(options_,'opt_gsa',1); options_gsa_ = options_.opt_gsa; +options_gsa_ = set_default_option(options_gsa_,'identification',0); +if options_gsa_.identification, + options_gsa_.redform=0; + options_gsa_ = set_default_option(options_gsa_,'morris',1); +end + % map stability options_gsa_ = set_default_option(options_gsa_,'stab',1); options_gsa_ = set_default_option(options_gsa_,'redform',0); options_gsa_ = set_default_option(options_gsa_,'pprior',1); +options_gsa_ = set_default_option(options_gsa_,'prior_range',1); options_gsa_ = set_default_option(options_gsa_,'ppost',0); options_gsa_ = set_default_option(options_gsa_,'ilptau',1); +options_gsa_ = set_default_option(options_gsa_,'morris',0); +options_gsa_ = set_default_option(options_gsa_,'glue',0); +options_gsa_ = set_default_option(options_gsa_,'morris_nliv',6); +options_gsa_ = set_default_option(options_gsa_,'morris_ntra',20); options_gsa_ = set_default_option(options_gsa_,'Nsam',2048); +options_gsa_ = set_default_option(options_gsa_,'load_redform',0); +options_gsa_ = set_default_option(options_gsa_,'load_rmse',0); options_gsa_ = set_default_option(options_gsa_,'load_stab',0); options_gsa_ = set_default_option(options_gsa_,'alpha2_stab',0.3); options_gsa_ = set_default_option(options_gsa_,'ksstat',0.1); %options_gsa_ = set_default_option(options_gsa_,'load_mh',0); -OutputDirectoryName = CheckPath('GSA'); +if options_gsa_.redform, + options_gsa_.pprior=1; + options_gsa_.ppost=0; +end +if options_gsa_.morris, + if ~(exist('Sampling_Function_2','file')==6 | exist('Sampling_Function_2','file')==2), + disp('Download pre-parsed mapping routines at:') + disp('http://eemc.jrc.ec.europa.eu/softwareDYNARE-Dowload.htm') + disp(' ' ) + error('Mapping routines missing!') + + end + + options_gsa_.pprior=1; + options_gsa_.ppost=0; + %options_gsa_.stab=1; + options_gsa_.glue=0; + options_gsa_.rmse=0; + options_gsa_.load_rmse=0; + options_gsa_.alpha2_stab=1; + options_.ksstat=1; + if options_gsa_.morris==2, + options_gsa_ = set_default_option(options_gsa_,'Nsam',256); + OutputDirectoryName = CheckPath('GSA/IDENTIF'); + else + OutputDirectoryName = CheckPath('GSA/SCREEN'); + end +else + OutputDirectoryName = CheckPath('GSA'); +end + +options_.opt_gsa = options_gsa_; + +if (options_gsa_.load_stab | options_gsa_.load_rmse | options_gsa_.load_redform) ... + & options_gsa_.pprior, + filetoload=[OutputDirectoryName '\' fname_ '_prior.mat']; + if isempty(ls(filetoload)), + disp([filetoload,' not found!']) + disp(['You asked to load a non existent analysis']) + %options_gsa_.load_stab=0; + return, + else + if isempty(strmatch('bkpprior',who('-file', filetoload),'exact')), + disp('Warning! Missing prior info for saved sample') % trap for files previous + disp('The saved files are generated with previous version of GSA package') % trap for files previous + else + load(filetoload,'bkpprior'), + if any(bayestopt_.pshape~=bkpprior.pshape) | ... + any(bayestopt_.p1~=bkpprior.p1) | ... + any(bayestopt_.p2~=bkpprior.p2) | ... + any(bayestopt_.p3(~isnan(bayestopt_.p3))~=bkpprior.p3(~isnan(bkpprior.p3))) | ... + any(bayestopt_.p4(~isnan(bayestopt_.p4))~=bkpprior.p4(~isnan(bkpprior.p4))), + disp('WARNING!') + disp('The saved sample has different priors w.r.t. to current ones!!') + disp('') + disp('Press ENTER to continue') + pause; + end + end + end +end if options_gsa_.stab & ~options_gsa_.ppost, - x0 = stab_map_(options_gsa_.Nsam, options_gsa_.load_stab, options_gsa_.ksstat, options_gsa_.alpha2_stab, ... - options_gsa_.redform, options_gsa_.pprior, options_gsa_.ilptau, OutputDirectoryName); + x0 = stab_map_(OutputDirectoryName); end % reduced form % redform_map(namendo, namlagendo, namexo, icomp, pprior, ilog, threshold) -options_gsa_ = set_default_option(options_gsa_,'load_redform',0); options_gsa_ = set_default_option(options_gsa_,'logtrans_redform',0); options_gsa_ = set_default_option(options_gsa_,'threshold_redform',[]); options_gsa_ = set_default_option(options_gsa_,'ksstat_redform',0.1); @@ -41,20 +126,22 @@ options_gsa_ = set_default_option(options_gsa_,'namendo',[]); options_gsa_ = set_default_option(options_gsa_,'namlagendo',[]); options_gsa_ = set_default_option(options_gsa_,'namexo',[]); -if options_gsa_.redform & ~isempty(options_gsa_.namendo) & ~options_gsa_.ppost, - redform_map(options_gsa_.namendo, options_gsa_.namlagendo, options_gsa_.namexo, ... - options_gsa_.load_redform, options_gsa_.pprior, options_gsa_.logtrans_redform, ... - options_gsa_.threshold_redform, options_gsa_.ksstat_redform, ... - options_gsa_.alpha2_redform, OutputDirectoryName); +options_.opt_gsa = options_gsa_; +if options_gsa_.redform & ~isempty(options_gsa_.namendo) ... + & ~options_gsa_.ppost, + if options_gsa_.morris, + redform_screen(OutputDirectoryName); + else + redform_map(OutputDirectoryName); + end end % RMSE mapping % function [rmse_MC, ixx] = filt_mc_(vvarvecm, loadSA, pfilt, alpha, alpha2) options_gsa_ = set_default_option(options_gsa_,'rmse',0); options_gsa_ = set_default_option(options_gsa_,'lik_only',0); options_gsa_ = set_default_option(options_gsa_,'var_rmse',options_.varobs); -options_gsa_ = set_default_option(options_gsa_,'load_rmse',0); options_gsa_ = set_default_option(options_gsa_,'pfilt_rmse',0.1); -options_gsa_ = set_default_option(options_gsa_,'istart_rmse',1); +options_gsa_ = set_default_option(options_gsa_,'istart_rmse',options_.presample+1); options_gsa_ = set_default_option(options_gsa_,'alpha_rmse',0.002); options_gsa_ = set_default_option(options_gsa_,'alpha2_rmse',1); options_.opt_gsa = options_gsa_; @@ -68,16 +155,27 @@ if options_gsa_.rmse, if isempty(a), dynare_MC([],OutputDirectoryName); options_gsa_.load_rmse=0; +% else +% if options_gsa_.load_rmse==0, +% disp('You already saved a MC filter/smoother analysis ') +% disp('Do you want to overwrite ?') +% pause; +% if options_gsa_.pprior +% delete([OutputDirectoryName,'\',fname_,'_prior_*.mat']) +% else +% delete([OutputDirectoryName,'\',fname_,'_mc_*.mat']) +% end +% dynare_MC([],OutputDirectoryName); +% options_gsa_.load_rmse=0; +% end + end end clear a; - filt_mc_(options_gsa_.var_rmse, options_gsa_.load_rmse, options_gsa_.pfilt_rmse, ... - options_gsa_.alpha_rmse, options_gsa_.alpha2_rmse, OutputDirectoryName, ... - options_gsa_.istart_rmse); + filt_mc_(OutputDirectoryName); end -options_gsa_ = set_default_option(options_gsa_,'glue',0); if options_gsa_.glue, dr_ = oo_.dr; if options_gsa_.ppost @@ -90,6 +188,11 @@ if options_gsa_.glue, load([OutputDirectoryName,'\',fname_,'_mc']); end end + if ~exist('x'), + disp(['No RMSE analysis is available for current options']) + disp(['No GLUE file prepared']) + return, + end nruns=size(x,1); gend = options_.nobs; rawdata = read_variables(options_.datafile,options_.varobs,[],options_.xls_sheet,options_.xls_range); @@ -114,12 +217,14 @@ if options_gsa_.glue, jxj = strmatch(vj,lgy_(dr_.order_var,:),'exact'); js = strmatch(vj,lgy_,'exact'); if ~options_gsa_.ppost - y0=zeros(gend+1,nruns); - nb = size(stock_filter,3); - y0 = squeeze(stock_filter(:,jxj,:)) + ... - kron(stock_ys(js,:),ones(size(stock_filter,1),1)); - Out(j).data = y0'; - Out(j).time = [1:size(y0,1)]; +% y0=zeros(gend+1,nruns); +% nb = size(stock_filter,3); +% y0 = squeeze(stock_filter(:,jxj,:)) + ... +% kron(stock_ys(js,:),ones(size(stock_filter,1),1)); +% Out(j).data = y0'; +% Out(j).time = [1:size(y0,1)]; + Out(j).data = jxj; + Out(j).time = [pwd,'\',OutputDirectoryName]; else Out(j).data = jxj; Out(j).time = [pwd,'\',DirectoryName]; @@ -132,12 +237,13 @@ if options_gsa_.glue, Lik(j).data = rmse_MC(:,j)'; if ~options_gsa_.ppost - y0 = squeeze(stock_smooth(:,jxj,:)) + ... - kron(stock_ys(js,:),ones(size(stock_smooth,1),1)); - Out1(j).name = vj; - Out1(j).ini = 'yes'; - Out1(j).time = [1:size(y0,1)]; - Out1(j).data = y0'; +% y0 = squeeze(stock_smooth(:,jxj,:)) + ... +% kron(stock_ys(js,:),ones(size(stock_smooth,1),1)); +% Out1(j).name = vj; +% Out1(j).ini = 'yes'; +% Out1(j).time = [1:size(y0,1)]; +% Out1(j).data = y0'; + Out1=Out; else Out1=Out; end @@ -150,10 +256,12 @@ if options_gsa_.glue, jsmoo=jsmoo+1; vj=deblank(M_.endo_names(dr_.order_var(j),:)); if ~options_gsa_.ppost - y0 = squeeze(stock_smooth(:,j,:)) + ... - kron(stock_ys(j,:),ones(size(stock_smooth,1),1)); - Out1(jsmoo).time = [1:size(y0,1)]; - Out1(jsmoo).data = y0'; +% y0 = squeeze(stock_smooth(:,j,:)) + ... +% kron(stock_ys(j,:),ones(size(stock_smooth,1),1)); +% Out1(jsmoo).time = [1:size(y0,1)]; +% Out1(jsmoo).data = y0'; + Out1(jsmoo).data = j; + Out1(jsmoo).time = [pwd,'\',OutputDirectoryName]; else Out1(jsmoo).data = j; Out1(jsmoo).time = [pwd,'\',DirectoryName]; @@ -179,22 +287,30 @@ if options_gsa_.glue, Rem.id = 'Original'; Rem.ind= [1:size(x,1)]; + Info.dynare=M_.fname; + Info.order_var=dr_.order_var; + Out=Out1; if options_gsa_.ppost - Info.dynare=M_.fname; - Info.order_var=dr_.order_var; - Out=Out1; - save([OutputDirectoryName,'\',fname_,'_post_glue'], 'Out', 'Sam', 'Lik', 'Obs', 'Rem','Info', 'Exo') +% Info.dynare=M_.fname; +% Info.order_var=dr_.order_var; +% Out=Out1; + Info.TypeofSample='post'; + save([OutputDirectoryName,'\',fname_,'_glue_post'], 'Out', 'Sam', 'Lik', 'Obs', 'Rem','Info', 'Exo') %save([fname_,'_post_glue_smooth'], 'Out', 'Sam', 'Lik', 'Obs', 'Rem','Info') else if options_gsa_.pprior - save([OutputDirectoryName,'\',fname_,'_prior_glue'], 'Out', 'Sam', 'Lik', 'Obs', 'Rem') - Out=Out1; - save([OutputDirectoryName,'\',fname_,'_prior_glue_smooth'], 'Out', 'Sam', 'Lik', 'Obs', 'Rem') + Info.TypeofSample='prior'; + save([OutputDirectoryName,'\',fname_,'_glue_prior'], 'Out', 'Sam', 'Lik', 'Obs', 'Rem','Info', 'Exo') +% save([OutputDirectoryName,'\',fname_,'_prior_glue'], 'Out', 'Sam', 'Lik', 'Obs', 'Rem') +% Out=Out1; +% save([OutputDirectoryName,'\',fname_,'_prior_glue_smooth'], 'Out', 'Sam', 'Lik', 'Obs', 'Rem') else - save([OutputDirectoryName,'\',fname_,'_mc_glue'], 'Out', 'Sam', 'Lik', 'Obs', 'Rem') - Out=Out1; - save([OutputDirectoryName,'\',fname_,'_mc_glue_smooth'], 'Out', 'Sam', 'Lik', 'Obs', 'Rem') + Info.TypeofSample='mc'; + save([OutputDirectoryName,'\',fname_,'_glue_mc'], 'Out', 'Sam', 'Lik', 'Obs', 'Rem','Info', 'Exo') +% save([OutputDirectoryName,'\',fname_,'_mc_glue'], 'Out', 'Sam', 'Lik', 'Obs', 'Rem') +% Out=Out1; +% save([OutputDirectoryName,'\',fname_,'_mc_glue_smooth'], 'Out', 'Sam', 'Lik', 'Obs', 'Rem') end end diff --git a/matlab/filt_mc_.m b/matlab/filt_mc_.m index 8bd219ff4..f3fc45469 100644 --- a/matlab/filt_mc_.m +++ b/matlab/filt_mc_.m @@ -1,29 +1,42 @@ -function [rmse_MC, ixx] = filt_mc_(vvarvecm, loadSA, pfilt, alpha, alpha2, OutDir, istart, alphaPC) +function [rmse_MC, ixx] = filt_mc_(OutDir) +% function [rmse_MC, ixx] = filt_mc_(OutDir) % copyright Marco Ratto 2006 +% inputs (from opt_gsa structure) +% vvarvecm = options_gsa_.var_rmse; +% loadSA = options_gsa_.load_rmse; +% pfilt = options_gsa_.pfilt_rmse; +% alpha = options_gsa_.alpha_rmse; +% alpha2 = options_gsa_.alpha2_rmse; +% istart = options_gsa_.istart_rmse; +% alphaPC = 0.5; +% +% Part of the Sensitivity Analysis Toolbox for DYNARE +% +% Written by Marco Ratto, 2006 +% Joint Research Centre, The European Commission, +% (http://eemc.jrc.ec.europa.eu/), +% marco.ratto@jrc.it +% +% Disclaimer: This software is not subject to copyright protection and is in the public domain. +% It is an experimental system. The Joint Research Centre of European Commission +% assumes no responsibility whatsoever for its use by other parties +% and makes no guarantees, expressed or implied, about its quality, reliability, or any other +% characteristic. We would appreciate acknowledgement if the software is used. +% Reference: +% M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006. +% + global bayestopt_ estim_params_ M_ options_ oo_ -if nargin<1 | isempty(vvarvecm), - vvarvecm = options_.varobs; -end -if nargin<2, - loadSA=0; -end -if nargin<3 | isempty(pfilt), - pfilt=0.1; % cut the best 10% of runs -end -if nargin<4 | isempty(alpha), - alpha=0.002; -end -if nargin<5 | isempty(alpha2), - alpha2=0.5; -end -if nargin<7 | isempty(istart), - istart=1; -end -if nargin<8, - alphaPC=0.5; -end - +options_gsa_=options_.opt_gsa; +vvarvecm = options_gsa_.var_rmse; +loadSA = options_gsa_.load_rmse; +pfilt = options_gsa_.pfilt_rmse; +alpha = options_gsa_.alpha_rmse; +alpha2 = options_gsa_.alpha2_rmse; +istart = options_gsa_.istart_rmse; +alphaPC = 0.5; + fname_ = M_.fname; lgy_ = M_.endo_names; dr_ = oo_.dr; @@ -35,13 +48,18 @@ disp('for the fit of EACH observed series ...') disp(' ') disp('Deleting old SA figures...') a=dir([OutDir,'\*.*']); +tmp1='0'; if options_.opt_gsa.ppost, - tmp=['_SA_fit_post']; + tmp=['_rmse_post']; else if options_.opt_gsa.pprior - tmp=['_SA_fit_prior']; + tmp=['_rmse_prior']; else - tmp=['_SA_fit_mc']; + tmp=['_rmse_mc']; + end + if options_gsa_.lik_only, + tmp1 = [tmp,'_post_SA']; + tmp = [tmp,'_lik_SA']; end end for j=1:length(a), @@ -49,28 +67,25 @@ for j=1:length(a), disp(a(j).name) delete([OutDir,'\',a(j).name]) end, + if strmatch([fname_,tmp1],a(j).name), + disp(a(j).name) + delete([OutDir,'\',a(j).name]) + end, end disp('done !') nshock=estim_params_.nvx + estim_params_.nvn + estim_params_.ncx + estim_params_.ncn; npar=estim_params_.np; -for j=1:npar+nshock, - if j>nshock - if isfield(oo_,'posterior_mode'), - xparam1(j)=oo_.posterior_mode.parameters.(bayestopt_.name{j}); - end - if isfield(oo_,'posterior_mean'), - xparam1_mean(j)=oo_.posterior_mean.parameters.(bayestopt_.name{j}); - end - else - if isfield(oo_,'posterior_mode'), - xparam1(j)=oo_.posterior_mode.shocks_std.(bayestopt_.name{j}); - end - if isfield(oo_,'posterior_mean'), - xparam1_mean(j)=oo_.posterior_mean.shocks_std.(bayestopt_.name{j}); - end - end +load(options_.mode_file,'xparam1') +if options_.opt_gsa.ppost, + c=load([fname_,'_mean'],'xparam1'); + xparam1_mean=c.xparam1; + clear c +elseif ~isempty(ls([fname_,'_mean.mat'])) + c=load([fname_,'_mean'],'xparam1'); + xparam1_mean=c.xparam1; + clear c end if options_.opt_gsa.ppost, @@ -94,11 +109,14 @@ if ~loadSA, steady_; ys_mean=oo_.steady_state; end - eval(options_.datafile) +% eval(options_.datafile) + obs = dat_fil_(options_.datafile); if ~options_.opt_gsa.ppost - load([OutDir,'\',fnamtmp]); + load([OutDir,'\',fnamtmp],'x','logpo2','stock_gend','stock_data'); + logpo2=-logpo2; else - load([DirectoryName '/' M_.fname '_data.mat']); + %load([DirectoryName '/' M_.fname '_data.mat']); + [stock_gend, stock_data] = read_data; filfilt = dir([DirectoryName '/' M_.fname '_filter*.mat']); filparam = dir([DirectoryName '/' M_.fname '_param*.mat']); x=[]; @@ -114,7 +132,6 @@ if ~loadSA, clear stock stock_logpo stock_ys; end end - logpo2=-logpo2; end nruns=size(x,1); nfilt=floor(pfilt*nruns); @@ -125,15 +142,17 @@ if ~loadSA, nobs=options_.nobs; for i=1:size(vvarvecm,1), vj=deblank(vvarvecm(i,:)); + eval(['vobs =obs.',vj,'(fobs:fobs-1+nobs);']) if options_.prefilter == 1 %eval([vj,'=',vj,'-bayestopt_.mean_varobs(i);']) - eval([vj,'=',vj,'-mean(',vj,',1);']) + %eval([vj,'=',vj,'-mean(',vj,',1);']) + vobs = vobs-mean(vobs,1); end jxj = strmatch(vj,lgy_(dr_.order_var,:),'exact'); js = strmatch(vj,lgy_,'exact'); if exist('xparam1','var') - eval(['rmse_mode(i) = sqrt(mean((',vj,'(fobs-1+istart:fobs-1+nobs)-oo_.steady_state(js)-oo_.FilteredVariables.',vj,'(istart:end-1)).^2));']) + eval(['rmse_mode(i) = sqrt(mean((vobs(istart:end)-oo_.steady_state(js)-oo_.FilteredVariables.',vj,'(istart:end-1)).^2));']) end y0=zeros(nobs+1,nruns); if options_.opt_gsa.ppost @@ -151,19 +170,31 @@ if ~loadSA, clear stock; end else - y0 = squeeze(stock_filter(:,jxj,:)) + ... - kron(stock_ys(js,:),ones(size(stock_filter,1),1)); + filfilt=ls([OutDir,'\',fnamtmp,'_*.mat']); + nbb=0; + for j=1:size(filfilt,1), + load([OutDir,'\',fnamtmp,'_',num2str(j),'.mat'],'stock_filter','stock_ys'); + nb = size(stock_filter,3); + y0(:,nbb+1:nbb+nb) = squeeze(stock_filter(jxj,:,:)) + ... + kron(stock_ys(:,js)',ones(nobs+1,1)); + %y0(:,:,size(y0,3):size(y0,3)+size(stock,3))=stock; + nbb=nbb+nb; + clear stock_filter; + end + +% y0 = squeeze(stock_filter(:,jxj,:)) + ... +% kron(stock_ys(js,:),ones(size(stock_filter,1),1)); end y0M=mean(y0,2); for j=1:nruns, - eval(['rmse_MC(j,i) = sqrt(mean((',vj,'(fobs-1+istart:fobs-1+nobs)-y0(istart:end-1,j)).^2));']) + rmse_MC(j,i) = sqrt(mean((vobs(istart:end)-y0(istart:end-1,j)).^2)); end if exist('xparam1_mean','var') %eval(['rmse_pmean(i) = sqrt(mean((',vj,'(fobs-1+istart:fobs-1+nobs)-y0M(istart:end-1)).^2));']) [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK] = DsgeSmoother(xparam1_mean,stock_gend,stock_data); y0 = ahat(jxj,:)' + ... kron(ys_mean(js,:),ones(size(ahat,2),1)); - eval(['rmse_pmean(i) = sqrt(mean((',vj,'(fobs-1+istart:fobs-1+nobs)-y0(istart:end-1)).^2));']) + rmse_pmean(i) = sqrt(mean((vobs(istart:end)-y0(istart:end-1)).^2)); end end clear stock_filter; @@ -171,7 +202,7 @@ if ~loadSA, for j=1:nruns, lnprior(j,1) = priordens(x(j,:),bayestopt_.pshape,bayestopt_.p1,bayestopt_.p2,bayestopt_.p3,bayestopt_.p4); end - likelihood=logpo2(:)+lnprior(:); + likelihood=logpo2(:)-lnprior(:); disp('... done!') if options_.opt_gsa.ppost @@ -180,7 +211,11 @@ if ~loadSA, if options_.opt_gsa.lik_only save([OutDir,'\',fnamtmp], 'likelihood', '-append') else - save([OutDir,'\',fnamtmp], 'likelihood', 'rmse_MC', 'rmse_mode','rmse_pmean','-append') + if exist('xparam1_mean','var') + save([OutDir,'\',fnamtmp], 'likelihood', 'rmse_MC', 'rmse_mode','rmse_pmean','-append') + else + save([OutDir,'\',fnamtmp], 'likelihood', 'rmse_MC', 'rmse_mode','-append') + end end end else @@ -189,7 +224,7 @@ else else load([OutDir,'\',fnamtmp],'x','logpo2','likelihood','rmse_MC','rmse_mode','rmse_pmean'); end - lnprior=likelihood(:)-logpo2(:); + lnprior=logpo2(:)-likelihood(:); nruns=size(x,1); nfilt=floor(pfilt*nruns); end @@ -197,21 +232,21 @@ end nfilt0=nfilt*ones(size(vvarvecm,1),1); logpo2=logpo2(:); if ~options_.opt_gsa.ppost - [dum, ipost]=sort(logpo2); - [dum, ilik]=sort(likelihood); + [dum, ipost]=sort(-logpo2); + [dum, ilik]=sort(-likelihood); end if ~options_.opt_gsa.ppost & options_.opt_gsa.lik_only if options_.opt_gsa.pprior - anam='SA_fit_prior_post'; + anam='rmse_prior_post'; else - anam='SA_fit_mc_post'; + anam='rmse_mc_post'; end stab_map_1(x, ipost(1:nfilt), ipost(nfilt+1:end), anam, 1,[],OutDir); stab_map_2(x(ipost(1:nfilt),:),alpha2,anam, OutDir); if options_.opt_gsa.pprior - anam='SA_fit_prior_lik'; + anam='rmse_prior_lik'; else - anam='SA_fit_mc_lik'; + anam='rmse_mc_lik'; end stab_map_1(x, ilik(1:nfilt), ilik(nfilt+1:end), anam, 1,[],OutDir); stab_map_2(x(ilik(1:nfilt),:),alpha2,anam, OutDir); @@ -258,18 +293,18 @@ for i=1:size(vvarvecm,1), title(vvarvecm(i,:)) if mod(i,9)==0 | i==size(vvarvecm,1) if options_.opt_gsa.ppost - saveas(gcf,[OutDir,'\',fname_,'_SA_fit_post_lnprior',int2str(ifig)]) - eval(['print -depsc2 ' OutDir '\' fname_ '_SA_fit_post_lnprior',int2str(ifig)]); - eval(['print -dpdf ' OutDir '\' fname_ '_SA_fit_post_lnprior',int2str(ifig)]); + saveas(gcf,[OutDir,'\',fname_,'_rmse_post_lnprior',int2str(ifig)]) + eval(['print -depsc2 ' OutDir '\' fname_ '_rmse_post_lnprior',int2str(ifig)]); + eval(['print -dpdf ' OutDir '\' fname_ '_rmse_post_lnprior',int2str(ifig)]); else if options_.opt_gsa.pprior - saveas(gcf,[OutDir,'\',fname_,'_SA_fit_prior_lnprior',int2str(ifig)]) - eval(['print -depsc2 ' OutDir '\' fname_ '_SA_fit_prior_lnprior',int2str(ifig)]); - eval(['print -dpdf ' OutDir '\' fname_ '_SA_fit_prior_lnprior',int2str(ifig)]); + saveas(gcf,[OutDir,'\',fname_,'_rmse_prior_lnprior',int2str(ifig)]) + eval(['print -depsc2 ' OutDir '\' fname_ '_rmse_prior_lnprior',int2str(ifig)]); + eval(['print -dpdf ' OutDir '\' fname_ '_rmse_prior_lnprior',int2str(ifig)]); else - saveas(gcf,[OutDir,'\',fname_,'_SA_fit_mc_lnprior',int2str(ifig)]) - eval(['print -depsc2 ' OutDir '\' fname_ '_SA_fit_mc_lnprior',int2str(ifig)]); - eval(['print -dpdf ' OutDir '\' fname_ '_SA_fit_mc_lnprior',int2str(ifig)]); + saveas(gcf,[OutDir,'\',fname_,'_rmse_mc_lnprior',int2str(ifig)]) + eval(['print -depsc2 ' OutDir '\' fname_ '_rmse_mc_lnprior',int2str(ifig)]); + eval(['print -dpdf ' OutDir '\' fname_ '_rmse_mc_lnprior',int2str(ifig)]); end end close(gcf) @@ -288,20 +323,23 @@ for i=1:size(vvarvecm,1), h=cumplot(likelihood(ixx(nfilt0(i)+1:end,i))); set(h,'color','green') title(vvarvecm(i,:)) + if options_.opt_gsa.ppost==0, + set(gca,'xlim',[min( likelihood(ixx(1:nfilt0(i),i)) ) max( likelihood(ixx(1:nfilt0(i),i)) )]) + end if mod(i,9)==0 | i==size(vvarvecm,1) if options_.opt_gsa.ppost - saveas(gcf,[OutDir,'\',fname_,'_SA_fit_post_lnlik',int2str(ifig)]) - eval(['print -depsc2 ' OutDir '\' fname_ '_SA_fit_post_lnlik',int2str(ifig)]); - eval(['print -dpdf ' OutDir '\' fname_ '_SA_fit_post_lnlik',int2str(ifig)]); + saveas(gcf,[OutDir,'\',fname_,'_rmse_post_lnlik',int2str(ifig)]) + eval(['print -depsc2 ' OutDir '\' fname_ '_rmse_post_lnlik',int2str(ifig)]); + eval(['print -dpdf ' OutDir '\' fname_ '_rmse_post_lnlik',int2str(ifig)]); else if options_.opt_gsa.pprior - saveas(gcf,[OutDir,'\',fname_,'_SA_fit_prior_lnlik',int2str(ifig)]) - eval(['print -depsc2 ' OutDir '\' fname_ '_SA_fit_prior_lnlik',int2str(ifig)]); - eval(['print -dpdf ' OutDir '\' fname_ '_SA_fit_prior_lnlik',int2str(ifig)]); + saveas(gcf,[OutDir,'\',fname_,'_rmse_prior_lnlik',int2str(ifig)]) + eval(['print -depsc2 ' OutDir '\' fname_ '_rmse_prior_lnlik',int2str(ifig)]); + eval(['print -dpdf ' OutDir '\' fname_ '_rmse_prior_lnlik',int2str(ifig)]); else - saveas(gcf,[OutDir,'\',fname_,'_SA_fit_mc_lnlik',int2str(ifig)]) - eval(['print -depsc2 ' OutDir '\' fname_ '_SA_fit_mc_lnlik',int2str(ifig)]); - eval(['print -dpdf ' OutDir '\' fname_ '_SA_fit_mc_lnlik',int2str(ifig)]); + saveas(gcf,[OutDir,'\',fname_,'_rmse_mc_lnlik',int2str(ifig)]) + eval(['print -depsc2 ' OutDir '\' fname_ '_rmse_mc_lnlik',int2str(ifig)]); + eval(['print -dpdf ' OutDir '\' fname_ '_rmse_mc_lnlik',int2str(ifig)]); end end close(gcf) @@ -320,20 +358,23 @@ for i=1:size(vvarvecm,1), h=cumplot(logpo2(ixx(nfilt0(i)+1:end,i))); set(h,'color','green') title(vvarvecm(i,:)) + if options_.opt_gsa.ppost==0, + set(gca,'xlim',[min( logpo2(ixx(1:nfilt0(i),i)) ) max( logpo2(ixx(1:nfilt0(i),i)) )]) + end if mod(i,9)==0 | i==size(vvarvecm,1) if options_.opt_gsa.ppost - saveas(gcf,[OutDir,'\',fname_,'_SA_fit_post_lnpost',int2str(ifig)]) - eval(['print -depsc2 ' OutDir '\' fname_ '_SA_fit_post_lnpost',int2str(ifig)]); - eval(['print -dpdf ' OutDir '\' fname_ '_SA_fit_post_lnpost',int2str(ifig)]); + saveas(gcf,[OutDir,'\',fname_,'_rmse_post_lnpost',int2str(ifig)]) + eval(['print -depsc2 ' OutDir '\' fname_ '_rmse_post_lnpost',int2str(ifig)]); + eval(['print -dpdf ' OutDir '\' fname_ '_rmse_post_lnpost',int2str(ifig)]); else if options_.opt_gsa.pprior - saveas(gcf,[OutDir,'\',fname_,'_SA_fit_prior_lnpost',int2str(ifig)]) - eval(['print -depsc2 ' OutDir '\' fname_ '_SA_fit_prior_lnpost',int2str(ifig)]); - eval(['print -dpdf ' OutDir '\' fname_ '_SA_fit_prior_lnpost',int2str(ifig)]); + saveas(gcf,[OutDir,'\',fname_,'_rmse_prior_lnpost',int2str(ifig)]) + eval(['print -depsc2 ' OutDir '\' fname_ '_rmse_prior_lnpost',int2str(ifig)]); + eval(['print -dpdf ' OutDir '\' fname_ '_rmse_prior_lnpost',int2str(ifig)]); else - saveas(gcf,[OutDir,'\',fname_,'_SA_fit_mc_lnpost',int2str(ifig)]) - eval(['print -depsc2 ' OutDir '\' fname_ '_SA_fit_mc_lnpost',int2str(ifig)]); - eval(['print -dpdf ' OutDir '\' fname_ '_SA_fit_mc_lnpost',int2str(ifig)]); + saveas(gcf,[OutDir,'\',fname_,'_rmse_mc_lnpost',int2str(ifig)]) + eval(['print -depsc2 ' OutDir '\' fname_ '_rmse_mc_lnpost',int2str(ifig)]); + eval(['print -dpdf ' OutDir '\' fname_ '_rmse_mc_lnpost',int2str(ifig)]); end end close(gcf) @@ -366,11 +407,11 @@ vvarvecm=vvarvecm(ivar,:); rmse_MC=rmse_MC(:,ivar); disp(' ') -if options_.opt_gsa.ppost==0 & options_.opt_gsa.pprior, +% if options_.opt_gsa.ppost==0 & options_.opt_gsa.pprior, disp(['Sample filtered the ',num2str(pfilt*100),'% best RMSE''s for each observed series ...' ]) -else - disp(['Sample filtered the best RMSE''s smaller than RMSE at the posterior mean ...' ]) -end +% else +% disp(['Sample filtered the best RMSE''s smaller than RMSE at the posterior mean ...' ]) +% end % figure, boxplot(rmse_MC) % set(gca,'xticklabel',vvarvecm) % saveas(gcf,[fname_,'_SA_RMSE']) @@ -462,18 +503,18 @@ for ix=1:ceil(length(nsnam)/6), set(findobj(get(h0,'children'),'type','text'),'interpreter','none') title([pnam{nsnam(j)}],'interpreter','none') if options_.opt_gsa.ppost - saveas(gcf,[OutDir,'\',fname_,'_SA_fit_post_',num2str(ix)]) - eval(['print -depsc2 ' OutDir '\' fname_ '_SA_fit_post_' int2str(ix)]); - eval(['print -dpdf ' OutDir '\' fname_ '_SA_fit_post_' int2str(ix)]); + saveas(gcf,[OutDir,'\',fname_,'_rmse_post_',num2str(ix)]) + eval(['print -depsc2 ' OutDir '\' fname_ '_rmse_post_' int2str(ix)]); + eval(['print -dpdf ' OutDir '\' fname_ '_rmse_post_' int2str(ix)]); else if options_.opt_gsa.pprior - saveas(gcf,[OutDir,'\',fname_,'_SA_fit_prior_',num2str(ix)]) - eval(['print -depsc2 ' OutDir '\' fname_ '_SA_fit_prior_' int2str(ix)]); - eval(['print -dpdf ' OutDir '\' fname_ '_SA_fit_prior_' int2str(ix)]); + saveas(gcf,[OutDir,'\',fname_,'_rmse_prior_',num2str(ix)]) + eval(['print -depsc2 ' OutDir '\' fname_ '_rmse_prior_' int2str(ix)]); + eval(['print -dpdf ' OutDir '\' fname_ '_rmse_prior_' int2str(ix)]); else - saveas(gcf,[OutDir,'\',fname_,'_SA_fit_mc_',num2str(ix)]) - eval(['print -depsc2 ' OutDir '\' fname_ '_SA_fit_mc_' int2str(ix)]); - eval(['print -dpdf ' OutDir '\' fname_ '_SA_fit_mc_' int2str(ix)]); + saveas(gcf,[OutDir,'\',fname_,'_rmse_mc_',num2str(ix)]) + eval(['print -depsc2 ' OutDir '\' fname_ '_rmse_mc_' int2str(ix)]); + eval(['print -dpdf ' OutDir '\' fname_ '_rmse_mc_' int2str(ix)]); end end end @@ -521,18 +562,18 @@ for ix=1:ceil(length(nsnam)/6), set(findobj(get(h0,'children'),'type','text'),'interpreter','none') title([pnam{nsnam(j)}],'interpreter','none') if options_.opt_gsa.ppost - saveas(gcf,[OutDir,'\',fname_,'_SA_fit_post_dens_',num2str(ix)]) - eval(['print -depsc2 ' OutDir '\' fname_ '_SA_fit_post_dens_' int2str(ix)]); - eval(['print -dpdf ' OutDir '\' fname_ '_SA_fit_post_dens_' int2str(ix)]); + saveas(gcf,[OutDir,'\',fname_,'_rmse_post_dens_',num2str(ix)]) + eval(['print -depsc2 ' OutDir '\' fname_ '_rmse_post_dens_' int2str(ix)]); + eval(['print -dpdf ' OutDir '\' fname_ '_rmse_post_dens_' int2str(ix)]); else if options_.opt_gsa.pprior - saveas(gcf,[OutDir,'\',fname_,'_SA_fit_prior_dens_',num2str(ix)]) - eval(['print -depsc2 ' OutDir '\' fname_ '_SA_fit_prior_dens_' int2str(ix)]); - eval(['print -dpdf ' OutDir '\' fname_ '_SA_fit_prior_dens_' int2str(ix)]); + saveas(gcf,[OutDir,'\',fname_,'_rmse_prior_dens_',num2str(ix)]) + eval(['print -depsc2 ' OutDir '\' fname_ '_rmse_prior_dens_' int2str(ix)]); + eval(['print -dpdf ' OutDir '\' fname_ '_rmse_prior_dens_' int2str(ix)]); else - saveas(gcf,[OutDir,'\',fname_,'_SA_fit_mc_dens_',num2str(ix)]) - eval(['print -depsc2 ' OutDir '\' fname_ '_SA_fit_mc_dens_' int2str(ix)]); - eval(['print -dpdf ' OutDir '\' fname_ '_SA_fit_mc_dens_' int2str(ix)]); + saveas(gcf,[OutDir,'\',fname_,'_rmse_mc_dens_',num2str(ix)]) + eval(['print -depsc2 ' OutDir '\' fname_ '_rmse_mc_dens_' int2str(ix)]); + eval(['print -dpdf ' OutDir '\' fname_ '_rmse_mc_dens_' int2str(ix)]); end end end @@ -576,7 +617,7 @@ close all % title([pnam{np(i)},'. K-S prob ', num2str(PP(np(i),j))],'interpreter','none') % xlabel('') % if mod(i,12)==0 | i==nsx(j), -% saveas(gcf,[fname_,'_SA_fit_',deblank(vvarvecm(j,:)),'_',int2str(nfig)]) +% saveas(gcf,[fname_,'_rmse_',deblank(vvarvecm(j,:)),'_',int2str(nfig)]) % close(gcf) % end % end @@ -605,12 +646,12 @@ disp('Starting bivariate analysis:') for i=1:size(vvarvecm,1) if options_.opt_gsa.ppost - fnam = ['SA_fit_post_',deblank(vvarvecm(i,:))]; + fnam = ['rmse_post_',deblank(vvarvecm(i,:))]; else if options_.opt_gsa.pprior - fnam = ['SA_fit_prior_',deblank(vvarvecm(i,:))]; + fnam = ['rmse_prior_',deblank(vvarvecm(i,:))]; else - fnam = ['SA_fit_mc_',deblank(vvarvecm(i,:))]; + fnam = ['rmse_mc_',deblank(vvarvecm(i,:))]; end end stab_map_2(x(ixx(1:nfilt0(i),i),:),alpha2,fnam, OutDir); diff --git a/matlab/ghx2transition.m b/matlab/ghx2transition.m new file mode 100644 index 000000000..2f10a5f3b --- /dev/null +++ b/matlab/ghx2transition.m @@ -0,0 +1,41 @@ +function [A,B] = ghx2transition(mm,iv,ic,aux) +% [A,B] = ghx2transition(mm,iv,ic,aux) +% +% Adapted by M. Ratto from kalman_transition_matrix.m +% (kalman_transition_matrix.m is part of DYNARE, copyright M. Juillard) +% +% Part of the Sensitivity Analysis Toolbox for DYNARE +% +% Written by Marco Ratto, 2006 +% Joint Research Centre, The European Commission, +% (http://eemc.jrc.ec.europa.eu/), +% marco.ratto@jrc.it +% +% Disclaimer: This software is not subject to copyright protection and is in the public domain. +% It is an experimental system. The Joint Research Centre of European Commission +% assumes no responsibility whatsoever for its use by other parties +% and makes no guarantees, expressed or implied, about its quality, reliability, or any other +% characteristic. We would appreciate acknowledgement if the software is used. +% Reference: +% M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006. +% + +global M_ + + [nr1, nc1] = size(mm); + ghx = mm(:, [1:(nc1-M_.exo_nbr)]); + ghu = mm(:, [(nc1-M_.exo_nbr+1):end] ); + n_iv = length(iv); + n_ir1 = size(aux,1); + nr = n_iv + n_ir1; + + A = zeros(nr,nr); + B = zeros(nr,M_.exo_nbr); + + i_n_iv = 1:n_iv; + A(i_n_iv,ic) = ghx(iv,:); + if n_ir1 > 0 + A(n_iv+1:end,:) = sparse(aux(:,1),aux(:,2),ones(n_ir1,1),n_ir1,nr); + end + + B(i_n_iv,:) = ghu(iv,:); diff --git a/matlab/gsa_sdp_dyn.m b/matlab/gsa_sdp_dyn.m new file mode 100644 index 000000000..fe79361cf --- /dev/null +++ b/matlab/gsa_sdp_dyn.m @@ -0,0 +1,74 @@ +function gsa_ = gsa_sdp_dyn(varargin) + +% function gsa_ = +% gsa_sdp_dyn(y, x, tvp, nvr, ialias, nvrT, ifig, fname, pnames, yi_anal) +% 1 2 3 4 5 6 7 8 9 10 +% INPUTS +% y, the output +% x, the inputs +% tvp: tipe of estimation +% 0 = load previous estimation +% 1 = RW for all parameters +% 2 = IRW for all parameters +% -1 = RW for all parameters, fixed nvr's (no estimation) +% -2 = IRW for all parameters, fixed nvr's (no estimation) +% an array of length(x) specifies the process for each parameter +% nvr initial nvr value for estimation or fixed nvr is no estimation +% ialias to be used with LPTAU samples < 1024 +% nvrT maximum nvr allowable +% ifig: 1 plot figures, 0 no plots +% fname, file name to save the analysis +% pnames: names of the params +% yi_anal, analytical value of the first order HDMR terms (optional) +% +% OUTPUT +% gsa_.univariate.f +% gsa_.univariate.fs +% gsa_.univariate.fses +% gsa_.univariate.si +% gsa_.univariate.si_std +% gsa_.multivariate.f +% gsa_.multivariate.fs +% gsa_.multivariate.fses +% gsa_.multivariate.si +% gsa_.multivariate.si_std +% gsa_.multivariate.stat +% gsa_.x0 +% gsa_.xx +% gsa_.y +% +% f, function estimates +% fs, sorted function estimates +% fses, sorted standard error of function estimates +% si, sensitivity indices +% si_std, st. error of sensitivity indices +% stat, euristic tstat for significance of f +% xx, transformed inputs used (rank-transformed) +% x0, original inputs +% y, output +% +% Part of the Sensitivity Analysis Toolbox for DYNARE +% +% Written by Marco Ratto, 2006 +% Joint Research Centre, The European Commission, +% (http://eemc.jrc.ec.europa.eu/), +% marco.ratto@jrc.it +% +% Disclaimer: This software is not subject to copyright protection and is in the public domain. +% It is an experimental system. The Joint Research Centre of European Commission +% assumes no responsibility whatsoever for its use by other parties +% and makes no guarantees, expressed or implied, about its quality, reliability, or any other +% characteristic. We would appreciate acknowledgement if the software is used. +% Reference: +% M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006. +% + +if exist('sdr','file')==6 | exist('sdr','file')==2, + gsa_ = gsa_sdp_fn(varargin{:}); +else + disp('Download the SDP mapping routines at:') + disp('http://eemc.jrc.ec.europa.eu/softwareDYNARE-Dowload.htm') + disp(' ' ) + error('SDP mapping routines missing!') + +end diff --git a/matlab/prior_draw_gsa.m b/matlab/prior_draw_gsa.m new file mode 100644 index 000000000..8eded0618 --- /dev/null +++ b/matlab/prior_draw_gsa.m @@ -0,0 +1,125 @@ +function pdraw = prior_draw_gsa(init,rdraw,cc) +% Draws from the prior distributions +% Adapted by M. Ratto from prior_draw (of DYNARE, copyright M. Juillard), +% for use with Sensitivity Toolbox for DYNARE +% +% +% INPUTS +% o init [integer] scalar equal to 1 (first call) or 0. +% o rdraw +% o cc [double] two columns matrix (same as in +% metropolis.m), constraints over the +% parameter space (upper and lower bounds). +% +% OUTPUTS +% o pdraw [double] draw from the joint prior density. +% +% ALGORITHM +% ... +% +% SPECIAL REQUIREMENTS +% MATLAB Statistics Toolbox +% +% +% Part of the Sensitivity Analysis Toolbox for DYNARE +% +% Written by Marco Ratto, 2006 +% Joint Research Centre, The European Commission, +% (http://eemc.jrc.ec.europa.eu/), +% marco.ratto@jrc.it +% +% Disclaimer: This software is not subject to copyright protection and is in the public domain. +% It is an experimental system. The Joint Research Centre of European Commission +% assumes no responsibility whatsoever for its use by other parties +% and makes no guarantees, expressed or implied, about its quality, reliability, or any other +% characteristic. We would appreciate acknowledgement if the software is used. +% Reference: +% M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006. +% + +global M_ options_ estim_params_ bayestopt_ +persistent fname npar bounds pshape pmean pstd a b p1 p2 p3 p4 condition + +if init + nvx = estim_params_.nvx; + nvn = estim_params_.nvn; + ncx = estim_params_.ncx; + ncn = estim_params_.ncn; + np = estim_params_.np ; + npar = nvx+nvn+ncx+ncn+np; + MhDirectoryName = CheckPath('metropolis'); + fname = [ MhDirectoryName '/' M_.fname]; + pshape = bayestopt_.pshape; + pmean = bayestopt_.pmean; + pstd = bayestopt_.pstdev; + p1 = bayestopt_.p1; + p2 = bayestopt_.p2; + p3 = bayestopt_.p3; + p4 = bayestopt_.p4; + a = zeros(npar,1); + b = zeros(npar,1); + if nargin == 2 + bounds = cc; + else + bounds = kron(ones(npar,1),[-Inf Inf]); + end + for i = 1:npar + switch pshape(i) + case 3% Gaussian prior + b(i) = pstd(i)^2/(pmean(i)-p3(i)); + a(i) = (pmean(i)-p3(i))/b(i); + case 1% Beta prior + mu = (p1(i)-p3(i))/(p4(i)-p3(i)); + stdd = p2(i)/(p4(i)-p3(i)); + a(i) = (1-mu)*mu^2/stdd^2 - mu; + b(i) = a(i)*(1/mu - 1); + case 2;%Gamma prior + mu = p1(i)-p3(i); + b(i) = p2(i)^2/mu; + a(i) = mu/b(i); + case {5,4,6} + % Nothing to do here + % + % 4: Inverse gamma, type 1, prior + % p2(i) = nu + % p1(i) = s + % 6: Inverse gamma, type 2, prior + % p2(i) = nu + % p1(i) = s + % 5: Uniform prior + % p3(i) and p4(i) are used. + otherwise + disp('prior_draw :: Error!') + disp('Unknown prior shape.') + return + end + pdraw = zeros(npar,1); + end + condition = 1; + pdraw = zeros(npar,1); + return +end + + +for i = 1:npar + switch pshape(i) + case 5% Uniform prior. + pdraw(:,i) = rdraw(:,i)*(p4(i)-p3(i)) + p3(i); + case 3% Gaussian prior. + pdraw(:,i) = norminv(rdraw(:,i),pmean(i),pstd(i)); + case 2% Gamma prior. + pdraw(:,i) = gaminv(rdraw(:,i),a(i),b(i))+p3(i); + case 1% Beta distribution (TODO: generalized beta distribution) + pdraw(:,i) = betainv(rdraw(:,i),a(i),b(i))*(p4(i)-p3(i))+p3(i); + case 4% INV-GAMMA1 distribution + % TO BE CHECKED + pdraw(:,i) = sqrt(1./gaminv(rdraw(:,i),p2(i)/2,1/p1(i))); + case 6% INV-GAMMA2 distribution + % TO BE CHECKED + pdraw(:,i) = 1./gaminv(rdraw(:,i),p2(i)/2,1/p1(i)); + otherwise + % Nothing to do here. + end +end + + diff --git a/matlab/read_data.m b/matlab/read_data.m new file mode 100644 index 000000000..f54f66cb6 --- /dev/null +++ b/matlab/read_data.m @@ -0,0 +1,39 @@ +function [gend, data] = read_data +% Part of the Sensitivity Analysis Toolbox for DYNARE +% +% Written by Marco Ratto, 2006 +% Joint Research Centre, The European Commission, +% (http://eemc.jrc.ec.europa.eu/), +% marco.ratto@jrc.it +% +% Disclaimer: This software is not subject to copyright protection and is in the public domain. +% It is an experimental system. The Joint Research Centre of European Commission +% assumes no responsibility whatsoever for its use by other parties +% and makes no guarantees, expressed or implied, about its quality, reliability, or any other +% characteristic. We would appreciate acknowledgement if the software is used. +% Reference: +% M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006. +% + +global options_ bayestopt_ + +rawdata = read_variables(options_.datafile,options_.varobs,[],options_.xls_sheet,options_.xls_range); + +options_ = set_default_option(options_,'nobs',size(rawdata,1)-options_.first_obs+1); +gend = options_.nobs; + +rawdata = rawdata(options_.first_obs:options_.first_obs+gend-1,:); +if options_.loglinear == 1 & ~options_.logdata + rawdata = log(rawdata); +end +if options_.prefilter == 1 + bayestopt_.mean_varobs = mean(rawdata,1); + data = transpose(rawdata-ones(gend,1)*bayestopt_.mean_varobs); +else + data = transpose(rawdata); +end + +if ~isreal(rawdata) + error(['There are complex values in the data. Probably a wrong' ... + ' transformation']) +end diff --git a/matlab/redform_map.m b/matlab/redform_map.m new file mode 100644 index 000000000..40a9d840b --- /dev/null +++ b/matlab/redform_map.m @@ -0,0 +1,314 @@ +function redform_map(dirname) +%function redform_map(dirname) +% inputs (from opt_gsa structure +% anamendo = options_gsa_.namendo; +% anamlagendo = options_gsa_.namlagendo; +% anamexo = options_gsa_.namexo; +% iload = options_gsa_.load_redform; +% pprior = options_gsa_.pprior; +% ilog = options_gsa_.logtrans_redform; +% threshold = options_gsa_.threshold_redform; +% ksstat = options_gsa_.ksstat_redform; +% alpha2 = options_gsa_.alpha2_redform; +% +% Part of the Sensitivity Analysis Toolbox for DYNARE +% +% Written by Marco Ratto, 2006 +% Joint Research Centre, The European Commission, +% (http://eemc.jrc.ec.europa.eu/), +% marco.ratto@jrc.it +% +% Disclaimer: This software is not subject to copyright protection and is in the public domain. +% It is an experimental system. The Joint Research Centre of European Commission +% assumes no responsibility whatsoever for its use by other parties +% and makes no guarantees, expressed or implied, about its quality, reliability, or any other +% characteristic. We would appreciate acknowledgement if the software is used. +% Reference: +% M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006. +% + +global M_ oo_ estim_params_ options_ bayestopt_ + +options_gsa_ = options_.opt_gsa; + +anamendo = options_gsa_.namendo; +anamlagendo = options_gsa_.namlagendo; +anamexo = options_gsa_.namexo; +iload = options_gsa_.load_redform; +pprior = options_gsa_.pprior; +ilog = options_gsa_.logtrans_redform; +threshold = options_gsa_.threshold_redform; +ksstat = options_gsa_.ksstat_redform; +alpha2 = options_gsa_.alpha2_redform; + +pnames = M_.param_names(estim_params_.param_vals(:,1),:); +if nargin==0, + dirname=''; +end + +if pprior + load([dirname,'\',M_.fname,'_prior']); + adir=[dirname '\redform_stab']; +else + load([dirname,'\',M_.fname,'_mc']); + adir=[dirname '\redform_mc']; +end +if isempty(dir(adir)) + mkdir(adir) +end +adir0=pwd; +%cd(adir) + +nspred=size(T,2)-M_.exo_nbr; +x0=lpmat(istable,:); +[kn, np]=size(x0); +offset = length(bayestopt_.pshape)-np; +pshape = bayestopt_.pshape(offset+1:end); +pd = [bayestopt_.p1(offset+1:end) bayestopt_.p2(offset+1:end) bayestopt_.p3(offset+1:end) bayestopt_.p4(offset+1:end)]; + +nsok = length(find(M_.lead_lag_incidence(M_.maximum_lag,:))); +clear lpmat lpmat0 egg iunstable yys +js=0; +for j=1:size(anamendo,1) + namendo=deblank(anamendo(j,:)); + iendo=strmatch(namendo,M_.endo_names(oo_.dr.order_var,:),'exact'); + ifig=0; + iplo=0; + for jx=1:size(anamexo,1) + namexo=deblank(anamexo(jx,:)); + iexo=strmatch(namexo,M_.exo_names,'exact'); + + if ~isempty(iexo), + %y0=squeeze(T(iendo,iexo+nspred,istable)); + y0=squeeze(T(iendo,iexo+nspred,:)); + if (max(y0)-min(y0))>1.e-10, + if mod(iplo,9)==0, + ifig=ifig+1; + figure('name',[namendo,' vs. shocks ',int2str(ifig)]), + iplo=0; + end + iplo=iplo+1; + js=js+1; + xdir0 = [adir,'\',namendo,'_vs_', namexo]; + if ilog==0, + if isempty(threshold) + si(:,js) = redform_private(x0, y0, pshape, pd, iload, pnames, namendo, namexo, xdir0); + else + iy=find( (y0>threshold(1)) & (y0=threshold(2))); + xdir = [xdir0,'_cut']; + si(:,js) = redform_private(x0(iy,:), y0(iy), pshape, pd, iload, pnames, namendo, namexo, xdir); + delete([xdir, '\*cut*.*']) + [proba, dproba] = stab_map_1(x0, iy, iyc, 'cut',0); + indsmirnov = find(dproba>ksstat); + stab_map_1(x0, iy, iyc, 'cut',1,indsmirnov,xdir); + stab_map_2(x0(iy,:),alpha2,'cut',xdir) + stab_map_2(x0(iyc,:),alpha2,'trim',xdir) + end + else + [yy, xdir] = log_trans_(y0,xdir0); + silog(:,js) = redform_private(x0, yy, pshape, pd, iload, pnames, namendo, namexo, xdir); + end + + subplot(3,3,iplo), + if ilog, + [saso, iso] = sort(-silog(:,js)); + bar([silog(iso(1:10),js)]) + logflag='log'; + else + [saso, iso] = sort(-si(:,js)); + bar(si(iso(1:10),js)) + logflag=''; + end + %set(gca,'xticklabel',pnames(iso(1:10),:),'fontsize',8) + set(gca,'xticklabel',' ','fontsize',10) + set(gca,'xlim',[0.5 10.5]) + for ip=1:10, + text(ip,-0.02,deblank(pnames(iso(ip),:)),'rotation',90,'HorizontalAlignment','right','interpreter','none') + end + title([logflag,' ',namendo,' vs. ',namexo],'interpreter','none') + if iplo==9, + saveas(gcf,[dirname,'\',M_.fname,'_redform_', namendo,'_vs_shocks_',logflag,num2str(ifig)]) + eval(['print -depsc2 ' dirname,'\',M_.fname,'_redform_', namendo,'_vs_shocks_',logflag,num2str(ifig)]); + eval(['print -dpdf ' dirname,'\',M_.fname,'_redform_', namendo,'_vs_shocks_',logflag,num2str(ifig)]); + close(gcf) + end + + end + end + end + if iplo<9 & iplo>0 & ifig, + saveas(gcf,[dirname,'\',M_.fname,'_redform_', namendo,'_vs_shocks_',logflag,num2str(ifig)]) + eval(['print -depsc2 ' dirname,'\',M_.fname,'_redform_', namendo,'_vs_shocks_',logflag,num2str(ifig)]); + eval(['print -dpdf ' dirname,'\',M_.fname,'_redform_', namendo,'_vs_shocks_',logflag,num2str(ifig)]); + close(gcf) + end + ifig=0; + iplo=0; + for je=1:size(anamlagendo,1) + namlagendo=deblank(anamlagendo(je,:)); + ilagendo=strmatch(namlagendo,M_.endo_names(oo_.dr.order_var(oo_.dr.nstatic+1:oo_.dr.nstatic+nsok),:),'exact'); + + if ~isempty(ilagendo), + %y0=squeeze(T(iendo,ilagendo,istable)); + y0=squeeze(T(iendo,ilagendo,:)); + if (max(y0)-min(y0))>1.e-10, + if mod(iplo,9)==0, + ifig=ifig+1; + figure('name',[namendo,' vs. lags ',int2str(ifig)]), + iplo=0; + end + iplo=iplo+1; + js=js+1; + xdir0 = [adir,'\',namendo,'_vs_', namlagendo]; + if ilog==0, + if isempty(threshold) + si(:,js) = redform_private(x0, y0, pshape, pd, iload, pnames, namendo, namlagendo, xdir0); + else + iy=find( (y0>threshold(1)) & (y0=threshold(2))); + xdir = [xdir0,'_cut']; + si(:,js) = redform_private(x0(iy,:), y0(iy), pshape, pd, iload, pnames, namendo, namlagendo, xdir); + delete([xdir, '\*cut*.*']) + [proba, dproba] = stab_map_1(x0, iy, iyc, 'cut',0); + indsmirnov = find(dproba>ksstat); + stab_map_1(x0, iy, iyc, 'cut',1,indsmirnov,xdir); + stab_map_2(x0(iy,:),alpha2,'cut',xdir) + stab_map_2(x0(iyc,:),alpha2,'trim',xdir) + end + else + [yy, xdir] = log_trans_(y0,xdir0); + silog(:,js) = redform_private(x0, yy, pshape, pd, iload, pnames, namendo, namlagendo, xdir); + end + + subplot(3,3,iplo), + if ilog, + [saso, iso] = sort(-silog(:,js)); + bar([silog(iso(1:10),js)]) + logflag='log'; + else + [saso, iso] = sort(-si(:,js)); + bar(si(iso(1:10),js)) + logflag=''; + end + %set(gca,'xticklabel',pnames(iso(1:10),:),'fontsize',8) + set(gca,'xticklabel',' ','fontsize',10) + set(gca,'xlim',[0.5 10.5]) + for ip=1:10, + text(ip,-0.02,deblank(pnames(iso(ip),:)),'rotation',90,'HorizontalAlignment','right','interpreter','none') + end + title([logflag,' ',namendo,' vs. ',namlagendo,'(-1)'],'interpreter','none') + if iplo==9, + saveas(gcf,[dirname,'\',M_.fname,'_redform_', namendo,'_vs_lags_',logflag,num2str(ifig)]) + eval(['print -depsc2 ' dirname,'\',M_.fname,'_redform_', namendo,'_vs_lags_',logflag,num2str(ifig)]); + eval(['print -dpdf ' dirname,'\',M_.fname,'_redform_', namendo,'_vs_lags_',logflag,num2str(ifig)]); + close(gcf) + end + + end + end + end + if iplo<9 & iplo>0 & ifig, + saveas(gcf,[dirname,'\',M_.fname,'_redform_', namendo,'_vs_lags_',logflag,num2str(ifig)]) + eval(['print -depsc2 ' dirname,'\',M_.fname,'_redform_', namendo,'_vs_lags_',logflag,num2str(ifig)]); + eval(['print -dpdf ' dirname,'\',M_.fname,'_redform_', namendo,'_vs_lags_',logflag,num2str(ifig)]); + close(gcf) + end +end + +if ilog==0, +figure, bar(si) +set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:np]) +set(gca,'xlim',[0.5 np+0.5]) +set(gca,'position',[0.13 0.2 0.775 0.7]) +for ip=1:np, + text(ip,-0.02,deblank(pnames(ip,:)),'rotation',90,'HorizontalAlignment','right','interpreter','none') +end +title('Reduced form GSA') + +saveas(gcf,[dirname,'\',M_.fname,'_redform_gsa']) +eval(['print -depsc2 ' dirname,'\',M_.fname,'_redform_gsa']); +eval(['print -dpdf ' dirname,'\',M_.fname,'_redform_gsa']); + +else +figure, bar(silog) +set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:np]) +set(gca,'xlim',[0.5 np+0.5]) +set(gca,'position',[0.13 0.2 0.775 0.7]) +for ip=1:np, + text(ip,-0.02,deblank(pnames(ip,:)),'rotation',90,'HorizontalAlignment','right','interpreter','none') +end +title('Reduced form GSA - Log-transformed elements') + +saveas(gcf,[dirname,'\',M_.fname,'_redform_gsa_log']) +eval(['print -depsc2 ' dirname,'\',M_.fname,'_redform_gsa_log']); +eval(['print -dpdf ' dirname,'\',M_.fname,'_redform_gsa_log']); + +end +function si = redform_private(x0, y0, pshape, pd, iload, pnames, namy, namx, xdir) +global bayestopt_ + +fname=[xdir,'\map']; +if iload==0, + figure, hist(y0,30), title([namy,' vs. ', namx]) + if isempty(dir(xdir)) + mkdir(xdir) + end + saveas(gcf,[xdir,'\', namy,'_vs_', namx]) + eval(['print -depsc2 ' xdir,'\', namy,'_vs_', namx]); + eval(['print -dpdf ' xdir,'\', namy,'_vs_', namx]); + close(gcf) + gsa_ = gsa_sdp_dyn(y0, x0, -2, [],[],[],1,fname, pnames); +else + gsa_ = gsa_sdp_dyn(y0, x0, 0, [],[],[],0,fname, pnames); +end +si = gsa_.multivariate.si; + +function [yy, xdir]=log_trans_(y0,xdir0) + +f=inline('skewness(log(y+lam))','lam','y'); +if ~(max(y0)<0 | min(y0)>0) + isig=1; + if skewness(y0)<0, + isig=-1; + y0=-y0; + end + n=hist(y0,10); + if n(1)>20*n(end), + try lam=fzero(f,[-min(y0)+10*eps -min(y0)+abs(median(y0))],[],y0); + catch + yl(1)=f(-min(y0)+10*eps,y0); + yl(2)=f(-min(y0)+abs(median(y0)),y0); + if abs(yl(1))0 + %yy=log(y0); + xdir=[xdir0,'_log']; + end + try lam=fzero(f,[-min(y0)+10*eps -min(y0)+median(y0)],[],y0); + catch + yl(1)=f(-min(y0)+10*eps,y0); + yl(2)=f(-min(y0)+abs(median(y0)),y0); + if abs(yl(1))0 & ifig, + saveas(gcf,[dirname,'\',M_.fname,'_', namendo,'_vs_shocks_',num2str(ifig)]) + eval(['print -depsc2 ' dirname,'\',M_.fname,'_', namendo,'_vs_shocks_',num2str(ifig)]); + eval(['print -dpdf ' dirname,'\',M_.fname,'_', namendo,'_vs_shocks_',num2str(ifig)]); + close(gcf) + end + + iplo=0; + ifig=0; + for je=1:size(anamlagendo,1) + namlagendo=deblank(anamlagendo(je,:)); + ilagendo=strmatch(namlagendo,M_.endo_names(oo_.dr.order_var(oo_.dr.nstatic+1:oo_.dr.nstatic+nsok),:),'exact'); + + if ~isempty(ilagendo), + y0=teff(T(iendo,ilagendo,:),kn,istable); + if ~isempty(y0), + if mod(iplo,9)==0, + ifig=ifig+1; + figure('name',[namendo,' vs. lagged endogenous ',int2str(ifig)]), + iplo=0; + end + iplo=iplo+1; + js=js+1; + subplot(3,3,iplo), + [SAmeas, SAMorris] = Morris_Measure_Groups(estim_params_.np, lpmat, y0,nliv); + SAM = squeeze(SAMorris(:,1)); + SA(:,js)=SAM./(max(SAM)+eps); + [saso, iso] = sort(-SA(:,js)); + bar(SA(iso(1:10),js)) + %set(gca,'xticklabel',pnames(iso(1:10),:),'fontsize',8) + set(gca,'xticklabel',' ','fontsize',10) + set(gca,'xlim',[0.5 10.5]) + for ip=1:10, + text(ip,-0.02,deblank(pnames(iso(ip),:)),'rotation',90,'HorizontalAlignment','right','interpreter','none') + end + + title([namendo,' vs. ',namlagendo,'(-1)'],'interpreter','none') + if iplo==9, + saveas(gcf,[dirname,'\',M_.fname,'_', namendo,'_vs_lags_',num2str(ifig)]) + eval(['print -depsc2 ' dirname,'\',M_.fname,'_', namendo,'_vs_lags_',num2str(ifig)]); + eval(['print -dpdf ' dirname,'\',M_.fname,'_', namendo,'_vs_lags_',num2str(ifig)]); + close(gcf) + end + end + end + end + if iplo<9 & iplo>0 & ifig, + saveas(gcf,[dirname,'\',M_.fname,'_', namendo,'_vs_lags_',num2str(ifig)]) + eval(['print -depsc2 ' dirname,'\',M_.fname,'_', namendo,'_vs_lags_',num2str(ifig)]); + eval(['print -dpdf ' dirname,'\',M_.fname,'_', namendo,'_vs_lags_',num2str(ifig)]); + close(gcf) + end +end + +figure, +%bar(SA) +boxplot(SA','whis',10,'symbol','r.') +set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:np]) +set(gca,'xlim',[0.5 np+0.5]) +set(gca,'ylim',[0 1]) +set(gca,'position',[0.13 0.2 0.775 0.7]) +for ip=1:np, + text(ip,-0.02,deblank(pnames(ip,:)),'rotation',90,'HorizontalAlignment','right','interpreter','none') +end +xlabel(' ') +ylabel('Elementary Effects') +title('Reduced form screening') + +saveas(gcf,[dirname,'\',M_.fname,'_redform_screen']) +eval(['print -depsc2 ' dirname,'\',M_.fname,'_redform_screen']); +eval(['print -dpdf ' dirname,'\',M_.fname,'_redform_screen']); + diff --git a/matlab/smirnov.m b/matlab/smirnov.m index 873504e93..a08b0869b 100644 --- a/matlab/smirnov.m +++ b/matlab/smirnov.m @@ -2,7 +2,20 @@ function [H,prob,d] = smirnov(x1 , x2 , alpha, iflag ) % Smirnov test for 2 distributions % [H,prob,d] = smirnov(x1 , x2 , alpha, iflag ) % -% Copyright (C) 2005 Marco Ratto +% Part of the Sensitivity Analysis Toolbox for DYNARE +% +% Written by Marco Ratto, 2006 +% Joint Research Centre, The European Commission, +% (http://eemc.jrc.ec.europa.eu/), +% marco.ratto@jrc.it +% +% Disclaimer: This software is not subject to copyright protection and is in the public domain. +% It is an experimental system. The Joint Research Centre of European Commission +% assumes no responsibility whatsoever for its use by other parties +% and makes no guarantees, expressed or implied, about its quality, reliability, or any other +% characteristic. We would appreciate acknowledgement if the software is used. +% Reference: +% M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006. % diff --git a/matlab/stab_map_.m b/matlab/stab_map_.m index 43307ec75..deec0e575 100644 --- a/matlab/stab_map_.m +++ b/matlab/stab_map_.m @@ -1,14 +1,11 @@ -function x0 = stab_map_(Nsam, fload, ksstat, alpha2, prepSA, pprior, ilptau, OutputDirectoryName) +function x0 = stab_map_(OutputDirectoryName) % -% function x0 = stab_map_(Nsam, fload, alpha2, prepSA, pprior) +% function x0 = stab_map_(OutputDirectoryName) % % Mapping of stability regions in the prior ranges applying % Monte Carlo filtering techniques. % -% M. Ratto, Global Sensitivity Analysis for Macroeconomic models -% I. Mapping stability, MIMEO, 2005. -% -% INPUTS +% INPUTS (from opt_gsa structure) % Nsam = MC sample size % fload = 0 to run new MC; 1 to load prevoiusly generated analysis % alpha2 = significance level for bivariate sensitivity analysis @@ -34,100 +31,120 @@ function x0 = stab_map_(Nsam, fload, ksstat, alpha2, prepSA, pprior, ilptau, Out % USES lptauSEQ, % stab_map_1, stab_map_2 % -% Copyright (C) 2005 Marco Ratto -% THIS PROGRAM WAS WRITTEN FOR MATLAB BY -% Marco Ratto, -% Unit of Econometrics and Statistics AF -% (http://www.jrc.cec.eu.int/uasa/), -% IPSC, Joint Research Centre -% The European Commission, -% TP 361, 21020 ISPRA(VA), ITALY +% Part of the Sensitivity Analysis Toolbox for DYNARE +% +% Written by Marco Ratto, 2006 +% Joint Research Centre, The European Commission, +% (http://eemc.jrc.ec.europa.eu/), % marco.ratto@jrc.it % -% ALL COPIES MUST BE PROVIDED FREE OF CHARGE AND MUST INCLUDE THIS COPYRIGHT -% NOTICE. +% Disclaimer: This software is not subject to copyright protection and is in the public domain. +% It is an experimental system. The Joint Research Centre of European Commission +% assumes no responsibility whatsoever for its use by other parties +% and makes no guarantees, expressed or implied, about its quality, reliability, or any other +% characteristic. We would appreciate acknowledgement if the software is used. +% Reference: +% M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006. % %global bayestopt_ estim_params_ dr_ options_ ys_ fname_ global bayestopt_ estim_params_ options_ oo_ M_ +opt_gsa=options_.opt_gsa; + +Nsam = opt_gsa.Nsam; +fload = opt_gsa.load_stab; +ksstat = opt_gsa.ksstat; +alpha2 = opt_gsa.alpha2_stab; +prepSA = (opt_gsa.redform | opt_gsa.identification); +pprior = opt_gsa.pprior; +ilptau = opt_gsa.ilptau; +nliv = opt_gsa.morris_nliv; +ntra = opt_gsa.morris_ntra; + dr_ = oo_.dr; -if isfield(dr_,'ghx'), +%if isfield(dr_,'ghx'), ys_ = oo_.dr.ys; - nspred = size(dr_.ghx,2); + nspred = dr_.nspred; %size(dr_.ghx,2); nboth = dr_.nboth; nfwrd = dr_.nfwrd; -end +%end fname_ = M_.fname; nshock = estim_params_.nvx; nshock = nshock + estim_params_.nvn; nshock = nshock + estim_params_.ncx; nshock = nshock + estim_params_.ncn; +lpmat0=[]; +pshape = bayestopt_.pshape(nshock+1:end); +p1 = bayestopt_.p1(nshock+1:end); +p2 = bayestopt_.p2(nshock+1:end); +p3 = bayestopt_.p3(nshock+1:end); +p4 = bayestopt_.p4(nshock+1:end); if nargin==0, - Nsam=2000; %2^13; %256; -end -if nargin<2, - fload=0; -end -if nargin<3, - ksstat=0.1; -end -if nargin<4, - alpha2=0.3; -end -if nargin<5, - prepSA=0; -end -if nargin<6, - pprior=1; -end -if nargin<7, - ilptau=1; -end -if nargin<8, OutputDirectoryName=''; end -options_.periods=0; -options_.nomoments=1; -options_.irf=0; -options_.noprint=1; - -if fload==0 | nargin<2 | isempty(fload), - if prepSA - T=zeros(size(dr_.ghx,1),size(dr_.ghx,2)+size(dr_.ghu,2),Nsam/2); - end +opt=options_; + options_.periods=0; + options_.nomoments=1; + options_.irf=0; + options_.noprint=1; + options_.simul=0; +if fload==0, +% if prepSA +% T=zeros(size(dr_.ghx,1),size(dr_.ghx,2)+size(dr_.ghu,2),Nsam/2); +% end if isfield(dr_,'ghx'), egg=zeros(length(dr_.eigval),Nsam); end yys=zeros(length(dr_.ys),Nsam); - - if estim_params_.np<52 & ilptau, - [lpmat] = lptauSEQ(Nsam,estim_params_.np); - if estim_params_.np>30 + if opt_gsa.morris + if opt_gsa.morris == 1 + [lpmat, OutFact] = Sampling_Function_2(nliv, estim_params_.np, ntra, ones(estim_params_.np, 1), zeros(estim_params_.np,1), []); + lpmat = lpmat.*(nliv-1)/nliv+1/nliv/2; + Nsam=size(lpmat,1); + elseif opt_gsa.morris==2 + lpmat = prep_ide(Nsam,estim_params_.np,5); + Nsam=size(lpmat,1); + end + else + if estim_params_.np<52 & ilptau>0, + [lpmat] = lptauSEQ(Nsam,estim_params_.np); % lptau + if estim_params_.np>30 | ilptau==2, % scrambled lptau for j=1:estim_params_.np, lpmat(:,j)=lpmat(randperm(Nsam),j); end end - else + else ilptau==0 %[lpmat] = rand(Nsam,estim_params_.np); for j=1:estim_params_.np, lpmat(:,j) = randperm(Nsam)'./(Nsam+1); %latin hypercube end + end - + end + prior_draw_gsa(1); if pprior, for j=1:nshock, lpmat0(:,j) = randperm(Nsam)'./(Nsam+1); %latin hypercube - lpmat0(:,j)=lpmat0(:,j).*(bayestopt_.ub(j)-bayestopt_.lb(j))+bayestopt_.lb(j); + if opt_gsa.prior_range + lpmat0(:,j)=lpmat0(:,j).*(bayestopt_.ub(j)-bayestopt_.lb(j))+bayestopt_.lb(j); + end end - for j=1:estim_params_.np, - lpmat(:,j)=lpmat(:,j).*(bayestopt_.ub(j+nshock)-bayestopt_.lb(j+nshock))+bayestopt_.lb(j+nshock); + if opt_gsa.prior_range + for j=1:estim_params_.np, + lpmat(:,j)=lpmat(:,j).*(bayestopt_.ub(j+nshock)-bayestopt_.lb(j+nshock))+bayestopt_.lb(j+nshock); + end + else + xx=prior_draw_gsa(0,[lpmat0 lpmat]); + lpmat0=xx(:,1:nshock); + lpmat=xx(:,nshock+1:end); + clear xx; end else % for j=1:nshock, @@ -183,7 +200,24 @@ if fload==0 | nargin<2 | isempty(fload), iwrong=zeros(1,Nsam); for j=1:Nsam, M_.params(estim_params_.param_vals(:,1)) = lpmat(j,:)'; - stoch_simul([]); + %try stoch_simul([]); + try + [Tt,Rr,SteadyState,info] = dynare_resolve(bayestopt_.restrict_var_list,... + bayestopt_.restrict_columns,... + bayestopt_.restrict_aux); + + if ~exist('T') + T=zeros(size(dr_.ghx,1),size(dr_.ghx,2)+size(dr_.ghu,2),Nsam); + end + catch + if isfield(oo_.dr,'eigval'), + oo_.dr=rmfield(oo_.dr,'eigval'); + end + if isfield(oo_.dr,'ghx'), + oo_.dr=rmfield(oo_.dr,'ghx'); + end + disp('No solution could be found'), + end dr_ = oo_.dr; if isfield(dr_,'ghx'), egg(:,j) = sort(dr_.eigval); @@ -191,9 +225,13 @@ if fload==0 | nargin<2 | isempty(fload), if prepSA jstab=jstab+1; T(:,:,jstab) = [dr_.ghx dr_.ghu]; + [A,B] = ghx2transition(squeeze(T(:,:,jstab)), ... + bayestopt_.restrict_var_list, ... + bayestopt_.restrict_columns, ... + bayestopt_.restrict_aux); end if ~exist('nspred'), - nspred = size(dr_.ghx,2); + nspred = dr_.nspred; %size(dr_.ghx,2); nboth = dr_.nboth; nfwrd = dr_.nfwrd; end @@ -212,7 +250,7 @@ if fload==0 | nargin<2 | isempty(fload), end else if exist('egg'), - egg(:,j)=ones(size(egg,1),1).*1.1; + egg(:,j)=ones(size(egg,1),1).*NaN; end iwrong(j)=j; end @@ -265,20 +303,31 @@ if fload==0 | nargin<2 | isempty(fload), % end % end % iunstable=iunstable(find(iunstable)); % unstable params + bkpprior.pshape=bayestopt_.pshape; + bkpprior.p1=bayestopt_.p1; + bkpprior.p2=bayestopt_.p2; + bkpprior.p3=bayestopt_.p3; + bkpprior.p4=bayestopt_.p4; if pprior, if ~prepSA - save([OutputDirectoryName '\' fname_ '_prior'],'lpmat','lpmat0','iunstable','istable','iindeterm','iwrong','egg','yys','nspred','nboth','nfwrd') + save([OutputDirectoryName '\' fname_ '_prior'], ... + 'bkpprior','lpmat','lpmat0','iunstable','istable','iindeterm','iwrong', ... + 'egg','yys','nspred','nboth','nfwrd') else - save([OutputDirectoryName '\' fname_ '_prior'],'lpmat','lpmat0','iunstable','istable','iindeterm','iwrong','egg','yys','T','nspred','nboth','nfwrd') + save([OutputDirectoryName '\' fname_ '_prior'], ... + 'bkpprior','lpmat','lpmat0','iunstable','istable','iindeterm','iwrong', ... + 'egg','yys','T','nspred','nboth','nfwrd') end else if ~prepSA save([OutputDirectoryName '\' fname_ '_mc'], ... - 'lpmat','lpmat0','iunstable','istable','iindeterm','iwrong','egg','yys','nspred','nboth','nfwrd') + 'lpmat','lpmat0','iunstable','istable','iindeterm','iwrong', ... + 'egg','yys','nspred','nboth','nfwrd') else save([OutputDirectoryName '\' fname_ '_mc'], ... - 'lpmat','lpmat0','iunstable','istable','iindeterm','iwrong','egg','yys','T','nspred','nboth','nfwrd') + 'lpmat','lpmat0','iunstable','istable','iindeterm','iwrong', ... + 'egg','yys','T','nspred','nboth','nfwrd') end end else @@ -289,6 +338,7 @@ else end load(filetoload,'lpmat','lpmat0','iunstable','istable','iindeterm','iwrong','egg','yys','nspred','nboth','nfwrd') Nsam = size(lpmat,1); + if prepSA & isempty(strmatch('T',who('-file', filetoload),'exact')), h = waitbar(0,'Please wait...'); @@ -297,15 +347,21 @@ else options_.irf=0; options_.noprint=1; stoch_simul([]); - T=zeros(size(dr_.ghx,1),size(dr_.ghx,2)+size(dr_.ghu,2),length(istable)); + %T=zeros(size(dr_.ghx,1),size(dr_.ghx,2)+size(dr_.ghu,2),length(istable)); ntrans=length(istable); for j=1:ntrans, M_.params(estim_params_.param_vals(:,1)) = lpmat(istable(j),:)'; - stoch_simul([]); + %stoch_simul([]); + [Tt,Rr,SteadyState,info] = dynare_resolve(bayestopt_.restrict_var_list,... + bayestopt_.restrict_columns,... + bayestopt_.restrict_aux); + if ~exist('T') + T=zeros(size(dr_.ghx,1),size(dr_.ghx,2)+size(dr_.ghu,2),ntrans); + end dr_ = oo_.dr; - T(:,:,j) = [dr_.ghx dr_.ghu]; + T(:,:,j) = [dr_.ghx dr_.ghu]; if ~exist('nspred') - nspred = size(dr_.ghx,2); + nspred = dr_.nspred; %size(dr_.ghx,2); nboth = dr_.nboth; nfwrd = dr_.nfwrd; end @@ -316,6 +372,8 @@ else end close(h) save(filetoload,'T','-append') + elseif prepSA + load(filetoload,'T') end end @@ -340,37 +398,54 @@ delete([OutputDirectoryName,'\',fname_,'_',aunstname,'_corr_*.*']); delete([OutputDirectoryName,'\',fname_,'_',aindname,'_corr_*.*']); if length(iunstable)>0 & length(iunstable)ksstat); - disp('The following parameters mostly drive acceptable behaviour') - disp(M_.param_names(estim_params_.param_vals(indstab,1),:)) + disp('Smirnov statistics in driving acceptable behaviour') + for j=1:estim_params_.np, + disp([M_.param_names(estim_params_.param_vals(j,1),:),' d-stat = ', num2str(dproba(j),3)]) + end + disp(' '); + if ~isempty(indstab) stab_map_1(lpmat, istable, iunstable, aname, 1, indstab, OutputDirectoryName); + end + ixun=iunstable(find(~ismember(iunstable,[iindeterm,iwrong]))); if ~isempty(iindeterm), - ixun=iunstable(find(~ismember(iunstable,[iindeterm,iwrong]))); [proba, dproba] = stab_map_1(lpmat, [1:Nsam], iindeterm, [aname, '_indet'],0); indindet=find(dproba>ksstat); - disp('The following parameters mostly drive indeterminacy') - disp(M_.param_names(estim_params_.param_vals(indindet,1),:)) - stab_map_1(lpmat, [1:Nsam], iindeterm, [aname, '_indet'], 1, indindet, OutputDirectoryName); - if ~isempty(ixun), - [proba, dproba] = stab_map_1(lpmat, [1:Nsam], ixun, [aname, '_unst'],0); - indunst=find(dproba>ksstat); - disp('The following parameters mostly drive instability') - disp(M_.param_names(estim_params_.param_vals(indunst,1),:)) - stab_map_1(lpmat, [1:Nsam], ixun, [aname, '_unst'], 1, indunst, OutputDirectoryName); + disp('Smirnov statistics in driving indeterminacy') + for j=1:estim_params_.np, + disp([M_.param_names(estim_params_.param_vals(j,1),:),' d-stat = ', num2str(dproba(j),3)]) + end + disp(' '); + if ~isempty(indindet) + stab_map_1(lpmat, istable, iindeterm, [aname, '_indet'], 1, indindet, OutputDirectoryName); end end + if ~isempty(ixun), + [proba, dproba] = stab_map_1(lpmat, [1:Nsam], ixun, [aname, '_unst'],0); + indunst=find(dproba>ksstat); + disp('Smirnov statistics in driving instability') + for j=1:estim_params_.np, + disp([M_.param_names(estim_params_.param_vals(j,1),:),' d-stat = ', num2str(dproba(j),3)]) + end + disp(' '); + if ~isempty(indunst) + stab_map_1(lpmat, istable, ixun, [aname, '_unst'], 1, indunst, OutputDirectoryName); + end + end + disp(' ') disp('Starting bivariate analysis:') @@ -378,14 +453,16 @@ if length(iunstable)>0 & length(iunstable)3, + stab_map_2(lpmat(iunstable,:),alpha2, auname, OutputDirectoryName); end - + if length(iindeterm)>3, + stab_map_2(lpmat(iindeterm,:),alpha2, aindname, OutputDirectoryName); + end + if length(ixun)>3, + stab_map_2(lpmat(ixun,:),alpha2, aunstname, OutputDirectoryName); + end + x0=0.5.*(bayestopt_.ub(1:nshock)-bayestopt_.lb(1:nshock))+bayestopt_.lb(1:nshock); x0 = [x0; lpmat(istable(1),:)']; if istable(end)~=Nsam @@ -403,3 +480,17 @@ else end end + + +options_.periods=opt.periods; +if isfield(opt,'nomoments'), + options_.nomoments=opt.nomoments; +end +options_.irf=opt.irf; +options_.noprint=opt.noprint; +if isfield(opt,'simul'), + options_.simul=opt.simul; +end + + + diff --git a/matlab/stab_map_1.m b/matlab/stab_map_1.m index 90557cf69..17aeaa794 100644 --- a/matlab/stab_map_1.m +++ b/matlab/stab_map_1.m @@ -14,6 +14,22 @@ function [proba, dproba] = stab_map_1(lpmat, ibehaviour, inonbehaviour, aname, i % Plots: dotted lines for BEHAVIOURAL % solid lines for NON BEHAVIOURAL % USES smirnov +% +% Part of the Sensitivity Analysis Toolbox for DYNARE +% +% Written by Marco Ratto, 2006 +% Joint Research Centre, The European Commission, +% (http://eemc.jrc.ec.europa.eu/), +% marco.ratto@jrc.it +% +% Disclaimer: This software is not subject to copyright protection and is in the public domain. +% It is an experimental system. The Joint Research Centre of European Commission +% assumes no responsibility whatsoever for its use by other parties +% and makes no guarantees, expressed or implied, about its quality, reliability, or any other +% characteristic. We would appreciate acknowledgement if the software is used. +% Reference: +% M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006. +% global estim_params_ bayestopt_ M_ options_ @@ -33,7 +49,7 @@ nshock = nshock + estim_params_.ncn; npar=size(lpmat,2); ishock= npar>estim_params_.np; -if nargin<6, +if nargin<6 | isempty(ipar), ipar=[1:npar]; end nparplot=length(ipar); @@ -68,4 +84,4 @@ for i=1:ceil(nparplot/12), eval(['print -dpdf ' dirname '\' fname_ '_' aname '_SA_' int2str(i)]); if options_.nograph, close(gcf), end end -end \ No newline at end of file +end diff --git a/matlab/stab_map_2.m b/matlab/stab_map_2.m index 711c52712..3941132a2 100644 --- a/matlab/stab_map_2.m +++ b/matlab/stab_map_2.m @@ -1,17 +1,22 @@ %function stab_map_2(x,alpha2,istab,fnam) function stab_map_2(x,alpha2,fnam, dirname) -% function stab_map_2(x,alpha2,fnam) +% function stab_map_2(x,alpha2,fnam, dirname) % -% Copyright (C) 2005 Marco Ratto -% THIS PROGRAM WAS WRITTEN FOR MATLAB BY -% Marco Ratto, -% Unit of Econometrics and Statistics AF -% (http://www.jrc.cec.eu.int/uasa/), -% IPSC, Joint Research Centre -% The European Commission, -% TP 361, 21020 ISPRA(VA), ITALY +% Part of the Sensitivity Analysis Toolbox for DYNARE +% +% Written by Marco Ratto, 2006 +% Joint Research Centre, The European Commission, +% (http://eemc.jrc.ec.europa.eu/), % marco.ratto@jrc.it % +% Disclaimer: This software is not subject to copyright protection and is in the public domain. +% It is an experimental system. The Joint Research Centre of European Commission +% assumes no responsibility whatsoever for its use by other parties +% and makes no guarantees, expressed or implied, about its quality, reliability, or any other +% characteristic. We would appreciate acknowledgement if the software is used. +% Reference: +% M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006. +% %global bayestopt_ estim_params_ dr_ options_ ys_ fname_ global bayestopt_ estim_params_ options_ oo_ M_ diff --git a/matlab/teff.m b/matlab/teff.m new file mode 100644 index 000000000..0072f21d3 --- /dev/null +++ b/matlab/teff.m @@ -0,0 +1,34 @@ +function [yt, j0]=teff(T,Nsam,istable) +% +% [yt, j0]=teff(T,Nsam,istable) +% +% Part of the Sensitivity Analysis Toolbox for DYNARE +% +% Written by Marco Ratto, 2006 +% Joint Research Centre, The European Commission, +% (http://eemc.jrc.ec.europa.eu/), +% marco.ratto@jrc.it +% +% Disclaimer: This software is not subject to copyright protection and is in the public domain. +% It is an experimental system. The Joint Research Centre of European Commission +% assumes no responsibility whatsoever for its use by other parties +% and makes no guarantees, expressed or implied, about its quality, reliability, or any other +% characteristic. We would appreciate acknowledgement if the software is used. +% Reference: +% M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006. +% + +tmax=max(T,[],3); +tmin=min(T,[],3); +[ir, ic]=(find( (tmax-tmin)>1.e-8)); +j0 = length(ir); +yt=zeros(Nsam, j0); + +for j=1:j0, + y0=squeeze(T(ir(j),ic(j),:)); + %y1=ones(size(lpmat,1),1)*NaN; + y1=ones(Nsam,1)*NaN; + y1(istable,1)=y0; + yt(:,j)=y1; +end +%clear y0 y1; diff --git a/matlab/trank.m b/matlab/trank.m new file mode 100644 index 000000000..a0c5bda99 --- /dev/null +++ b/matlab/trank.m @@ -0,0 +1,25 @@ +function yr = trank(y); +% yr = trank(y); +% yr is the rank transformation of y +% +% Part of the Sensitivity Analysis Toolbox for DYNARE +% +% Written by Marco Ratto, 2006 +% Joint Research Centre, The European Commission, +% (http://eemc.jrc.ec.europa.eu/), +% marco.ratto@jrc.it +% +% Disclaimer: This software is not subject to copyright protection and is in the public domain. +% It is an experimental system. The Joint Research Centre of European Commission +% assumes no responsibility whatsoever for its use by other parties +% and makes no guarantees, expressed or implied, about its quality, reliability, or any other +% characteristic. We would appreciate acknowledgement if the software is used. +% Reference: +% M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006. +% + +[nr, nc] = size(y); +for j=1:nc, + [dum, is]=sort(y(:,j)); + yr(is,j)=[1:nr]'./nr; +end