dynare_sensitivity

 Frontend to the Sensitivity Analysis Toolbox for DYNARE

 Reference:
 M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006.

0001 function x0=dynare_sensitivity(options_gsa)
0002 % Frontend to the Sensitivity Analysis Toolbox for DYNARE
0003 %
0004 % Reference:
0005 % M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006.
0006 
0007 % Copyright (C) 2008-2011 Dynare Team
0008 %
0009 % This file is part of Dynare.
0010 %
0011 % Dynare is free software: you can redistribute it and/or modify
0012 % it under the terms of the GNU General Public License as published by
0013 % the Free Software Foundation, either version 3 of the License, or
0014 % (at your option) any later version.
0015 %
0016 % Dynare is distributed in the hope that it will be useful,
0017 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0018 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0019 % GNU General Public License for more details.
0020 %
0021 % You should have received a copy of the GNU General Public License
0022 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0023 
0024 global M_ options_ oo_ bayestopt_ estim_params_
0025 
0026 fname_ = M_.fname;
0027 if ~isfield(M_,'dname'),
0028     M_.dname = M_.fname;
0029 end
0030 lgy_ = M_.endo_names;
0031 x0=[];
0032 
0033 options_gsa = set_default_option(options_gsa,'datafile',[]);
0034 if isfield(options_gsa,'nograph'),
0035     options_.nograph=options_gsa.nograph;
0036 end
0037 if isfield(options_gsa,'mode_file'),
0038     options_.mode_file=options_gsa.mode_file;
0039 end
0040 
0041 options_.order = 1;
0042 
0043 if ~isempty(options_gsa.datafile) || isempty(bayestopt_),
0044     options_.datafile = options_gsa.datafile;
0045     if isfield(options_gsa,'first_obs'),
0046         options_.first_obs=options_gsa.first_obs;
0047     end
0048     if isfield(options_gsa,'nobs'),
0049         options_.nobs=options_gsa.nobs;
0050     end
0051     if isfield(options_gsa,'presample'),
0052         options_.presample=options_gsa.presample;
0053     end
0054     if isfield(options_gsa,'prefilter'),
0055         options_.prefilter=options_gsa.prefilter;
0056     end
0057     if isfield(options_gsa,'loglinear'),
0058         options_.loglinear=options_gsa.loglinear;
0059     end
0060     options_.mode_compute = 0;
0061     options_.filtered_vars = 1;
0062     options_.plot_priors = 0;
0063 %     [data,rawdata,xparam1,data_info]=dynare_estimation_init([],fname_,1);
0064     [dataset_,xparam1, M_, options_, oo_, estim_params_,bayestopt_]=dynare_estimation_init(M_.endo_names,fname_,1, M_, options_, oo_, estim_params_, bayestopt_);
0065     % computes a first linear solution to set up various variables
0066 else
0067     if isempty(options_.qz_criterium)
0068         options_.qz_criterium = 1+1e-6;
0069     end    
0070 end
0071 [make,my,day,punk,M_,options_,oo_] = dynare_resolve(M_,options_,oo_);
0072 
0073 options_gsa = set_default_option(options_gsa,'identification',0);
0074 if options_gsa.identification,
0075     options_gsa.redform=0;
0076     options_gsa = set_default_option(options_gsa,'morris',1);
0077     options_gsa = set_default_option(options_gsa,'trans_ident',0);
0078     options_gsa = set_default_option(options_gsa,'load_ident_files',0);
0079     options_gsa = set_default_option(options_gsa,'ar',1);
0080     options_gsa = set_default_option(options_gsa,'useautocorr',0);
0081     options_.ar = options_gsa.ar;
0082     if options_gsa.morris==2,
0083         if isfield(options_,'options_ident'),
0084             options_.options_ident.load_ident_files = options_gsa.load_ident_files;
0085             options_.options_ident.useautocorr = options_gsa.useautocorr;
0086             options_.options_ident.ar = options_gsa.ar;            
0087             options_ident=options_.options_ident;
0088         else
0089             options_ident=[];
0090             options_ident = set_default_option(options_ident,'load_ident_files',options_gsa.load_ident_files);
0091             options_ident = set_default_option(options_ident,'useautocorr',options_gsa.useautocorr);
0092             options_ident = set_default_option(options_ident,'ar',options_gsa.ar);
0093             options_.options_ident = options_ident;
0094         end
0095     end
0096 end
0097 
0098 % map stability
0099 options_gsa = set_default_option(options_gsa,'stab',1);
0100 options_gsa = set_default_option(options_gsa,'redform',0);
0101 options_gsa = set_default_option(options_gsa,'pprior',1);
0102 options_gsa = set_default_option(options_gsa,'prior_range',1);
0103 options_gsa = set_default_option(options_gsa,'ppost',0);
0104 options_gsa = set_default_option(options_gsa,'neighborhood_width',0);
0105 options_gsa = set_default_option(options_gsa,'ilptau',1);
0106 options_gsa = set_default_option(options_gsa,'morris',0);
0107 options_gsa = set_default_option(options_gsa,'glue',0);
0108 options_gsa = set_default_option(options_gsa,'morris_nliv',6);
0109 options_gsa = set_default_option(options_gsa,'morris_ntra',20);
0110 options_gsa = set_default_option(options_gsa,'Nsam',2048);
0111 options_gsa = set_default_option(options_gsa,'load_redform',0);
0112 options_gsa = set_default_option(options_gsa,'load_rmse',0);
0113 options_gsa = set_default_option(options_gsa,'load_stab',0);
0114 options_gsa = set_default_option(options_gsa,'alpha2_stab',0.3);
0115 options_gsa = set_default_option(options_gsa,'ksstat',0.1);
0116 options_gsa = set_default_option(options_gsa,'pvalue_ks',0.001);
0117 options_gsa = set_default_option(options_gsa,'pvalue_corr',0.001);
0118 %options_gsa = set_default_option(options_gsa,'load_mh',0);
0119 
0120 if options_gsa.redform,
0121     options_gsa.pprior=1;
0122     options_gsa.ppost=0;
0123 end
0124 
0125 if options_gsa.morris==1 || options_gsa.morris==3,
0126     if ~options_gsa.identification,
0127         options_gsa.redform=1;
0128     end
0129     options_gsa.pprior=1;
0130     options_gsa.ppost=0;
0131     %options_gsa.stab=1;
0132     options_gsa.glue=0;
0133     options_gsa.rmse=0;
0134     options_gsa.load_rmse=0;
0135     options_gsa.alpha2_stab=1;
0136     options_gsa.ksstat=1;
0137     if options_gsa.morris==3,
0138         options_gsa = set_default_option(options_gsa,'Nsam',256);
0139         OutputDirectoryName = CheckPath('gsa/identif',M_.dname);
0140     else
0141         OutputDirectoryName = CheckPath('gsa/screen',M_.dname);
0142     end
0143 else
0144     OutputDirectoryName = CheckPath('gsa',M_.dname);
0145 end
0146 
0147 options_.opt_gsa = options_gsa;
0148 
0149 if (options_gsa.load_stab || options_gsa.load_rmse || options_gsa.load_redform) ...
0150         && options_gsa.pprior,
0151     filetoload=[OutputDirectoryName '/' fname_ '_prior.mat'];
0152     if ~exist(filetoload),
0153         disp([filetoload,' not found!'])
0154         disp(['You asked to load a non existent analysis'])
0155         %options_gsa.load_stab=0;
0156         return,
0157     else
0158         if isempty(strmatch('bkpprior',who('-file', filetoload),'exact')),
0159             disp('Warning! Missing prior info for saved sample') % trap for files previous
0160             disp('The saved files are generated with previous version of GSA package') % trap for files previous
0161         else
0162             load(filetoload,'bkpprior'),
0163             if any(bayestopt_.pshape~=bkpprior.pshape) || ...
0164                     any(bayestopt_.p1~=bkpprior.p1) || ...
0165                     any(bayestopt_.p2~=bkpprior.p2) || ...
0166                     any(bayestopt_.p3(~isnan(bayestopt_.p3))~=bkpprior.p3(~isnan(bkpprior.p3))) || ...
0167                     any(bayestopt_.p4(~isnan(bayestopt_.p4))~=bkpprior.p4(~isnan(bkpprior.p4))),
0168                 disp('WARNING!')
0169                 disp('The saved sample has different priors w.r.t. to current ones!!')
0170                 disp('')
0171                 disp('Press ENTER to continue')
0172                 pause;
0173             end
0174         end
0175     end
0176 end
0177 
0178 if options_gsa.stab && ~options_gsa.ppost,
0179     x0 = stab_map_(OutputDirectoryName);
0180 end
0181 
0182 % reduced form
0183 % redform_map(namendo, namlagendo, namexo, icomp, pprior, ilog, threshold)
0184 options_gsa = set_default_option(options_gsa,'logtrans_redform',0);
0185 options_gsa = set_default_option(options_gsa,'threshold_redform',[]);
0186 options_gsa = set_default_option(options_gsa,'ksstat_redform',0.1);
0187 options_gsa = set_default_option(options_gsa,'alpha2_redform',0.3);
0188 options_gsa = set_default_option(options_gsa,'namendo',[]);
0189 options_gsa = set_default_option(options_gsa,'namlagendo',[]);
0190 options_gsa = set_default_option(options_gsa,'namexo',[]);
0191 
0192 options_.opt_gsa = options_gsa;
0193 if options_gsa.identification,
0194     map_ident_(OutputDirectoryName);
0195 end
0196 
0197 if options_gsa.redform && ~isempty(options_gsa.namendo) ...
0198         && ~options_gsa.ppost,
0199     if strmatch(':',options_gsa.namendo,'exact'),
0200         options_gsa.namendo=M_.endo_names;
0201     end
0202     if strmatch(':',options_gsa.namexo,'exact'),
0203         options_gsa.namexo=M_.exo_names;
0204     end
0205     if strmatch(':',options_gsa.namlagendo,'exact'),
0206         options_gsa.namlagendo=M_.endo_names;
0207     end
0208     options_.opt_gsa = options_gsa;
0209     if options_gsa.morris==1,
0210         redform_screen(OutputDirectoryName);
0211     else
0212         % check existence of the SS_ANOVA toolbox
0213         if ~(exist('gsa_sdp','file')==6 || exist('gsa_sdp','file')==2),
0214             disp('Download Mapping routines at:')
0215             disp('http://eemc.jrc.ec.europa.eu/softwareDYNARE-Dowload.htm')
0216             disp(' ' )
0217             error('Mapping routines missing!')
0218         end
0219 
0220         redform_map(OutputDirectoryName);
0221     end
0222 end
0223 % RMSE mapping
0224 % function [rmse_MC, ixx] = filt_mc_(vvarvecm, loadSA, pfilt, alpha, alpha2)
0225 options_gsa = set_default_option(options_gsa,'rmse',0);
0226 options_gsa = set_default_option(options_gsa,'lik_only',0);
0227 options_gsa = set_default_option(options_gsa,'var_rmse',options_.varobs);
0228 options_gsa = set_default_option(options_gsa,'pfilt_rmse',0.1);
0229 options_gsa = set_default_option(options_gsa,'istart_rmse',options_.presample+1);
0230 options_gsa = set_default_option(options_gsa,'alpha_rmse',0.002);
0231 options_gsa = set_default_option(options_gsa,'alpha2_rmse',1);
0232 options_.opt_gsa = options_gsa;
0233 if options_gsa.rmse,
0234     if ~options_gsa.ppost
0235         if options_gsa.pprior
0236             a=whos('-file',[OutputDirectoryName,'/',fname_,'_prior'],'logpo2');
0237         else
0238             a=whos('-file',[OutputDirectoryName,'/',fname_,'_mc'],'logpo2');
0239         end
0240         if exist('OCTAVE_VERSION'),
0241             aflag=0;
0242             for ja=1:length(a),
0243                 aflag=aflag+strcmp('logpo2',a(ja).name);
0244             end
0245             if aflag==0,
0246                 a=[];
0247             else
0248                 a=1;
0249             end
0250         end
0251         if isempty(a),
0252 %             dynare_MC([],OutputDirectoryName,data,rawdata,data_info);
0253             prior_posterior_statistics('gsa',dataset_);
0254             if options_.bayesian_irf
0255                 PosteriorIRF('gsa');
0256             end
0257             options_gsa.load_rmse=0;
0258             %   else
0259             %     if options_gsa.load_rmse==0,
0260             %       disp('You already saved a MC filter/smoother analysis ')
0261             %       disp('Do you want to overwrite ?')
0262             %       pause;
0263             %       if options_gsa.pprior
0264             %         delete([OutputDirectoryName,'/',fname_,'_prior_*.mat'])
0265             %       else
0266             %         delete([OutputDirectoryName,'/',fname_,'_mc_*.mat'])
0267             %       end
0268             %       dynare_MC([],OutputDirectoryName);
0269             %       options_gsa.load_rmse=0;
0270             %     end
0271             
0272         end
0273     end
0274     clear a;
0275 %     filt_mc_(OutputDirectoryName,data_info);
0276     filt_mc_(OutputDirectoryName);
0277 end
0278 
0279 
0280 if options_gsa.glue,
0281     dr_ = oo_.dr;
0282     if options_gsa.ppost
0283         load([OutputDirectoryName,'/',fname_,'_post']);
0284         DirectoryName = CheckPath('metropolis',M_.dname);
0285     else
0286         if options_gsa.pprior
0287             load([OutputDirectoryName,'/',fname_,'_prior']);
0288         else
0289             load([OutputDirectoryName,'/',fname_,'_mc']);
0290         end
0291     end
0292     if ~exist('x'),
0293         disp(['No RMSE analysis is available for current options'])
0294         disp(['No GLUE file prepared'])
0295         return,
0296     end
0297     nruns=size(x,1);
0298     gend = options_.nobs;
0299     rawdata = read_variables(options_.datafile,options_.varobs,[],options_.xls_sheet,options_.xls_range);
0300     rawdata = rawdata(options_.first_obs:options_.first_obs+gend-1,:);
0301     if options_.loglinear == 1
0302         rawdata = log(rawdata);
0303     end
0304     if options_.prefilter == 1
0305         %data = transpose(rawdata-ones(gend,1)*bayestopt_.mean_varobs);
0306         data = transpose(rawdata-ones(gend,1)*mean(rawdata,1));
0307     else
0308         data = transpose(rawdata);
0309     end
0310     
0311     Obs.data = data;
0312     Obs.time = [1:gend];
0313     Obs.num  = gend;
0314     for j=1:size(options_.varobs,1)
0315         Obs.name{j} = deblank(options_.varobs(j,:));
0316         vj=deblank(options_.varobs(j,:));
0317         
0318         jxj = strmatch(vj,lgy_(dr_.order_var,:),'exact');
0319         js = strmatch(vj,lgy_,'exact');
0320         if ~options_gsa.ppost
0321             %       y0=zeros(gend+1,nruns);
0322             %       nb = size(stock_filter,3);
0323             %       y0 = squeeze(stock_filter(:,jxj,:)) + ...
0324             %         kron(stock_ys(js,:),ones(size(stock_filter,1),1));
0325             %       Out(j).data = y0';
0326             %       Out(j).time = [1:size(y0,1)];
0327             Out(j).data = jxj;
0328             Out(j).time = [pwd,'/',OutputDirectoryName];
0329         else
0330             Out(j).data = jxj;
0331             Out(j).time = [pwd,'/',DirectoryName];
0332         end
0333         Out(j).name = vj;
0334         Out(j).ini  = 'yes';
0335         Lik(j).name = ['rmse_',vj];
0336         Lik(j).ini  = 'yes';
0337         Lik(j).isam = 1;
0338         Lik(j).data = rmse_MC(:,j)';
0339         
0340         if ~options_gsa.ppost
0341             %       y0 = squeeze(stock_smooth(:,jxj,:)) + ...
0342             %         kron(stock_ys(js,:),ones(size(stock_smooth,1),1));
0343             %       Out1(j).name = vj;
0344             %       Out1(j).ini  = 'yes';
0345             %       Out1(j).time = [1:size(y0,1)];
0346             %       Out1(j).data = y0';
0347             Out1=Out;
0348         else
0349             Out1=Out;
0350         end
0351         ismoo(j)=jxj;
0352         
0353     end
0354     jsmoo = size(options_.varobs,1);
0355     for j=1:M_.endo_nbr,
0356         if ~ismember(j,ismoo),
0357             jsmoo=jsmoo+1;
0358             vj=deblank(M_.endo_names(dr_.order_var(j),:));
0359             if ~options_gsa.ppost        
0360                 %         y0 = squeeze(stock_smooth(:,j,:)) + ...
0361                 %           kron(stock_ys(j,:),ones(size(stock_smooth,1),1));
0362                 %         Out1(jsmoo).time = [1:size(y0,1)];
0363                 %         Out1(jsmoo).data = y0';
0364                 Out1(jsmoo).data = j;
0365                 Out1(jsmoo).time = [pwd,'/',OutputDirectoryName];
0366             else
0367                 Out1(jsmoo).data = j;
0368                 Out1(jsmoo).time = [pwd,'/',DirectoryName];
0369             end
0370             Out1(jsmoo).name = vj;
0371             Out1(jsmoo).ini  = 'yes';
0372         end
0373     end
0374     tit(M_.exo_names_orig_ord,:) = M_.exo_names;
0375     for j=1:M_.exo_nbr,
0376         Exo(j).name = deblank(tit(j,:));    
0377     end
0378     if ~options_gsa.ppost
0379         Lik(size(options_.varobs,1)+1).name = 'logpo';
0380         Lik(size(options_.varobs,1)+1).ini  = 'yes';
0381         Lik(size(options_.varobs,1)+1).isam = 1;
0382         Lik(size(options_.varobs,1)+1).data = -logpo2;
0383     end
0384     Sam.name = bayestopt_.name;
0385     Sam.dim  = [size(x) 0];
0386     Sam.data = [x];
0387     
0388     Rem.id = 'Original';
0389     Rem.ind= [1:size(x,1)];
0390     
0391     Info.dynare=M_.fname;
0392     Info.order_var=dr_.order_var;
0393     Out=Out1;
0394     if options_gsa.ppost
0395         %     Info.dynare=M_.fname;
0396         %     Info.order_var=dr_.order_var;
0397         %     Out=Out1;
0398         Info.TypeofSample='post';
0399         save([OutputDirectoryName,'/',fname_,'_glue_post.mat'], 'Out', 'Sam', 'Lik', 'Obs', 'Rem','Info', 'Exo')
0400         %save([fname_,'_post_glue_smooth'], 'Out', 'Sam', 'Lik', 'Obs', 'Rem','Info')
0401         
0402     else
0403         if options_gsa.pprior
0404             Info.TypeofSample='prior';
0405             save([OutputDirectoryName,'/',fname_,'_glue_prior.mat'], 'Out', 'Sam', 'Lik', 'Obs', 'Rem','Info', 'Exo')
0406             %       save([OutputDirectoryName,'/',fname_,'_prior_glue'], 'Out', 'Sam', 'Lik', 'Obs', 'Rem')
0407             %       Out=Out1;
0408             %       save([OutputDirectoryName,'/',fname_,'_prior_glue_smooth'], 'Out', 'Sam', 'Lik', 'Obs', 'Rem')
0409         else
0410             Info.TypeofSample='mc';
0411             save([OutputDirectoryName,'/',fname_,'_glue_mc.mat'], 'Out', 'Sam', 'Lik', 'Obs', 'Rem','Info', 'Exo')
0412             %       save([OutputDirectoryName,'/',fname_,'_mc_glue'], 'Out', 'Sam', 'Lik', 'Obs', 'Rem')
0413             %       Out=Out1;
0414             %       save([OutputDirectoryName,'/',fname_,'_mc_glue_smooth'], 'Out', 'Sam', 'Lik', 'Obs', 'Rem')
0415         end
0416     end
0417     
0418 end