dynare_identification.m 23.4 KB
Newer Older
ratto's avatar
ratto committed
1
function [pdraws, TAU, GAM, LRE, gp, H, JJ] = dynare_identification(options_ident, pdraws0)
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
%function [pdraws, TAU, GAM, LRE, gp, H, JJ] = dynare_identification(options_ident, pdraws0)
%
% INPUTS
%    o options_ident    [structure] identification options
%    o pdraws0          [matrix] optional: matrix of MC sample of model params. 
%    
% OUTPUTS
%    o pdraws           [matrix] matrix of MC sample of model params used
%    o TAU,             [matrix] MC sample of entries in the model solution (stacked vertically)
%    o GAM,             [matrix] MC sample of entries in the moments (stacked vertically)
%    o LRE,             [matrix] MC sample of entries in LRE model (stacked vertically)
%    o gp,              [matrix] derivatives of the Jacobian (LRE model)
%    o H,               [matrix] derivatives of the model solution
%    o JJ               [matrix] derivatives of the  moments
%    
% SPECIAL REQUIREMENTS
%    None
19
20
21

% main 
%
Sébastien Villemot's avatar
Sébastien Villemot committed
22
% Copyright (C) 2010-2012 Dynare Team
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.

global M_ options_ oo_ bayestopt_ estim_params_

Marco Ratto's avatar
Marco Ratto committed
41
if exist('OCTAVE_VERSION')
42
    warning('off'),
Marco Ratto's avatar
Marco Ratto committed
43
else
44
    warning off,
Marco Ratto's avatar
Marco Ratto committed
45
46
end

47
fname_ = M_.fname;
48
49
50
51
52
if ~isfield(M_,'dname'),
    M_.dname = M_.fname;
end
options_ident = set_default_option(options_ident,'gsa_sample_file',0);
options_ident = set_default_option(options_ident,'parameter_set','prior_mean');
53
options_ident = set_default_option(options_ident,'load_ident_files',0);
54
options_ident = set_default_option(options_ident,'useautocorr',0);
55
options_ident = set_default_option(options_ident,'ar',1);
56
options_ident = set_default_option(options_ident,'prior_mc',1);
57
options_ident = set_default_option(options_ident,'prior_range',0);
58
59
60
options_ident = set_default_option(options_ident,'periods',300);
options_ident = set_default_option(options_ident,'replic',100);
options_ident = set_default_option(options_ident,'advanced',0);
61
options_ident = set_default_option(options_ident,'normalize_jacobians',1);
62
options_ident = set_default_option(options_ident,'lik_init',1);
63
options_ident = set_default_option(options_ident,'analytic_derivation',1);
64
65
66
67
68
69
70
71
72
73
if isfield(options_ident,'nograph'),
    options_.nograph=options_ident.nograph;
end
if isfield(options_ident,'nodisplay'),
    options_.nodisplay=options_ident.nodisplay;
end
if isfield(options_ident,'graph_format'),
    options_.graph_format=options_ident.graph_format;
end

74
if options_ident.gsa_sample_file,
Marco Ratto's avatar
Marco Ratto committed
75
    GSAFolder = checkpath('gsa',M_.dname);
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
    if options_ident.gsa_sample_file==1,
        load([GSAFolder,filesep,fname_,'_prior'],'lpmat','lpmat0','istable');
    elseif options_ident.gsa_sample_file==2,
        load([GSAFolder,filesep,fname_,'_mc'],'lpmat','lpmat0','istable');
    else
        load([GSAFolder,filesep,options_ident.gsa_sample_file],'lpmat','lpmat0','istable');
    end
    if isempty(istable),
        istable=1:size(lpmat,1);
    end
    if ~isempty(lpmat0),
        lpmatx=lpmat0(istable,:);
    else
        lpmatx=[];
    end
    pdraws0 = [lpmatx lpmat(istable,:)];
    clear lpmat lpmat0 istable;
93
elseif nargin==1,
94
95
96
97
    pdraws0=[];
end
external_sample=0;
if nargin==2 || ~isempty(pdraws0),
98
    options_ident.prior_mc=size(pdraws0,1);
99
100
    options_ident.load_ident_files = 0;
    external_sample=1;
101
end
102
103
104
105
106
107
if isempty(estim_params_),
    options_ident.prior_mc=1;
    options_ident.prior_range=0;
    prior_exist=0;
else
    prior_exist=1;
108
    parameters = options_ident.parameter_set;
109
end
110

Marco Ratto's avatar
Marco Ratto committed
111
% options_ident.load_ident_files=1;
112
iload = options_ident.load_ident_files;
Marco Ratto's avatar
Marco Ratto committed
113
%options_ident.advanced=1;
114
advanced = options_ident.advanced;
115
nlags = options_ident.ar;
116
117
periods = options_ident.periods;
replic = options_ident.replic;
118
useautocorr = options_ident.useautocorr;
Marco Ratto's avatar
Marco Ratto committed
119
options_.order=1;
120
121
122
options_.ar=nlags;
options_.prior_mc = options_ident.prior_mc;
options_.options_ident = options_ident;
Marco Ratto's avatar
Marco Ratto committed
123
options_.Schur_vec_tol = 1.e-8;
124
options_.nomoments=0;
Marco Ratto's avatar
Marco Ratto committed
125
options_.analytic_derivation=1;
126

127
options_ = set_default_option(options_,'datafile','');
128
options_.mode_compute = 0;
129
options_.plot_priors = 0;
130
options_.smoother=1;
131
[dataset_,xparam1, M_, options_, oo_, estim_params_,bayestopt_]=dynare_estimation_init(M_.endo_names,fname_,1, M_, options_, oo_, estim_params_, bayestopt_);
132
options_ident.analytic_derivation_mode = options_.analytic_derivation_mode;
133
134
135
136
137
138
if isempty(dataset_),
    dataset_.info.ntobs = periods;
    dataset_.info.nvobs = rows(options_.varobs);
    dataset_.info.varobs = options_.varobs;
    dataset_.rawdata = [];
    dataset_.missing.state = 0;
139
140
141
142
    for jdata=1:periods,
        temp1{jdata}=[1:dataset_.info.nvobs]';
    end
    dataset_.missing.aindex = temp1;
143
144
    dataset_.missing.vindex = [];
    dataset_.missing.number_of_observations = [];
145
    dataset_.missing.no_more_missing_observations = 1;
146
147
148
149
150
151
152
153
154
    dataset_.descriptive.mean = [];
    dataset_.data = [];

%     data_info.gend = periods;
%     data_info.data = [];
%     data_info.data_index = [];
%     data_info.number_of_observations = periods*size(options_.varobs,1);
%     data_info.no_more_missing_observations = 0;
%     data_info.missing_value = 0;
155
end
156
157
158

% results = prior_sampler(0,M_,bayestopt_,options_,oo_);

159
if prior_exist
160
    if any(bayestopt_.pshape > 0)
161
162
163
164
165
        if options_ident.prior_range
            prior_draw(1,1);
        else
            prior_draw(1);
        end
166
    else
167
        options_ident.prior_mc=1;
168
    end
169
end
170

171
172
SampleSize = options_ident.prior_mc;

Marco Ratto's avatar
Marco Ratto committed
173
if ~(exist('sylvester3','file')==2),
174
175
176
177
178

    dynareroot = strrep(which('dynare'),'dynare.m','');
    addpath([dynareroot 'gensylv'])
end

179
IdentifDirectoryName = CheckPath('identification',M_.dname);
180
if prior_exist,
181

182
183
184
185
    indx = [];
    if ~isempty(estim_params_.param_vals),
        indx = estim_params_.param_vals(:,1);
    end
186
187
188
189
190
191
192
193
    indexo=[];
    if ~isempty(estim_params_.var_exo)
        indexo = estim_params_.var_exo(:,1);
    end

    nparam = length(bayestopt_.name);
    np = estim_params_.np;
    name = bayestopt_.name;
194
    name_tex = char(M_.exo_names_tex(indexo,:),M_.param_names_tex(indx,:));
195

196
197
198
199
200
201
202
203
204
205
206
    offset = estim_params_.nvx;
    offset = offset + estim_params_.nvn;
    offset = offset + estim_params_.ncx;
    offset = offset + estim_params_.ncn;
else
    indx = [1:M_.param_nbr];
    indexo = [1:M_.exo_nbr];
    offset = M_.exo_nbr;
    np = M_.param_nbr;
    nparam = np+offset;
    name = [cellstr(M_.exo_names); cellstr(M_.param_names)];
207
    name_tex = [cellstr(M_.exo_names_tex); cellstr(M_.param_names_tex)];
208
end
209

Marco Ratto's avatar
Marco Ratto committed
210
options_ident = set_default_option(options_ident,'max_dim_cova_group',min([2,nparam-1]));
211
options_ident.max_dim_cova_group = min([options_ident.max_dim_cova_group,nparam-1]);
Marco Ratto's avatar
Marco Ratto committed
212
213


214
MaxNumberOfBytes=options_.MaxNumberOfBytes;
215
store_options_ident = options_ident;
216
217
218
disp(' ')
disp(['==== Identification analysis ====' ]),
disp(' ')
219

220
if iload <=0,
221
    
Marco Ratto's avatar
Marco Ratto committed
222
223
    [I,J]=find(M_.lead_lag_incidence');
    if prior_exist,
224
225
226
227
228
229
230
231
232
233
234
235
%         if exist([fname_,'_mean.mat'],'file'),
% %             disp('Testing posterior mean')
%             load([fname_,'_mean'],'xparam1')
%             pmean = xparam1';
%             clear xparam1
%         end
%         if exist([fname_,'_mode.mat'],'file'),
% %             disp('Testing posterior mode')
%             load([fname_,'_mode'],'xparam1')
%             pmode = xparam1';
%             clear xparam1
%         end
Marco Ratto's avatar
Marco Ratto committed
236
        params = set_prior(estim_params_,M_,options_)';
237
238
239
        if all(bayestopt_.pshape == 0)
            parameters = 'ML_Starting_value';
            disp('Testing ML Starting value')
240
        else
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
        switch parameters
            case 'posterior_mode'
                disp('Testing posterior mode')
                params(1,:) = get_posterior_parameters('mode');
            case 'posterior_mean'
                disp('Testing posterior mean')
                params(1,:) = get_posterior_parameters('mean');
            case 'posterior_median'
                disp('Testing posterior median')
                params(1,:) = get_posterior_parameters('median');
            case 'prior_mode'
                disp('Testing prior mode')
                params(1,:) = bayestopt_.p5(:);
            case 'prior_mean'
                disp('Testing prior mean')
                params(1,:) = bayestopt_.p1;
            otherwise
                disp('The option parameter_set has to be equal to:')
                disp('                   ''posterior_mode'', ')
                disp('                   ''posterior_mean'', ')
                disp('                   ''posterior_median'', ')
                disp('                   ''prior_mode'' or')
                disp('                   ''prior_mean''.')
                error;
        end
266
        end
Marco Ratto's avatar
Marco Ratto committed
267
268
    else
        params = [sqrt(diag(M_.Sigma_e))', M_.params'];
269
        parameters = 'Current_params';
270
        disp('Testing current parameter values')
Marco Ratto's avatar
Marco Ratto committed
271
    end
272
    [idehess_point, idemoments_point, idemodel_point, idelre_point, derivatives_info_point, info] = ...
273
        identification_analysis(params,indx,indexo,options_ident,dataset_, prior_exist, name_tex,1);
274
275
276
277
278
    if info(1)~=0,
        disp(' ')
        disp('----------- ')
        disp('Parameter error:')
        disp(['The model does not solve for ', parameters, ' with error code info = ', int2str(info(1))]),
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
        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
307
308
        disp('----------- ')
        disp(' ')
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
        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
330
    else
331
332
333
334
    idehess_point.params=params;
%     siH = idemodel_point.siH;
%     siJ = idemoments_point.siJ;
%     siLRE = idelre_point.siLRE;
Marco Ratto's avatar
Marco Ratto committed
335
336
337
%     normH = max(abs(siH)')';
%     normJ = max(abs(siJ)')';
%     normLRE = max(abs(siLRE)')';
338
    save([IdentifDirectoryName '/' M_.fname '_identif.mat'], 'idehess_point', 'idemoments_point','idemodel_point', 'idelre_point','store_options_ident')
339
    disp_identification(params, idemodel_point, idemoments_point, name, advanced);
340
341
342
    if ~options_.nograph,
        plot_identification(params,idemoments_point,idehess_point,idemodel_point,idelre_point,advanced,parameters,name,IdentifDirectoryName);
    end
343
    end
344

345
    if SampleSize > 1,
Marco Ratto's avatar
Marco Ratto committed
346
347
        disp(' ')
        disp('Monte Carlo Testing')
348
        h = dyn_waitbar(0,'Monte Carlo identification checks ...');
Marco Ratto's avatar
Marco Ratto committed
349
350
351
352
353
354
355
356
357
        iteration = 0;
        loop_indx = 0;
        file_index = 0;
        run_index = 0;
        options_MC=options_ident;
        options_MC.advanced=0;
    else
        iteration = 1;
        pdraws = [];
358
    end
359
360
    while iteration < SampleSize,
        loop_indx = loop_indx+1;
361
        if external_sample,
Marco Ratto's avatar
Marco Ratto committed
362
            params = pdraws0(iteration+1,:);
363
        else
Marco Ratto's avatar
Marco Ratto committed
364
365
            params = prior_draw();
        end
366
        [dum1, ideJ, ideH, ideGP, dum2 , info] = ...
367
            identification_analysis(params,indx,indexo,options_MC,dataset_, prior_exist, name_tex,0);
368
        if iteration==0 && info(1)==0,
Marco Ratto's avatar
Marco Ratto committed
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
            MAX_tau   = min(SampleSize,ceil(MaxNumberOfBytes/(size(ideH.siH,1)*nparam)/8));
            stoH = zeros([size(ideH.siH,1),nparam,MAX_tau]);
            stoJJ = zeros([size(ideJ.siJ,1),nparam,MAX_tau]);
            stoLRE = zeros([size(ideGP.siLRE,1),np,MAX_tau]);
            TAU = zeros(size(ideH.siH,1),SampleSize);
            GAM = zeros(size(ideJ.siJ,1),SampleSize);
            LRE = zeros(size(ideGP.siLRE,1),SampleSize);
            pdraws = zeros(SampleSize,nparam);
            idemoments.indJJ = ideJ.indJJ;
            idemodel.indH = ideH.indH;
            idelre.indLRE = ideGP.indLRE;
            idemoments.ind0 = zeros(SampleSize,nparam);
            idemodel.ind0 = zeros(SampleSize,nparam);
            idelre.ind0 = zeros(SampleSize,np);
            idemoments.jweak = zeros(SampleSize,nparam);
            idemodel.jweak = zeros(SampleSize,nparam);
            idelre.jweak = zeros(SampleSize,np);
            idemoments.jweak_pair = zeros(SampleSize,nparam*(nparam+1)/2);
            idemodel.jweak_pair = zeros(SampleSize,nparam*(nparam+1)/2);
            idelre.jweak_pair = zeros(SampleSize,np*(np+1)/2);
            idemoments.cond = zeros(SampleSize,1);
            idemodel.cond = zeros(SampleSize,1);
            idelre.cond = zeros(SampleSize,1);
            idemoments.Mco = zeros(SampleSize,nparam);
            idemodel.Mco = zeros(SampleSize,nparam);
            idelre.Mco = zeros(SampleSize,np);
395
396
            idemoments.S = zeros(SampleSize,min(8,nparam));
            idemoments.V = zeros(SampleSize,nparam,min(8,nparam));
Marco Ratto's avatar
Marco Ratto committed
397
            delete([IdentifDirectoryName '/' M_.fname '_identif_*.mat'])
398
        end
399
        if info(1)==0,
400
401
            iteration = iteration + 1;
            run_index = run_index + 1;
Marco Ratto's avatar
Marco Ratto committed
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
            TAU(:,iteration)=ideH.TAU;
            LRE(:,iteration)=ideGP.LRE;
            GAM(:,iteration)=ideJ.GAM;
            idemoments.cond(iteration,1)=ideJ.cond;
            idemodel.cond(iteration,1)=ideH.cond;
            idelre.cond(iteration,1)=ideGP.cond;
            idemoments.ino(iteration,1)=ideJ.ino;
            idemodel.ino(iteration,1)=ideH.ino;
            idelre.ino(iteration,1)=ideGP.ino;
            idemoments.ind0(iteration,:)=ideJ.ind0;
            idemodel.ind0(iteration,:)=ideH.ind0;
            idelre.ind0(iteration,:)=ideGP.ind0;
            idemoments.jweak(iteration,:)=ideJ.jweak;
            idemodel.jweak(iteration,:)=ideH.jweak;
            idelre.jweak(iteration,:)=ideGP.jweak;
            idemoments.jweak_pair(iteration,:)=ideJ.jweak_pair;
            idemodel.jweak_pair(iteration,:)=ideH.jweak_pair;
            idelre.jweak_pair(iteration,:)=ideGP.jweak_pair;
            idemoments.Mco(iteration,:)=ideJ.Mco;
            idemodel.Mco(iteration,:)=ideH.Mco;
            idelre.Mco(iteration,:)=ideGP.Mco;
423
424
            idemoments.S(iteration,:)=ideJ.S;
            idemoments.V(iteration,:,:)=ideJ.V;
Marco Ratto's avatar
Marco Ratto committed
425
426
427
428
429
430
431
432
433
434
            stoLRE(:,:,run_index) = ideGP.siLRE;
            stoH(:,:,run_index) = ideH.siH;
            stoJJ(:,:,run_index) = ideJ.siJ;
            pdraws(iteration,:) = params;
            if run_index==MAX_tau || iteration==SampleSize,
                file_index = file_index + 1;
                if run_index<MAX_tau,
                    stoH = stoH(:,:,1:run_index);
                    stoJJ = stoJJ(:,:,1:run_index);
                    stoLRE = stoLRE(:,:,1:run_index);
435
                end
Marco Ratto's avatar
Marco Ratto committed
436
437
438
439
440
                save([IdentifDirectoryName '/' M_.fname '_identif_' int2str(file_index) '.mat'], 'stoH', 'stoJJ', 'stoLRE')
                run_index = 0;
                stoH = zeros(size(stoH));
                stoJJ = zeros(size(stoJJ));
                stoLRE = zeros(size(stoLRE));
441
                
Marco Ratto's avatar
Marco Ratto committed
442
443
444
            end
            
            if SampleSize > 1,
445
446
447
448
449
%                 if exist('OCTAVE_VERSION') || options_.console_mode,
%                     console_waitbar(0,iteration/SampleSize);
%                 else
                    dyn_waitbar(iteration/SampleSize,h,['MC identification checks ',int2str(iteration),'/',int2str(SampleSize)])
%                 end
450
451
            end
        end
Marco Ratto's avatar
Marco Ratto committed
452
        
453
    end
454
455
    
    
456
    if SampleSize > 1,
457
458
459
460
461
462
        if exist('OCTAVE_VERSION') || options_.console_mode,
            fprintf('\n');
            diary on;
        else
            close(h),
        end
Marco Ratto's avatar
Marco Ratto committed
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
        normTAU=std(TAU')';
        normLRE=std(LRE')';
        normGAM=std(GAM')';
        normaliz1=std(pdraws);
        iter=0;
        for ifile_index = 1:file_index,
            load([IdentifDirectoryName '/' M_.fname '_identif_' int2str(ifile_index) '.mat'], 'stoH', 'stoJJ', 'stoLRE')
            for irun=1:size(stoH,3),
                iter=iter+1;
                siJnorm(iter,:) = vnorm(stoJJ(:,:,irun)./repmat(normGAM,1,nparam)).*normaliz1;
                siHnorm(iter,:) = vnorm(stoH(:,:,irun)./repmat(normTAU,1,nparam)).*normaliz1;
                siLREnorm(iter,:) = vnorm(stoLRE(:,:,irun)./repmat(normLRE,1,nparam-offset)).*normaliz1(offset+1:end);
            end
            
        end
Marco Ratto's avatar
Marco Ratto committed
478
479
480
481
482
483
        idemoments.siJnorm = siJnorm;
        idemodel.siHnorm = siHnorm;
        idelre.siLREnorm = siLREnorm;
        save([IdentifDirectoryName '/' M_.fname '_identif.mat'], 'pdraws', 'idemodel', 'idemoments', 'idelre', ... %'indJJ', 'indH', 'indLRE', ...
            'TAU', 'GAM', 'LRE','-append')
    else
484
485
486
        siJnorm = idemoments_point.siJnorm;
        siHnorm = idemodel_point.siHnorm;
        siLREnorm = idelre_point.siLREnorm;
487
    end
488
    
489
else
Marco Ratto's avatar
Marco Ratto committed
490
    load([IdentifDirectoryName '/' M_.fname '_identif'])
491
492
493
%     identFiles = dir([IdentifDirectoryName '/' M_.fname '_identif_*']);
    parameters = store_options_ident.parameter_set;
    options_ident.parameter_set = parameters;
494
    options_ident.prior_mc=size(pdraws,1);
495
    SampleSize = options_ident.prior_mc;
496
    options_.options_ident = options_ident;
ratto's avatar
ratto committed
497
end  
498

499
if nargout>3 && iload,
500
501
502
    filnam = dir([IdentifDirectoryName '/' M_.fname '_identif_*.mat']);
    H=[];
    JJ = [];
503
    gp = [];
504
505
506
507
    for j=1:length(filnam),
        load([IdentifDirectoryName '/' M_.fname '_identif_',int2str(j),'.mat']);
        H = cat(3,H, stoH(:,abs(iload),:));
        JJ = cat(3,JJ, stoJJ(:,abs(iload),:));
ratto's avatar
ratto committed
508
        gp = cat(3,gp, stoLRE(:,abs(iload),:));
509
        
510
511
512
    end
end

Marco Ratto's avatar
Marco Ratto committed
513
if iload,
514
    disp(['Testing ',parameters])
515
    disp_identification(idehess_point.params, idemodel_point, idemoments_point, name,advanced);
516
517
518
    if ~options_.nograph,
        plot_identification(idehess_point.params,idemoments_point,idehess_point,idemodel_point,idelre_point,advanced,parameters,name,IdentifDirectoryName);
    end
519
end
520
if SampleSize > 1,
521
    fprintf('\n')
Marco Ratto's avatar
Marco Ratto committed
522
523
    disp('Testing MC sample')
    disp_identification(pdraws, idemodel, idemoments, name);
524
525
526
    if ~options_.nograph,
        plot_identification(pdraws,idemoments,idehess_point,idemodel,idelre,advanced,'MC sample - ',name, IdentifDirectoryName);
    end
Marco Ratto's avatar
Marco Ratto committed
527
528
    if advanced,
        jcrit=find(idemoments.ino);
529
530
        if length(jcrit)<SampleSize,
            if isempty(jcrit),
531
                [dum,jmax]=max(idemoments.cond);
532
533
534
                fprintf('\n')
                tittxt = 'Draw with HIGHEST condition number';
                fprintf('\n')
535
                disp(['Testing ',tittxt, '. Press ENTER']), pause(5),
536
537
                if ~iload,
                    [idehess_max, idemoments_max, idemodel_max, idelre_max, derivatives_info_max] = ...
538
                        identification_analysis(pdraws(jmax,:),indx,indexo,options_ident,dataset_, prior_exist, name_tex,1);
539
540
                    save([IdentifDirectoryName '/' M_.fname '_identif.mat'], 'idehess_max', 'idemoments_max','idemodel_max', 'idelre_max', 'jmax', '-append');
                end
541
                disp_identification(pdraws(jmax,:), idemodel_max, idemoments_max, name,1);
542
                close all,
543
544
545
                if ~options_.nograph,
                    plot_identification(pdraws(jmax,:),idemoments_max,idehess_max,idemodel_max,idelre_max,1,tittxt,name,IdentifDirectoryName);
                end
546
                [dum,jmin]=min(idemoments.cond);
547
548
549
                fprintf('\n')
                tittxt = 'Draw with SMALLEST condition number';
                fprintf('\n')
550
                disp(['Testing ',tittxt, '. Press ENTER']), pause(5),
551
552
                if ~iload,
                    [idehess_min, idemoments_min, idemodel_min, idelre_min, derivatives_info_min] = ...
553
                        identification_analysis(pdraws(jmin,:),indx,indexo,options_ident,dataset_, prior_exist, name_tex,1);
554
555
                    save([IdentifDirectoryName '/' M_.fname '_identif.mat'], 'idehess_min', 'idemoments_min','idemodel_min', 'idelre_min', 'jmin', '-append');
                end
556
                disp_identification(pdraws(jmin,:), idemodel_min, idemoments_min, name,1);
557
                close all,
558
559
560
                if ~options_.nograph,
                    plot_identification(pdraws(jmin,:),idemoments_min,idehess_min,idemodel_min,idelre_min,1,tittxt,name,IdentifDirectoryName);
                end
561
562
563
564
            else
                for j=1:length(jcrit),
                    tittxt = ['Rank deficient draw n. ',int2str(j)];
                    fprintf('\n')
565
                    disp(['Testing ',tittxt, '. Press ENTER']), pause(5),
566
567
                    if ~iload,
                        [idehess_(j), idemoments_(j), idemodel_(j), idelre_(j), derivatives_info_(j)] = ...
568
                            identification_analysis(pdraws(jcrit(j),:),indx,indexo,options_ident,dataset_, prior_exist, name_tex,1);
569
                    end
570
                    disp_identification(pdraws(jcrit(j),:), idemodel_(j), idemoments_(j), name,1);
571
                    close all,
572
573
574
                    if ~options_.nograph,
                        plot_identification(pdraws(jcrit(j),:),idemoments_(j),idehess_(j),idemodel_(j),idelre_(j),1,tittxt,name,IdentifDirectoryName);
                    end
575
576
577
578
                end
                if ~iload,
                    save([IdentifDirectoryName '/' M_.fname '_identif.mat'], 'idehess_', 'idemoments_','idemodel_', 'idelre_', 'jcrit', '-append');
                end
579
            end
580
        end
581
    end
582
end
Marco Ratto's avatar
Marco Ratto committed
583
584

if exist('OCTAVE_VERSION')
585
    warning('on'),
Marco Ratto's avatar
Marco Ratto committed
586
else
587
    warning on,
Marco Ratto's avatar
Marco Ratto committed
588
end
589
590
591
592
593

disp(' ')
disp(['==== Identification analysis completed ====' ]),
disp(' ')
disp(' ')