stab_map_

 function x0 = stab_map_(OutputDirectoryName)

 Mapping of stability regions in the prior ranges applying
 Monte Carlo filtering techniques.

 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
 [abs(corrcoef) > alpha2]
 prepSA = 1: save transition matrices for mapping reduced form
        = 0: no transition matrix saved (default)
 pprior = 1: sample from prior ranges (default): sample saved in
            _prior.mat   file
        = 0: sample from posterior ranges: sample saved in
            _mc.mat file
 OUTPUT:
 x0: one parameter vector for which the model is stable.

 GRAPHS
 1) Pdf's of marginal distributions under the stability (dotted
     lines) and unstability (solid lines) regions
 2) Cumulative distributions of:
   - stable subset (dotted lines)
   - unacceptable subset (solid lines)
 3) Bivariate plots of significant correlation patterns
  ( abs(corrcoef) > alpha2) under the stable and unacceptable subsets

 USES qmc_sequence, stab_map_1, stab_map_2

 Written by Marco Ratto
 Joint Research Centre, The European Commission,
 (http://eemc.jrc.ec.europa.eu/),
 marco.ratto@jrc.it

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

0001 function x0 = stab_map_(OutputDirectoryName)
0002 %
0003 % function x0 = stab_map_(OutputDirectoryName)
0004 %
0005 % Mapping of stability regions in the prior ranges applying
0006 % Monte Carlo filtering techniques.
0007 %
0008 % INPUTS (from opt_gsa structure)
0009 % Nsam = MC sample size
0010 % fload = 0 to run new MC; 1 to load prevoiusly generated analysis
0011 % alpha2 =  significance level for bivariate sensitivity analysis
0012 % [abs(corrcoef) > alpha2]
0013 % prepSA = 1: save transition matrices for mapping reduced form
0014 %        = 0: no transition matrix saved (default)
0015 % pprior = 1: sample from prior ranges (default): sample saved in
0016 %            _prior.mat   file
0017 %        = 0: sample from posterior ranges: sample saved in
0018 %            _mc.mat file
0019 % OUTPUT:
0020 % x0: one parameter vector for which the model is stable.
0021 %
0022 % GRAPHS
0023 % 1) Pdf's of marginal distributions under the stability (dotted
0024 %     lines) and unstability (solid lines) regions
0025 % 2) Cumulative distributions of:
0026 %   - stable subset (dotted lines)
0027 %   - unacceptable subset (solid lines)
0028 % 3) Bivariate plots of significant correlation patterns
0029 %  ( abs(corrcoef) > alpha2) under the stable and unacceptable subsets
0030 %
0031 % USES qmc_sequence, stab_map_1, stab_map_2
0032 %
0033 % Written by Marco Ratto
0034 % Joint Research Centre, The European Commission,
0035 % (http://eemc.jrc.ec.europa.eu/),
0036 % marco.ratto@jrc.it
0037 %
0038 % Reference:
0039 % M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006.
0040 
0041 % Copyright (C) 2012 Dynare Team
0042 %
0043 % This file is part of Dynare.
0044 %
0045 % Dynare is free software: you can redistribute it and/or modify
0046 % it under the terms of the GNU General Public License as published by
0047 % the Free Software Foundation, either version 3 of the License, or
0048 % (at your option) any later version.
0049 %
0050 % Dynare is distributed in the hope that it will be useful,
0051 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0052 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0053 % GNU General Public License for more details.
0054 %
0055 % You should have received a copy of the GNU General Public License
0056 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0057 
0058 %global bayestopt_ estim_params_ dr_ options_ ys_ fname_
0059 global bayestopt_ estim_params_ options_ oo_ M_
0060 
0061 opt_gsa=options_.opt_gsa;
0062 
0063 Nsam   = opt_gsa.Nsam;
0064 fload  = opt_gsa.load_stab;
0065 ksstat = opt_gsa.ksstat;
0066 alpha2 = opt_gsa.alpha2_stab;
0067 pvalue_ks = opt_gsa.pvalue_ks;
0068 pvalue_corr = opt_gsa.pvalue_corr;
0069 prepSA = (opt_gsa.redform | opt_gsa.identification);
0070 pprior = opt_gsa.pprior;
0071 neighborhood_width = opt_gsa.neighborhood_width;
0072 ilptau = opt_gsa.ilptau;
0073 nliv   = opt_gsa.morris_nliv;
0074 ntra   = opt_gsa.morris_ntra;
0075 
0076 dr_ = oo_.dr;
0077 %if isfield(dr_,'ghx'),
0078 ys_ = oo_.dr.ys;
0079 nspred = dr_.nspred; %size(dr_.ghx,2);
0080 nboth = dr_.nboth;
0081 nfwrd = dr_.nfwrd;
0082 %end
0083 fname_ = M_.fname;
0084 
0085 np = estim_params_.np;
0086 nshock = estim_params_.nvx;
0087 nshock = nshock + estim_params_.nvn;
0088 nshock = nshock + estim_params_.ncx;
0089 nshock = nshock + estim_params_.ncn;
0090 lpmat0=[];
0091 xparam1=[];
0092 
0093 pshape = bayestopt_.pshape(nshock+1:end);
0094 p1 = bayestopt_.p1(nshock+1:end);
0095 p2 = bayestopt_.p2(nshock+1:end);
0096 p3 = bayestopt_.p3(nshock+1:end);
0097 p4 = bayestopt_.p4(nshock+1:end);
0098 
0099 if nargin==0,
0100     OutputDirectoryName='';
0101 end
0102 
0103 opt=options_;
0104 options_.periods=0;
0105 options_.nomoments=1;
0106 options_.irf=0;
0107 options_.noprint=1;
0108 options_.simul=0;
0109 if fload==0,
0110     %   if prepSA
0111     %     T=zeros(size(dr_.ghx,1),size(dr_.ghx,2)+size(dr_.ghu,2),Nsam/2);
0112     %   end
0113 
0114     if isfield(dr_,'ghx'),
0115         egg=zeros(length(dr_.eigval),Nsam);
0116     end
0117     yys=zeros(length(dr_.ys),Nsam);
0118 
0119     if opt_gsa.morris == 1
0120         [lpmat, OutFact] = Sampling_Function_2(nliv, np+nshock, ntra, ones(np+nshock, 1), zeros(np+nshock,1), []);
0121         lpmat = lpmat.*(nliv-1)/nliv+1/nliv/2;
0122         Nsam=size(lpmat,1);
0123         lpmat0 = lpmat(:,1:nshock);
0124         lpmat = lpmat(:,nshock+1:end);
0125     elseif opt_gsa.morris==3,
0126         lpmat = prep_ide(Nsam,np,5);
0127         Nsam=size(lpmat,1);
0128     else
0129         if np<52 & ilptau>0,
0130             [lpmat] = qmc_sequence(np, int64(0), 0, Nsam)';
0131             if np>30 | ilptau==2, % scrambled lptau
0132                 for j=1:np,
0133                     lpmat(:,j)=lpmat(randperm(Nsam),j);
0134                 end
0135             end
0136         else %ilptau==0
0137             %[lpmat] = rand(Nsam,np);
0138             for j=1:np,
0139                 lpmat(:,j) = randperm(Nsam)'./(Nsam+1); %latin hypercube
0140             end
0141 
0142         end
0143     end
0144     %   try
0145     dummy=prior_draw_gsa(1);
0146     %   catch
0147     %     if pprior,
0148     %       if opt_gsa.prior_range==0;
0149     %         error('Some unknown prior is specified or ML estimation,: use prior_range=1 option!!');
0150     %       end
0151     %     end
0152     %
0153     %   end
0154     if pprior,
0155         for j=1:nshock,
0156             if opt_gsa.morris~=1,
0157                 lpmat0(:,j) = randperm(Nsam)'./(Nsam+1); %latin hypercube
0158             end
0159             if opt_gsa.prior_range
0160                 lpmat0(:,j)=lpmat0(:,j).*(bayestopt_.ub(j)-bayestopt_.lb(j))+bayestopt_.lb(j);
0161             end
0162         end
0163         if opt_gsa.prior_range
0164             %       if opt_gsa.identification,
0165             %         deltx=min(0.001, 1/Nsam/2);
0166             %         for j=1:np,
0167             %           xdelt(:,:,j)=prior_draw_gsa(0,[lpmat0 lpmat]+deltx);
0168             %         end
0169             %       end
0170             for j=1:np,
0171                 lpmat(:,j)=lpmat(:,j).*(bayestopt_.ub(j+nshock)-bayestopt_.lb(j+nshock))+bayestopt_.lb(j+nshock);
0172             end
0173         else
0174             xx=prior_draw_gsa(0,[lpmat0 lpmat]);
0175             %       if opt_gsa.identification,
0176             %         deltx=min(0.001, 1/Nsam/2);
0177             %         ldum=[lpmat0 lpmat];
0178             %         ldum = prior_draw_gsa(0,ldum+deltx);
0179             %         for j=1:nshock+np,
0180             %           xdelt(:,:,j)=xx;
0181             %           xdelt(:,j,j)=ldum(:,j);
0182             %         end
0183             %         clear ldum
0184             %       end
0185             lpmat0=xx(:,1:nshock);
0186             lpmat=xx(:,nshock+1:end);
0187             clear xx;
0188         end
0189     else
0190         %         for j=1:nshock,
0191         %             xparam1(j) = oo_.posterior_mode.shocks_std.(bayestopt_.name{j});
0192         %             sd(j) = oo_.posterior_std.shocks_std.(bayestopt_.name{j});
0193         %             lpmat0(:,j) = randperm(Nsam)'./(Nsam+1); %latin hypercube
0194         %             lb = max(bayestopt_.lb(j), xparam1(j)-2*sd(j));
0195         %             ub1=xparam1(j)+(xparam1(j) - lb); % define symmetric range around the mode!
0196         %             ub = min(bayestopt_.ub(j),ub1);
0197         %             if ub<ub1,
0198         %                 lb=xparam1(j)-(ub-xparam1(j)); % define symmetric range around the mode!
0199         %             end
0200         %             lpmat0(:,j) = lpmat0(:,j).*(ub-lb)+lb;
0201         %         end
0202         %         %
0203         %         for j=1:np,
0204         %             xparam1(j+nshock) = oo_.posterior_mode.parameters.(bayestopt_.name{j+nshock});
0205         %             sd(j+nshock) = oo_.posterior_std.parameters.(bayestopt_.name{j+nshock});
0206         %             lb = max(bayestopt_.lb(j+nshock),xparam1(j+nshock)-2*sd(j+nshock));
0207         %             ub1=xparam1(j+nshock)+(xparam1(j+nshock) - lb); % define symmetric range around the mode!
0208         %             ub = min(bayestopt_.ub(j+nshock),ub1);
0209         %             if ub<ub1,
0210         %                 lb=xparam1(j+nshock)-(ub-xparam1(j+nshock)); % define symmetric range around the mode!
0211         %             end
0212         %             %ub = min(bayestopt_.ub(j+nshock),xparam1(j+nshock)+2*sd(j+nshock));
0213         %             if np>30 & np<52
0214         %                 lpmat(:,j) = lpmat(randperm(Nsam),j).*(ub-lb)+lb;
0215         %             else
0216         %                 lpmat(:,j) = lpmat(:,j).*(ub-lb)+lb;
0217         %             end
0218         %         end
0219         %load([fname_,'_mode'])
0220         eval(['load ' options_.mode_file '.mat;']);
0221         if neighborhood_width>0,
0222             for j=1:nshock,
0223                 lpmat0(:,j) = randperm(Nsam)'./(Nsam+1); %latin hypercube
0224                 ub=min([bayestopt_.ub(j) xparam1(j)*(1+neighborhood_width)]);
0225                 lb=max([bayestopt_.lb(j) xparam1(j)*(1-neighborhood_width)]);
0226                 lpmat0(:,j)=lpmat0(:,j).*(ub-lb)+lb;
0227             end
0228             for j=1:np,
0229                 ub=min([bayestopt_.ub(j+nshock) xparam1(j+nshock)*(1+neighborhood_width)]);
0230                 lb=max([bayestopt_.lb(j+nshock) xparam1(j+nshock)*(1-neighborhood_width)]);
0231                 lpmat(:,j)=lpmat(:,j).*(ub-lb)+lb;
0232             end
0233         else
0234             d = chol(inv(hh));
0235             lp=randn(Nsam*2,nshock+np)*d+kron(ones(Nsam*2,1),xparam1');
0236             for j=1:Nsam*2,
0237                 lnprior(j) = any(lp(j,:)'<=bayestopt_.lb | lp(j,:)'>=bayestopt_.ub);
0238             end
0239             ireal=[1:2*Nsam];
0240             ireal=ireal(find(lnprior==0));
0241             lp=lp(ireal,:);
0242             Nsam=min(Nsam, length(ireal));
0243             lpmat0=lp(1:Nsam,1:nshock);
0244             lpmat=lp(1:Nsam,nshock+1:end);
0245             clear lp lnprior ireal;
0246         end
0247     end
0248     %
0249     h = dyn_waitbar(0,'Please wait...');
0250     istable=[1:Nsam];
0251     jstab=0;
0252     iunstable=[1:Nsam];
0253     iindeterm=zeros(1,Nsam);
0254     iwrong=zeros(1,Nsam);
0255     for j=1:Nsam,
0256         M_.params(estim_params_.param_vals(:,1)) = lpmat(j,:)';
0257         %try stoch_simul([]);
0258         try
0259             [Tt,Rr,SteadyState,infox{j},M_,options_,oo_] = dynare_resolve(M_,options_,oo_,'restrict');
0260             if infox{j}(1)==0 && ~exist('T'),
0261                 dr_=oo_.dr;
0262                 T=zeros(size(dr_.ghx,1),size(dr_.ghx,2)+size(dr_.ghu,2),Nsam);
0263                 egg=zeros(length(dr_.eigval),Nsam);
0264             end
0265             if infox{j},
0266 %                 disp('no solution'),
0267                 if isfield(oo_.dr,'ghx'),
0268                     oo_.dr=rmfield(oo_.dr,'ghx');
0269                 end
0270             end
0271         catch
0272             if isfield(oo_.dr,'eigval'),
0273                 oo_.dr=rmfield(oo_.dr,'eigval');
0274             end
0275             if isfield(oo_.dr,'ghx'),
0276                 oo_.dr=rmfield(oo_.dr,'ghx');
0277             end
0278             disp('No solution could be found'),
0279         end
0280         dr_ = oo_.dr;
0281         if isfield(dr_,'ghx'),
0282             egg(:,j) = sort(dr_.eigval);
0283             iunstable(j)=0;
0284             if prepSA
0285                 jstab=jstab+1;
0286                 T(:,:,jstab) = [dr_.ghx dr_.ghu];
0287                 %         [A,B] = ghx2transition(squeeze(T(:,:,jstab)), ...
0288                 %           bayestopt_.restrict_var_list, ...
0289                 %           bayestopt_.restrict_columns, ...
0290                 %           bayestopt_.restrict_aux);
0291             end
0292             if ~exist('nspred'),
0293                 nspred = dr_.nspred; %size(dr_.ghx,2);
0294                 nboth = dr_.nboth;
0295                 nfwrd = dr_.nfwrd;
0296             end
0297         else
0298             istable(j)=0;
0299             if isfield(dr_,'eigval')
0300                 egg(:,j) = sort(dr_.eigval);
0301                 if exist('nspred')
0302                     if any(isnan(egg(1:nspred,j)))
0303                         iwrong(j)=j;
0304                     else
0305                         if (nboth | nfwrd) & abs(egg(nspred+1,j))<=options_.qz_criterium,
0306                             iindeterm(j)=j;
0307                         end
0308                     end
0309                 end
0310             else
0311                 if exist('egg'),
0312                     egg(:,j)=ones(size(egg,1),1).*NaN;
0313                 end
0314                 iwrong(j)=j;
0315             end
0316         end
0317         ys_=real(dr_.ys);
0318         yys(:,j) = ys_;
0319         ys_=yys(:,1);
0320         dyn_waitbar(j/Nsam,h,['MC iteration ',int2str(j),'/',int2str(Nsam)])
0321     end
0322     dyn_waitbar_close(h);
0323     if prepSA,
0324         T=T(:,:,1:jstab);
0325     end
0326     istable=istable(find(istable));  % stable params
0327     iunstable=iunstable(find(iunstable));   % unstable params
0328     iindeterm=iindeterm(find(iindeterm));  % indeterminacy
0329     iwrong=iwrong(find(iwrong));  % dynare could not find solution
0330 
0331     %     % map stable samples
0332     %     istable=[1:Nsam];
0333     %     for j=1:Nsam,
0334     %         if any(isnan(egg(1:nspred,j)))
0335     %             istable(j)=0;
0336     %         else
0337     %             if abs(egg(nspred,j))>=options_.qz_criterium; %(1-(options_.qz_criterium-1)); %1-1.e-5;
0338     %                 istable(j)=0;
0339     %                 %elseif (dr_.nboth | dr_.nfwrd) & abs(egg(nspred+1,j))<=options_.qz_criterium; %1+1.e-5;
0340     %             elseif (nboth | nfwrd) & abs(egg(nspred+1,j))<=options_.qz_criterium; %1+1.e-5;
0341     %                 istable(j)=0;
0342     %             end
0343     %         end
0344     %     end
0345     %     istable=istable(find(istable));  % stable params
0346     %
0347     %     % map unstable samples
0348     %     iunstable=[1:Nsam];
0349     %     for j=1:Nsam,
0350     %         %if abs(egg(dr_.npred+1,j))>1+1.e-5 & abs(egg(dr_.npred,j))<1-1.e-5;
0351     %         %if (dr_.nboth | dr_.nfwrd),
0352     %         if ~any(isnan(egg(1:5,j)))
0353     %             if (nboth | nfwrd),
0354     %                 if abs(egg(nspred+1,j))>options_.qz_criterium & abs(egg(nspred,j))<options_.qz_criterium; %(1-(options_.qz_criterium-1));
0355     %                     iunstable(j)=0;
0356     %                 end
0357     %             else
0358     %                 if abs(egg(nspred,j))<options_.qz_criterium; %(1-(options_.qz_criterium-1));
0359     %                     iunstable(j)=0;
0360     %                 end
0361     %             end
0362     %         end
0363     %     end
0364     %     iunstable=iunstable(find(iunstable));   % unstable params
0365     bkpprior.pshape=bayestopt_.pshape;
0366     bkpprior.p1=bayestopt_.p1;
0367     bkpprior.p2=bayestopt_.p2;
0368     bkpprior.p3=bayestopt_.p3;
0369     bkpprior.p4=bayestopt_.p4;
0370     if pprior,
0371         if ~prepSA
0372             save([OutputDirectoryName '/' fname_ '_prior.mat'], ...
0373                 'bkpprior','lpmat','lpmat0','iunstable','istable','iindeterm','iwrong', ...
0374                 'egg','yys','nspred','nboth','nfwrd')
0375         else
0376             save([OutputDirectoryName '/' fname_ '_prior.mat'], ...
0377                 'bkpprior','lpmat','lpmat0','iunstable','istable','iindeterm','iwrong', ...
0378                 'egg','yys','T','nspred','nboth','nfwrd')
0379         end
0380 
0381     else
0382         if ~prepSA
0383             save([OutputDirectoryName '/' fname_ '_mc.mat'], ...
0384                 'lpmat','lpmat0','iunstable','istable','iindeterm','iwrong', ...
0385                 'egg','yys','nspred','nboth','nfwrd')
0386         else
0387             save([OutputDirectoryName '/' fname_ '_mc.mat'], ...
0388                 'lpmat','lpmat0','iunstable','istable','iindeterm','iwrong', ...
0389                 'egg','yys','T','nspred','nboth','nfwrd')
0390         end
0391     end
0392 else
0393     if pprior,
0394         filetoload=[OutputDirectoryName '/' fname_ '_prior.mat'];
0395     else
0396         filetoload=[OutputDirectoryName '/' fname_ '_mc.mat'];
0397     end
0398     load(filetoload,'lpmat','lpmat0','iunstable','istable','iindeterm','iwrong','egg','yys','nspred','nboth','nfwrd')
0399     Nsam = size(lpmat,1);
0400     if pprior==0,
0401         eval(['load ' options_.mode_file '.mat;']);
0402     end
0403 
0404 
0405     if prepSA & isempty(strmatch('T',who('-file', filetoload),'exact')),
0406         h = dyn_waitbar(0,'Please wait...');
0407         options_.periods=0;
0408         options_.nomoments=1;
0409         options_.irf=0;
0410         options_.noprint=1;
0411         stoch_simul([]);
0412         %T=zeros(size(dr_.ghx,1),size(dr_.ghx,2)+size(dr_.ghu,2),length(istable));
0413         ntrans=length(istable);
0414         for j=1:ntrans,
0415             M_.params(estim_params_.param_vals(:,1)) = lpmat(istable(j),:)';
0416             %stoch_simul([]);
0417             [Tt,Rr,SteadyState,info,M_,options_,oo_] = dynare_resolve(M_,options_,oo_,'restrict');
0418             % This syntax is not compatible with the current version of dynare_resolve [stepan].
0419             %[Tt,Rr,SteadyState,info] = dynare_resolve(bayestopt_.restrict_var_list,...
0420             %    bayestopt_.restrict_columns,...
0421             %    bayestopt_.restrict_aux);
0422             if ~exist('T')
0423                 T=zeros(size(dr_.ghx,1),size(dr_.ghx,2)+size(dr_.ghu,2),ntrans);
0424             end
0425             dr_ = oo_.dr;
0426             T(:,:,j) = [dr_.ghx dr_.ghu];
0427             if ~exist('nspred')
0428                 nspred = dr_.nspred; %size(dr_.ghx,2);
0429                 nboth = dr_.nboth;
0430                 nfwrd = dr_.nfwrd;
0431             end
0432             ys_=real(dr_.ys);
0433             yys(:,j) = ys_;
0434             ys_=yys(:,1);
0435             dyn_waitbar(j/ntrans,h,['MC iteration ',int2str(j),'/',int2str(ntrans)])
0436         end
0437         dyn_waitbar_close(h);
0438         save(filetoload,'T','-append')
0439     elseif prepSA
0440         load(filetoload,'T')
0441     end
0442 end
0443 
0444 if pprior
0445     aname='prior_stab';
0446     auname='prior_unacceptable';
0447     aunstname='prior_unstable';
0448     aindname='prior_indeterm';
0449     awrongname='prior_wrong';
0450     asname='prior_stable';
0451 else
0452     aname='mc_stab';
0453     auname='mc_unacceptable';
0454     aunstname='mc_unstable';
0455     aindname='mc_indeterm';
0456     awrongname='mc_wrong';
0457     asname='mc_stable';
0458 end
0459 delete([OutputDirectoryName,'/',fname_,'_',aname,'_*.*']);
0460 %delete([OutputDirectoryName,'/',fname_,'_',aname,'_SA_*.*']);
0461 delete([OutputDirectoryName,'/',fname_,'_',asname,'_corr_*.*']);
0462 delete([OutputDirectoryName,'/',fname_,'_',auname,'_corr_*.*']);
0463 delete([OutputDirectoryName,'/',fname_,'_',aunstname,'_corr_*.*']);
0464 delete([OutputDirectoryName,'/',fname_,'_',aindname,'_corr_*.*']);
0465 
0466 if length(iunstable)>0 & length(iunstable)<Nsam,
0467     fprintf(['%4.1f%% of the prior support is stable.\n'],length(istable)/Nsam*100)
0468     fprintf(['%4.1f%% of the prior support is unstable.\n'],(length(iunstable)-length(iwrong)-length(iindeterm) )/Nsam*100)
0469     if ~isempty(iindeterm),
0470         fprintf(['%4.1f%% of the prior support gives indeterminacy.'],length(iindeterm)/Nsam*100)
0471     end
0472     if ~isempty(iwrong),
0473         disp(' ');
0474         disp(['For ',num2str(length(iwrong)/Nsam*100,'%1.3f'),'\% of the prior support dynare could not find a solution.'])
0475     end
0476     disp(' ');
0477     % Blanchard Kahn
0478     [proba, dproba] = stab_map_1(lpmat, istable, iunstable, aname,0);
0479 %     indstab=find(dproba>ksstat);
0480     indstab=find(proba<pvalue_ks);
0481     disp('Smirnov statistics in driving acceptable behaviour')
0482     for j=1:length(indstab),
0483         disp([M_.param_names(estim_params_.param_vals(indstab(j),1),:),'   d-stat = ', num2str(dproba(indstab(j)),'%1.3f'),'   p-value = ', num2str(proba(indstab(j)),'%1.3f')])
0484     end
0485     disp(' ');
0486     if ~isempty(indstab)
0487         stab_map_1(lpmat, istable, iunstable, aname, 1, indstab, OutputDirectoryName);
0488     end
0489     ixun=iunstable(find(~ismember(iunstable,[iindeterm,iwrong])));
0490     if ~isempty(iindeterm),
0491         [proba, dproba] = stab_map_1(lpmat, [1:Nsam], iindeterm, [aname, '_indet'],0);
0492 %         indindet=find(dproba>ksstat);
0493         indindet=find(proba<pvalue_ks);
0494         disp('Smirnov statistics in driving indeterminacy')
0495         for j=1:length(indindet),
0496             disp([M_.param_names(estim_params_.param_vals(indindet(j),1),:),'   d-stat = ', num2str(dproba(indindet(j)),'%1.3f'),'   p-value = ', num2str(proba(indindet(j)),'%1.3f')])
0497         end
0498         disp(' ');
0499         if ~isempty(indindet)
0500             stab_map_1(lpmat, [1:Nsam], iindeterm, [aname, '_indet'], 1, indindet, OutputDirectoryName);
0501         end
0502     end
0503 
0504     if ~isempty(ixun),
0505         [proba, dproba] = stab_map_1(lpmat, [1:Nsam], ixun, [aname, '_unst'],0);
0506 %         indunst=find(dproba>ksstat);
0507         indunst=find(proba<pvalue_ks);
0508         disp('Smirnov statistics in driving instability')
0509         for j=1:length(indunst),
0510             disp([M_.param_names(estim_params_.param_vals(indunst(j),1),:),'   d-stat = ', num2str(dproba(indunst(j)),'%1.3f'),'   p-value = ', num2str(proba(indunst(j)),'%1.3f')])
0511         end
0512         disp(' ');
0513         if ~isempty(indunst)
0514             stab_map_1(lpmat, [1:Nsam], ixun, [aname, '_unst'], 1, indunst, OutputDirectoryName);
0515         end
0516     end
0517 
0518     if ~isempty(iwrong),
0519         [proba, dproba] = stab_map_1(lpmat, [1:Nsam], iwrong, [aname, '_wrong'],0);
0520 %         indwrong=find(dproba>ksstat);
0521         indwrong=find(proba<pvalue_ks);
0522         disp('Smirnov statistics in driving no solution')
0523         for j=1:length(indwrong),
0524             disp([M_.param_names(estim_params_.param_vals(indwrong(j),1),:),'   d-stat = ', num2str(dproba(indwrong(j)),'%1.3f'),'   p-value = ', num2str(proba(indwrong(j)),'%1.3f')])
0525         end
0526         disp(' ');
0527         if ~isempty(indwrong)
0528             stab_map_1(lpmat, [1:Nsam], iwrong, [aname, '_wrong'], 1, indwrong, OutputDirectoryName);
0529         end
0530     end
0531 
0532     disp(' ')
0533     disp('Starting bivariate analysis:')
0534 
0535     c0=corrcoef(lpmat(istable,:));
0536     c00=tril(c0,-1);
0537 
0538     stab_map_2(lpmat(istable,:),alpha2, pvalue_corr, asname, OutputDirectoryName,xparam1);
0539     if length(iunstable)>10,
0540         stab_map_2(lpmat(iunstable,:),alpha2, pvalue_corr, auname, OutputDirectoryName,xparam1);
0541     end
0542     if length(iindeterm)>10,
0543         stab_map_2(lpmat(iindeterm,:),alpha2, pvalue_corr, aindname, OutputDirectoryName,xparam1);
0544     end
0545     if length(ixun)>10,
0546         stab_map_2(lpmat(ixun,:),alpha2, pvalue_corr, aunstname, OutputDirectoryName,xparam1);
0547     end
0548     if length(iwrong)>10,
0549         stab_map_2(lpmat(iwrong,:),alpha2, pvalue_corr, awrongname, OutputDirectoryName,xparam1);
0550     end
0551 
0552     x0=0.5.*(bayestopt_.ub(1:nshock)-bayestopt_.lb(1:nshock))+bayestopt_.lb(1:nshock);
0553     x0 = [x0; lpmat(istable(1),:)'];
0554     if istable(end)~=Nsam
0555         M_.params(estim_params_.param_vals(:,1)) = lpmat(istable(1),:)';
0556         [oo_.dr,info,M_,options_,oo_] = resol(0,M_,options_,oo_);
0557         %     stoch_simul([]);
0558     end
0559 else
0560     if length(iunstable)==0,
0561         disp('All parameter values in the specified ranges are stable!')
0562         x0=0.5.*(bayestopt_.ub(1:nshock)-bayestopt_.lb(1:nshock))+bayestopt_.lb(1:nshock);
0563         x0 = [x0; lpmat(istable(1),:)'];
0564     else
0565         disp('All parameter values in the specified ranges are not acceptable!')
0566         x0=[];
0567     end
0568 
0569 end
0570 
0571 
0572 options_.periods=opt.periods;
0573 if isfield(opt,'nomoments'),
0574     options_.nomoments=opt.nomoments;
0575 end
0576 options_.irf=opt.irf;
0577 options_.noprint=opt.noprint;
0578 if isfield(opt,'simul'),
0579     options_.simul=opt.simul;
0580 end
0581 
0582 
0583