Merge pull request #327 from rattoma/master

Identification and sensitivity fixes

matlab/dynare_identification.m

@@ -127,6 +127,7 @@ options_.analytic_derivation=1;
 options_ = set_default_option(options_,'datafile','');
 options_.mode_compute = 0;
 options_.plot_priors = 0;
+options_.smoother=1;
 [dataset_,xparam1, M_, options_, oo_, estim_params_,bayestopt_]=dynare_estimation_init(M_.endo_names,fname_,1, M_, options_, oo_, estim_params_, bayestopt_);
 options_ident.analytic_derivation_mode = options_.analytic_derivation_mode;
 if isempty(dataset_),
@@ -156,7 +157,7 @@ end
 % results = prior_sampler(0,M_,bayestopt_,options_,oo_);
 
 if prior_exist
-    if (~isnan(bayestopt_.pshape))
+    if any(bayestopt_.pshape > 0)
         if options_ident.prior_range
             prior_draw(1,1);
         else
@@ -233,9 +234,9 @@ if iload <=0,
 %             clear xparam1
 %         end
         params = set_prior(estim_params_,M_,options_)';
-        if isnan(bayestopt_.pshape)
-        parameters = 'ML_Starting_value';
-        disp('Testing ML Starting value')
+        if all(bayestopt_.pshape == 0)
+            parameters = 'ML_Starting_value';
+            disp('Testing ML Starting value')
         else
         switch parameters
             case 'posterior_mode'
@@ -275,10 +276,57 @@ if iload <=0,
         disp('----------- ')
         disp('Parameter error:')
         disp(['The model does not solve for ', parameters, ' with error code info = ', int2str(info(1))]),
-        disp('Identification stopped ')
+        disp(' ')
+        if info(1)==1,
+        disp('info==1 %! The model doesn''t determine the current variables uniquely.')
+        elseif info(1)==2,
+        disp('info==2 %! MJDGGES returned an error code.')
+        elseif info(1)==3,
+        disp('info==3 %! Blanchard & Kahn conditions are not satisfied: no stable equilibrium. ')
+        elseif info(1)==4,
+        disp('info==4 %! Blanchard & Kahn conditions are not satisfied: indeterminacy. ')
+        elseif info(1)==5,
+        disp('info==5 %! Blanchard & Kahn conditions are not satisfied: indeterminacy due to rank failure. ')
+        elseif info(1)==6,
+        disp('info==6 %! The jacobian evaluated at the deterministic steady state is complex.')
+        elseif info(1)==19,
+        disp('info==19 %! The steadystate routine thrown an exception (inconsistent deep parameters). ')
+        elseif info(1)==20,
+        disp('info==20 %! Cannot find the steady state, info(2) contains the sum of square residuals (of the static equations). ')
+        elseif info(1)==21,
+        disp('info==21 %! The steady state is complex, info(2) contains the sum of square of imaginary parts of the steady state.')
+        elseif info(1)==22,
+        disp('info==22 %! The steady has NaNs. ')
+        elseif info(1)==23,
+        disp('info==23 %! M_.params has been updated in the steadystate routine and has complex valued scalars. ')
+        elseif info(1)==24,
+        disp('info==24 %! M_.params has been updated in the steadystate routine and has some NaNs. ')
+        elseif info(1)==30,
+        disp('info==30 %! Ergodic variance can''t be computed. ')
+        end
         disp('----------- ')
         disp(' ')
-        return,
+        if any(bayestopt_.pshape)
+            disp('Try sampling up to 50 parameter sets from the prior.')
+            kk=0;
+            while kk<50 && info(1),
+                kk=kk+1;
+                params = prior_draw();
+                [idehess_point, idemoments_point, idemodel_point, idelre_point, derivatives_info_point, info] = ...
+                    identification_analysis(params,indx,indexo,options_ident,dataset_, prior_exist, name_tex,1);
+            end
+        end
+        if info(1)
+            disp(' ')
+            disp('----------- ')
+            disp('Identification stopped:')
+            if any(bayestopt_.pshape)
+                disp('The model did not solve for any of 50 attempts of random samples from the prior')
+            end
+            disp('----------- ')
+            disp(' ')
+            return,           
+        end
     else
     idehess_point.params=params;
 %     siH = idemodel_point.siH;

matlab/dynare_sensitivity.m

@@ -128,7 +128,7 @@ 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,'alpha2_stab',0);
 options_gsa = set_default_option(options_gsa,'ksstat',0.1);
 options_gsa = set_default_option(options_gsa,'pvalue_ks',0.001);
 options_gsa = set_default_option(options_gsa,'pvalue_corr',0.001);
@@ -218,6 +218,11 @@ end
 
 if options_gsa.stab && ~options_gsa.ppost,
     x0 = stab_map_(OutputDirectoryName,options_gsa);
+    if isempty(x0),
+        disp(' ')
+        disp('Sensitivity computations stopped: no parameter set provided a unique solution')
+        return
+    end
 end
 
 % reduced form

matlab/gsa/filt_mc_.m

@@ -202,7 +202,7 @@ if ~loadSA,
         clear stock_filter;
     end
     for j=1:nruns,
-        lnprior(j,1) = priordens(x(j,:)',bayestopt_.pshape,bayestopt_.p1,bayestopt_.p2,bayestopt_.p3,bayestopt_.p4);
+        lnprior(j,1) = priordens(x(j,:)',bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7,bayestopt_.p3,bayestopt_.p4);
     end
     likelihood=logpo2(:)-lnprior(:);
     disp('... done!')

matlab/gsa/map_ident_.m

@@ -54,7 +54,11 @@ if opt_gsa.load_ident_files==0,
   mss = yys(bayestopt_.mfys,:);
   mss = teff(mss(:,istable),Nsam,istable);
   yys = teff(yys(oo_.dr.order_var,istable),Nsam,istable);
-  [vdec, cc, ac] = mc_moments(T, lpmatx, oo_.dr);
+  if exist('T'),
+      [vdec, cc, ac] = mc_moments(T, lpmatx, oo_.dr);
+  else
+      return,
+  end
 
 
   if opt_gsa.morris==2,

matlab/gsa/redform_map.m

@@ -73,6 +73,11 @@ if ~exist('T')
     else
         load([dirname,'/',M_.fname,'_mc'],'T');
     end
+    if ~exist('T'),
+        disp('The model is too large!')
+        disp('Reduced form mapping stopped!')
+        return
+    end
 end
 if isempty(dir(adir))
     mkdir(adir)

matlab/gsa/stab_map_.m

@@ -256,22 +256,36 @@ if fload==0,
         M_.params(estim_params_.param_vals(:,1)) = lpmat(j,:)';
         %try stoch_simul([]);
         try
-            [Tt,Rr,SteadyState,infox{j},M_,options_,oo_] = dynare_resolve(M_,options_,oo_,'restrict');
-            if infox{j}(1)==0 && ~exist('T'),
+            [Tt,Rr,SteadyState,info,M_,options_,oo_] = dynare_resolve(M_,options_,oo_,'restrict');
+            infox(j,1)=info(1);
+            if infox(j,1)==0 && ~exist('T'),
                 dr_=oo_.dr;
-                T=zeros(size(dr_.ghx,1),size(dr_.ghx,2)+size(dr_.ghu,2),Nsam);
+                if prepSA,
+                    try
+                        T=zeros(size(dr_.ghx,1),size(dr_.ghx,2)+size(dr_.ghu,2),Nsam);
+                    catch ME
+                        if strcmp('MATLAB:nomem',ME.identifier),
+                            prepSA=0;
+                            disp('The model is too large for storing state space matrices ...')
+                            disp('for mapping reduced form or for identification')
+                        end
+                        T=[];
+                    end
+                else
+                    T=[];
+                end
                 egg=zeros(length(dr_.eigval),Nsam);
             end
-            if infox{j}(1),
+            if infox(j,1),
 %                 disp('no solution'),
                 if isfield(oo_.dr,'ghx'),
                     oo_.dr=rmfield(oo_.dr,'ghx');
                 end
-                if (infox{j}(1)<3 || infox{j}(1)>5) && isfield(oo_.dr,'eigval'),
+                if (infox(j,1)<3 || infox(j,1)>5) && isfield(oo_.dr,'eigval'),
                     oo_.dr=rmfield(oo_.dr,'eigval');
                 end
             end
-        catch
+        catch ME
             if isfield(oo_.dr,'eigval'),
                 oo_.dr=rmfield(oo_.dr,'eigval');
             end
@@ -323,8 +337,10 @@ if fload==0,
         dyn_waitbar(j/Nsam,h,['MC iteration ',int2str(j),'/',int2str(Nsam)])
     end
     dyn_waitbar_close(h);
-    if prepSA,
+    if prepSA && jstab,
         T=T(:,:,1:jstab);
+    else
+        T=[];
     end
     istable=istable(find(istable));  % stable params
     iunstable=iunstable(find(iunstable));   % unstable params
@@ -374,22 +390,22 @@ if fload==0,
         if ~prepSA
             save([OutputDirectoryName '/' fname_ '_prior.mat'], ...
                 'bkpprior','lpmat','lpmat0','iunstable','istable','iindeterm','iwrong', ...
-                'egg','yys','nspred','nboth','nfwrd')
+                'egg','yys','nspred','nboth','nfwrd','infox')
         else
             save([OutputDirectoryName '/' fname_ '_prior.mat'], ...
                 'bkpprior','lpmat','lpmat0','iunstable','istable','iindeterm','iwrong', ...
-                'egg','yys','T','nspred','nboth','nfwrd')
+                'egg','yys','T','nspred','nboth','nfwrd','infox')
         end
 
     else
         if ~prepSA
             save([OutputDirectoryName '/' fname_ '_mc.mat'], ...
                 'lpmat','lpmat0','iunstable','istable','iindeterm','iwrong', ...
-                'egg','yys','nspred','nboth','nfwrd')
+                'egg','yys','nspred','nboth','nfwrd','infox')
         else
             save([OutputDirectoryName '/' fname_ '_mc.mat'], ...
                 'lpmat','lpmat0','iunstable','istable','iindeterm','iwrong', ...
-                'egg','yys','T','nspred','nboth','nfwrd')
+                'egg','yys','T','nspred','nboth','nfwrd','infox')
         end
     end
 else
@@ -398,7 +414,7 @@ else
     else
         filetoload=[OutputDirectoryName '/' fname_ '_mc.mat'];
     end
-    load(filetoload,'lpmat','lpmat0','iunstable','istable','iindeterm','iwrong','egg','yys','nspred','nboth','nfwrd')
+    load(filetoload,'lpmat','lpmat0','iunstable','istable','iindeterm','iwrong','egg','yys','nspred','nboth','nfwrd','infox')
     Nsam = size(lpmat,1);
     if pprior==0,
         eval(['load ' options_.mode_file '.mat;']);
@@ -467,14 +483,46 @@ delete([OutputDirectoryName,'/',fname_,'_',aunstname,'_corr_*.*']);
 delete([OutputDirectoryName,'/',fname_,'_',aindname,'_corr_*.*']);
 
 if length(iunstable)>0 && length(iunstable)<Nsam,
-    fprintf(['%4.1f%% of the prior support is stable.\n'],length(istable)/Nsam*100)
-    fprintf(['%4.1f%% of the prior support is unstable.\n'],(length(iunstable)-length(iwrong)-length(iindeterm) )/Nsam*100)
+    fprintf(['%4.1f%% of the prior support gives unique saddle-path solution.\n'],length(istable)/Nsam*100)
+    fprintf(['%4.1f%% of the prior support gives explosive dynamics.\n'],(length(iunstable)-length(iwrong)-length(iindeterm) )/Nsam*100)
     if ~isempty(iindeterm),
         fprintf(['%4.1f%% of the prior support gives indeterminacy.'],length(iindeterm)/Nsam*100)
     end
     if ~isempty(iwrong),
         disp(' ');
         disp(['For ',num2str(length(iwrong)/Nsam*100,'%1.3f'),'\% of the prior support dynare could not find a solution.'])
+        disp(' ');
+        if any(infox==1),
+            disp(['    For ',num2str(length(find(infox==1))/Nsam*100,'%1.3f'),'\% The model doesn''t determine the current variables uniquely.'])
+        end
+        if any(infox==2),
+            disp(['    For ',num2str(length(find(infox==2))/Nsam*100,'%1.3f'),'\% MJDGGES returned an error code.'])
+        end
+        if any(infox==6),
+            disp(['    For ',num2str(length(find(infox==6))/Nsam*100,'%1.3f'),'\% The jacobian evaluated at the deterministic steady state is complex.'])
+        end
+        if any(infox==19),
+            disp(['    For ',num2str(length(find(infox==19))/Nsam*100,'%1.3f'),'\% The steadystate routine thrown an exception (inconsistent deep parameters).'])
+        end
+        if any(infox==20),
+            disp(['    For ',num2str(length(find(infox==20))/Nsam*100,'%1.3f'),'\% Cannot find the steady state.'])
+        end
+        if any(infox==21),
+            disp(['    For ',num2str(length(find(infox==21))/Nsam*100,'%1.3f'),'\% The steady state is complex.'])
+        end
+        if any(infox==22),
+            disp(['    For ',num2str(length(find(infox==22))/Nsam*100,'%1.3f'),'\% The steady has NaNs.'])
+        end
+        if any(infox==23),
+            disp(['    For ',num2str(length(find(infox==23))/Nsam*100,'%1.3f'),'\% M_.params has been updated in the steadystate routine and has complex valued scalars.'])
+        end
+        if any(infox==24),
+            disp(['    For ',num2str(length(find(infox==24))/Nsam*100,'%1.3f'),'\% M_.params has been updated in the steadystate routine and has some NaNs.'])
+        end
+        if any(infox==30),
+            disp(['    For ',num2str(length(find(infox==30))/Nsam*100,'%1.3f'),'\% Ergodic variance can''t be computed.'])
+        end
+
     end
     disp(' ');
     % Blanchard Kahn
@@ -561,7 +609,7 @@ if length(iunstable)>0 && length(iunstable)<Nsam,
     end
 else
     if length(iunstable)==0,
-        disp('All parameter values in the specified ranges are stable!')
+        disp('All parameter values in the specified ranges give unique saddle-path solution!')
         x0=0.5.*(bayestopt_.ub(1:nshock)-bayestopt_.lb(1:nshock))+bayestopt_.lb(1:nshock);
         x0 = [x0; lpmat(istable(1),:)'];
     else

matlab/identification_analysis.m

@@ -48,7 +48,9 @@ persistent indH indJJ indLRE
 nparam=length(params);
 np=length(indx);
 offset=nparam-np;
-M_ = set_all_parameters(params,estim_params_,M_);
+if ~isempty(estim_params_),
+    M_ = set_all_parameters(params,estim_params_,M_);
+end
 
 nlags = options_ident.ar;
 useautocorr = options_ident.useautocorr;