dynare_sensitivity.m 17 KB
Newer Older
1
2
3
4
5
6
function x0=dynare_sensitivity(options_gsa)
% Frontend to the Sensitivity Analysis Toolbox for DYNARE
%
% Reference:
% M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006.

7
% Copyright (C) 2008-2011 Dynare Team
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
%
% 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_

fname_ = M_.fname;
Marco Ratto's avatar
Marco Ratto committed
27
28
29
if ~isfield(M_,'dname'),
    M_.dname = M_.fname;
end
30
31
32
33
34
35
36
lgy_ = M_.endo_names;
x0=[];

options_gsa = set_default_option(options_gsa,'datafile',[]);
if isfield(options_gsa,'nograph'),
    options_.nograph=options_gsa.nograph;
end
37
38
39
40
41
42
if isfield(options_gsa,'nodisplay'),
    options_.nodisplay=options_gsa.nodisplay;
end
if isfield(options_gsa,'graph_format'),
    options_.graph_format=options_gsa.graph_format;
end
43
44
45
46
if isfield(options_gsa,'mode_file'),
    options_.mode_file=options_gsa.mode_file;
end

ratto's avatar
ratto committed
47
options_.order = 1;
48

49
if ~isempty(options_gsa.datafile) || isempty(bayestopt_),
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
    options_.datafile = options_gsa.datafile;
    if isfield(options_gsa,'first_obs'),
        options_.first_obs=options_gsa.first_obs;
    end
    if isfield(options_gsa,'nobs'),
        options_.nobs=options_gsa.nobs;
    end
    if isfield(options_gsa,'presample'),
        options_.presample=options_gsa.presample;
    end
    if isfield(options_gsa,'prefilter'),
        options_.prefilter=options_gsa.prefilter;
    end
    if isfield(options_gsa,'loglinear'),
        options_.loglinear=options_gsa.loglinear;
    end
66
67
68
    if isfield(options_gsa,'lik_init'),
        options_.lik_init=options_gsa.lik_init;
    end
69
70
71
    if isfield(options_gsa,'kalman_algo'),
        options_.kalman_algo=options_gsa.kalman_algo;
    end
72
    options_.mode_compute = 0;
73
74
    options_.filtered_vars = 1;
    options_.plot_priors = 0;
Marco Ratto's avatar
Marco Ratto committed
75
76
%     [data,rawdata,xparam1,data_info]=dynare_estimation_init([],fname_,1);
    [dataset_,xparam1, M_, options_, oo_, estim_params_,bayestopt_]=dynare_estimation_init(M_.endo_names,fname_,1, M_, options_, oo_, estim_params_, bayestopt_);
77
    % computes a first linear solution to set up various variables
78
79
80
81
else
    if isempty(options_.qz_criterium)
        options_.qz_criterium = 1+1e-6;
    end    
82
end
83
[make,my,day,punk,M_,options_,oo_] = dynare_resolve(M_,options_,oo_);
84
85
86
87
88
89
90

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);
    options_gsa = set_default_option(options_gsa,'trans_ident',0);
    options_gsa = set_default_option(options_gsa,'load_ident_files',0);
91
92
    options_gsa = set_default_option(options_gsa,'ar',1);
    options_gsa = set_default_option(options_gsa,'useautocorr',0);
93
    options_.ar = options_gsa.ar;
94
95
96
97
98
    if options_gsa.morris==0,
        disp('The option morris = 0 is no longer supported (Type I errors)')
        disp('This option is reset at morris = 2 (local identification analysis).')
        options_gsa.morris=2;
    end
99
    if options_gsa.morris==2,
100
101
102
103
104
105
106
107
108
109
110
111
        if isfield(options_,'options_ident'),
            options_.options_ident.load_ident_files = options_gsa.load_ident_files;
            options_.options_ident.useautocorr = options_gsa.useautocorr;
            options_.options_ident.ar = options_gsa.ar;            
            options_ident=options_.options_ident;
        else
            options_ident=[];
            options_ident = set_default_option(options_ident,'load_ident_files',options_gsa.load_ident_files);
            options_ident = set_default_option(options_ident,'useautocorr',options_gsa.useautocorr);
            options_ident = set_default_option(options_ident,'ar',options_gsa.ar);
            options_.options_ident = options_ident;
        end
112
113
114
115
116
117
118
119
120
    end
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);
121
options_gsa = set_default_option(options_gsa,'neighborhood_width',0);
122
123
124
125
126
127
128
129
130
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);
131
options_gsa = set_default_option(options_gsa,'alpha2_stab',0);
132
options_gsa = set_default_option(options_gsa,'ksstat',0.1);
Marco Ratto's avatar
Marco Ratto committed
133
134
options_gsa = set_default_option(options_gsa,'pvalue_ks',0.001);
options_gsa = set_default_option(options_gsa,'pvalue_corr',0.001);
135
%options_gsa = set_default_option(options_gsa,'load_mh',0);
136
137
138
139
% REDFORM mapping
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.001);
140
options_gsa = set_default_option(options_gsa,'alpha2_redform',0);
141
142
143
144
145
146
147
148
149
150
options_gsa = set_default_option(options_gsa,'namendo',[]);
options_gsa = set_default_option(options_gsa,'namlagendo',[]);
options_gsa = set_default_option(options_gsa,'namexo',[]);
% RMSE mapping
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,'pfilt_rmse',0.1);
options_gsa = set_default_option(options_gsa,'istart_rmse',options_.presample+1);
options_gsa = set_default_option(options_gsa,'alpha_rmse',0.001);
151
options_gsa = set_default_option(options_gsa,'alpha2_rmse',0);
152

153
if options_gsa.redform && options_gsa.neighborhood_width==0 && isempty(options_gsa.threshold_redform),
154
155
156
157
    options_gsa.pprior=1;
    options_gsa.ppost=0;
end

158
159
160
161
162
163
164
if options_gsa.morris>2,
    disp('The option morris = 3 is no longer supported')
    disp('the option is reset at morris = 1 .')
    options_gsa.morris=1;
end
   
if options_gsa.morris==1,
165
166
167
168
169
170
171
172
173
174
175
    if ~options_gsa.identification,
        options_gsa.redform=1;
    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_gsa.ksstat=1;
176
177
    options_gsa.pvalue_ks=0;
    options_gsa.pvalue_corr=0;
178
179
180
181
%     if options_gsa.morris==3,
%         options_gsa = set_default_option(options_gsa,'Nsam',256);
%         OutputDirectoryName = CheckPath('gsa/identif',M_.dname);
%     else
182
        OutputDirectoryName = CheckPath('gsa/screen',M_.dname);
183
%     end
184
else
185
    OutputDirectoryName = CheckPath('gsa',M_.dname);
186
187
end

188
% options_.opt_gsa = options_gsa;
189

190
191
if (options_gsa.load_stab || options_gsa.load_rmse || options_gsa.load_redform) ...
        && options_gsa.pprior,
192
    filetoload=[OutputDirectoryName '/' fname_ '_prior.mat'];
Marco Ratto's avatar
Marco Ratto committed
193
    if ~exist(filetoload),
194
195
196
197
198
199
200
201
202
203
        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'),
204
205
206
207
            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))) || ...
208
209
210
211
212
213
214
215
216
217
218
                    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

219
if options_gsa.stab && ~options_gsa.ppost,
220
    x0 = stab_map_(OutputDirectoryName,options_gsa);
221
222
223
224
225
    if isempty(x0),
        disp(' ')
        disp('Sensitivity computations stopped: no parameter set provided a unique solution')
        return
    end
226
227
228
229
230
231
232
end

% reduced form
% redform_map(namendo, namlagendo, namexo, icomp, pprior, ilog, threshold)

options_.opt_gsa = options_gsa;
if options_gsa.identification,
233
    map_ident_(OutputDirectoryName,options_gsa);
234
235
end

236
237
if options_gsa.redform && ~isempty(options_gsa.namendo) ...
        && ~options_gsa.ppost,
238
239
240
241
242
243
244
245
246
    if strmatch(':',options_gsa.namendo,'exact'),
        options_gsa.namendo=M_.endo_names;
    end
    if strmatch(':',options_gsa.namexo,'exact'),
        options_gsa.namexo=M_.exo_names;
    end
    if strmatch(':',options_gsa.namlagendo,'exact'),
        options_gsa.namlagendo=M_.endo_names;
    end
247
%     options_.opt_gsa = options_gsa;
248
    if options_gsa.morris==1,
249
        redform_screen(OutputDirectoryName,options_gsa);
250
251
    else
        % check existence of the SS_ANOVA toolbox
252
253
        if isempty(options_gsa.threshold_redform) && ...
         ~(exist('gsa_sdp','file')==6 || exist('gsa_sdp','file')==2),
254
255
256
257
258
259
            disp('Download Mapping routines at:')
            disp('http://eemc.jrc.ec.europa.eu/softwareDYNARE-Dowload.htm')
            disp(' ' )
            error('Mapping routines missing!')
        end

260
        redform_map(OutputDirectoryName,options_gsa);
261
262
263
264
265
266
267
268
269
270
271
272
    end
end
% RMSE mapping
% function [rmse_MC, ixx] = filt_mc_(vvarvecm, loadSA, pfilt, alpha, alpha2)
options_.opt_gsa = options_gsa;
if options_gsa.rmse,
    if ~options_gsa.ppost
        if options_gsa.pprior
            a=whos('-file',[OutputDirectoryName,'/',fname_,'_prior'],'logpo2');
        else
            a=whos('-file',[OutputDirectoryName,'/',fname_,'_mc'],'logpo2');
        end
Marco Ratto's avatar
Marco Ratto committed
273
        if exist('OCTAVE_VERSION'),
Marco Ratto's avatar
Marco Ratto committed
274
275
276
277
278
279
280
281
282
283
            aflag=0;
            for ja=1:length(a),
                aflag=aflag+strcmp('logpo2',a(ja).name);
            end
            if aflag==0,
                a=[];
            else
                a=1;
            end
        end
284
        if isempty(a),
285
286
287
288
289
290
           if options_gsa.lik_only,
               options_.smoother=0;
               options_.filter_step_ahead=[];
               options_.forecast=0;
               options_.filtered_vars=0;               
           end
291
%             dynare_MC([],OutputDirectoryName,data,rawdata,data_info);
292
            prior_posterior_statistics('gsa',dataset_);
293
294
295
            if options_.bayesian_irf
                PosteriorIRF('gsa');
            end
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
            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;
314
%     filt_mc_(OutputDirectoryName,data_info);
315
    filt_mc_(OutputDirectoryName,options_gsa);
316
end
317
options_.opt_gsa = options_gsa;
318
319
320
321
322
323


if options_gsa.glue,
    dr_ = oo_.dr;
    if options_gsa.ppost
        load([OutputDirectoryName,'/',fname_,'_post']);
324
        DirectoryName = CheckPath('metropolis',M_.dname);
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
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
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
    else
        if options_gsa.pprior
            load([OutputDirectoryName,'/',fname_,'_prior']);
        else
            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);
    rawdata = rawdata(options_.first_obs:options_.first_obs+gend-1,:);
    if options_.loglinear == 1
        rawdata = log(rawdata);
    end
    if options_.prefilter == 1
        %data = transpose(rawdata-ones(gend,1)*bayestopt_.mean_varobs);
        data = transpose(rawdata-ones(gend,1)*mean(rawdata,1));
    else
        data = transpose(rawdata);
    end
    
    Obs.data = data;
    Obs.time = [1:gend];
    Obs.num  = gend;
    for j=1:size(options_.varobs,1)
        Obs.name{j} = deblank(options_.varobs(j,:));
        vj=deblank(options_.varobs(j,:));
        
        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)];
            Out(j).data = jxj;
            Out(j).time = [pwd,'/',OutputDirectoryName];
        else
            Out(j).data = jxj;
            Out(j).time = [pwd,'/',DirectoryName];
        end
        Out(j).name = vj;
        Out(j).ini  = 'yes';
        Lik(j).name = ['rmse_',vj];
        Lik(j).ini  = 'yes';
        Lik(j).isam = 1;
        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';
            Out1=Out;
        else
            Out1=Out;
        end
        ismoo(j)=jxj;
        
    end
    jsmoo = size(options_.varobs,1);
    for j=1:M_.endo_nbr,
        if ~ismember(j,ismoo),
            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';
                Out1(jsmoo).data = j;
                Out1(jsmoo).time = [pwd,'/',OutputDirectoryName];
            else
                Out1(jsmoo).data = j;
                Out1(jsmoo).time = [pwd,'/',DirectoryName];
            end
            Out1(jsmoo).name = vj;
            Out1(jsmoo).ini  = 'yes';
        end
    end
    tit(M_.exo_names_orig_ord,:) = M_.exo_names;
    for j=1:M_.exo_nbr,
        Exo(j).name = deblank(tit(j,:));    
    end
    if ~options_gsa.ppost
        Lik(size(options_.varobs,1)+1).name = 'logpo';
        Lik(size(options_.varobs,1)+1).ini  = 'yes';
        Lik(size(options_.varobs,1)+1).isam = 1;
        Lik(size(options_.varobs,1)+1).data = -logpo2;
    end
    Sam.name = bayestopt_.name;
    Sam.dim  = [size(x) 0];
    Sam.data = [x];
    
    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;
        Info.TypeofSample='post';
439
        save([OutputDirectoryName,'/',fname_,'_glue_post.mat'], 'Out', 'Sam', 'Lik', 'Obs', 'Rem','Info', 'Exo')
440
441
442
443
444
        %save([fname_,'_post_glue_smooth'], 'Out', 'Sam', 'Lik', 'Obs', 'Rem','Info')
        
    else
        if options_gsa.pprior
            Info.TypeofSample='prior';
445
            save([OutputDirectoryName,'/',fname_,'_glue_prior.mat'], 'Out', 'Sam', 'Lik', 'Obs', 'Rem','Info', 'Exo')
446
447
448
449
450
            %       save([OutputDirectoryName,'/',fname_,'_prior_glue'], 'Out', 'Sam', 'Lik', 'Obs', 'Rem')
            %       Out=Out1;
            %       save([OutputDirectoryName,'/',fname_,'_prior_glue_smooth'], 'Out', 'Sam', 'Lik', 'Obs', 'Rem')
        else
            Info.TypeofSample='mc';
451
            save([OutputDirectoryName,'/',fname_,'_glue_mc.mat'], 'Out', 'Sam', 'Lik', 'Obs', 'Rem','Info', 'Exo')
452
453
454
455
456
457
458
            %       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
    
end