identification_analysis.m 17.7 KB
Newer Older
1
function [ide_hess, ide_moments, ide_model, ide_lre, derivatives_info, info, options_ident] = identification_analysis(params,indx,indexo,options_ident,dataset_,dataset_info, prior_exist,name_tex,init,tittxt,bounds)
2
% function [ide_hess, ide_moments, ide_model, ide_lre, derivatives_info, info] = identification_analysis(params,indx,indexo,options_ident,data_info, prior_exist,name_tex,init,analyis_type)
3
4
5
6
7
8
9
% given the parameter vector params, wraps all identification analyses
%
% INPUTS
%    o params             [array] parameter values for identification checks
%    o indx               [array] index of estimated parameters
%    o indexo             [array] index of estimated shocks
%    o options_ident      [structure] identification options
10
11
%    o dataset_           [structure] the dataset after required transformation
%    o dataset_info       [structure] Various informations about the dataset (descriptive statistics and missing observations) info for Kalman Filter
12
%    o prior_exist        [integer]
13
14
%                           =1 when prior exists and indentification is checked only for estimated params and shocks
%                           =0 when prior is not defined and indentification is checked for all params and shocks
15
%    o name_tex           [char] list of tex names
16
%    o init               [integer] flag  for initialization of persistent vars
17
18
%    o tittxt             [string]  string indicating the title text for
%                                   graphs and figures
19
%
20
21
22
23
24
25
26
% OUTPUTS
%    o ide_hess           [structure] identification results using Asymptotic Hessian
%    o ide_moments        [structure] identification results using theoretical moments
%    o ide_model          [structure] identification results using reduced form solution
%    o ide_lre            [structure] identification results using LRE model
%    o derivatives_info   [structure] info about analytic derivs
%    o info               output from dynare resolve
27
%
28
29
30
% SPECIAL REQUIREMENTS
%    None

Stéphane Adjemian's avatar
Stéphane Adjemian committed
31
% Copyright (C) 2008-2018 Dynare Team
32
33
34
35
36
37
38
39
40
41
42
43
44
%
% 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.
%
45
% You should have received a copy of the GNU General Public License
46
47
48
49
50
51
52
53
% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.

global oo_ M_ options_ bayestopt_ estim_params_
persistent indH indJJ indLRE

nparam=length(params);
np=length(indx);
offset=nparam-np;
54
if ~isempty(estim_params_)
55
56
    M_ = set_all_parameters(params,estim_params_,M_);
end
57
58
59
60
61
62
63

nlags = options_ident.ar;
useautocorr = options_ident.useautocorr;
advanced = options_ident.advanced;
replic = options_ident.replic;
periods = options_ident.periods;
max_dim_cova_group = options_ident.max_dim_cova_group;
64
normalize_jacobians = options_ident.normalize_jacobians;
65
kron_flag = options_ident.analytic_derivation_mode;
66

67
68
69
70
71
72
73
74
[I,J]=find(M_.lead_lag_incidence');

ide_hess = struct();
ide_moments = struct();
ide_model = struct();
ide_lre = struct();
derivatives_info = struct();

75
[A,B,ys,info,M_,options_,oo_] = dynare_resolve(M_,options_,oo_);
76
if info(1)==0
77
78
79
    oo0=oo_;
    tau=[oo_.dr.ys(oo_.dr.order_var); vec(A); dyn_vech(B*M_.Sigma_e*B')];
    yy0=oo_.dr.ys(I);
80
    [residual, g1 ] = feval([M_.fname,'.dynamic'],yy0, ...
81
82
                            repmat(oo_.exo_steady_state',[M_.maximum_exo_lag+M_.maximum_exo_lead+1]), M_.params, ...
                            oo_.dr.ys, 1);
83
84
    vg1 = [oo_.dr.ys(oo_.dr.order_var); vec(g1)];

85
    [JJ, H, gam, gp, dA, dOm, dYss] = getJJ(A, B, estim_params_, M_,oo0,options_,kron_flag,indx,indexo,bayestopt_.mf2,nlags,useautocorr);
86
87
88
    derivatives_info.DT=dA;
    derivatives_info.DOm=dOm;
    derivatives_info.DYss=dYss;
89
    if init
90
        indJJ = (find(max(abs(JJ'),[],1)>1.e-8));
91
        if isempty(indJJ) && any(any(isnan(JJ)))
92
93
94
            error('There are NaN in the JJ matrix. Please check whether your model has units roots and you forgot to set diffuse_filter=1.' )
        elseif any(any(isnan(gam)))
            error('There are NaN''s in the theoretical moments: make sure that for non-stationary models stationary transformations of non-stationary observables are used for checking identification. [TIP: use first differences].')
95
        end
96
        while length(indJJ)<nparam && nlags<10
97
            disp('The number of moments with non-zero derivative is smaller than the number of parameters')
98
            disp(['Try increasing ar = ', int2str(nlags+1)])
99
            nlags=nlags+1;
100
            [JJ, H, gam, gp, dA, dOm, dYss] = getJJ(A, B, estim_params_, M_,oo0,options_,kron_flag,indx,indexo,bayestopt_.mf2,nlags,useautocorr);
101
102
103
            derivatives_info.DT=dA;
            derivatives_info.DOm=dOm;
            derivatives_info.DYss=dYss;
104
            options_.ar=nlags;
105
            options_ident.ar=nlags;
106
            indJJ = (find(max(abs(JJ'),[],1)>1.e-8));
107
        end
108
        if length(indJJ)<nparam
109
            disp('The number of moments with non-zero derivative is smaller than the number of parameters')
110
111
            disp('up to 10 lags: check your model')
            disp('Either further increase ar or reduce the list of estimated parameters')
112
            error('identification_analysis: there are not enough moments and too many parameters'),
113
        end
114
115
        indH = (find(max(abs(H'),[],1)>1.e-8));
        indLRE = (find(max(abs(gp'),[],1)>1.e-8));
116
117
118
119
120
    end
    TAU(:,1)=tau(indH);
    LRE(:,1)=vg1(indLRE);
    GAM(:,1)=gam(indJJ);
    siJ = (JJ(indJJ,:));
121
    siH = (H(indH,:));
122
123
124
    siLRE = (gp(indLRE,:));
    ide_strength_J=NaN(1,nparam);
    ide_strength_J_prior=NaN(1,nparam);
125
    if init
126
        ide_uncert_unnormaliz = NaN(1,nparam);
127
        if prior_exist
128
            offset_prior=0;
129
            if ~isempty(estim_params_.var_exo)
130
                normaliz_prior_std = bayestopt_.p2(1:estim_params_.nvx)'; % normalize with prior standard deviation
131
                offset_prior=offset_prior+estim_params_.nvx+estim_params_.nvn;
132
            else
133
                normaliz_prior_std=[];
134
            end
135
            if ~isempty(estim_params_.corrx)
136
137
                normaliz_prior_std = [normaliz_prior_std bayestopt_.p2(offset_prior+1:offset_prior+estim_params_.ncx)']; % normalize with prior standard deviation
                offset_prior=offset_prior+estim_params_.ncx+estim_params_.ncn;
138
            end
139
            if ~isempty(estim_params_.param_vals)
140
                normaliz_prior_std = [normaliz_prior_std bayestopt_.p2(offset_prior+1:offset_prior+estim_params_.np)']; % normalize with prior standard deviation
141
142
            end
            %                         normaliz = max([normaliz; normaliz1]);
143
%             normaliz1(isinf(normaliz1)) = 1;
144

145
        else
146
            normaliz_prior_std = NaN(1,nparam);
147
        end
148
        try
149
150
151
            options_.irf = 0;
            options_.noprint = 1;
            options_.order = 1;
152
            options_.SpectralDensity.trigger = 0;
153
            options_.periods = periods+100;
154
            if options_.kalman_algo > 2
155
156
                options_.kalman_algo = 1;
            end
157
            analytic_derivation = options_.analytic_derivation;
158
            options_.analytic_derivation = -2;
159
            info = stoch_simul(options_.varobs);
160
            dataset_ = dseries(oo_.endo_simul(options_.varobs_id,100+1:end)',dates('1Q1'), options_.varobs);
161
            derivatives_info.no_DLIK=1;
162
            bounds = prior_bounds(bayestopt_, options_.prior_trunc); %reset bounds as lb and ub must only be operational during mode-finding
163
164
            [fval,info,cost_flag,DLIK,AHess,ys,trend_coeff,M_,options_,bayestopt_,oo_] = dsge_likelihood(params',dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_,derivatives_info);
            %                 fval = DsgeLikelihood(xparam1,data_info,options_,M_,estim_params_,bayestopt_,oo_);
165
            options_.analytic_derivation = analytic_derivation;
166
            AHess=-AHess;
167
            if min(eig(AHess))<-1.e-10
168
                error('identification_analysis: Analytic Hessian is not positive semi-definite!')
169
            end
170
            %             chol(AHess);
171
172
173
174
            ide_hess.AHess= AHess;
            deltaM = sqrt(diag(AHess));
            iflag=any((deltaM.*deltaM)==0);
            tildaM = AHess./((deltaM)*(deltaM'));
175
            if iflag || rank(AHess)>rank(tildaM)
176
177
178
179
180
181
                [ide_hess.cond, ide_hess.ind0, ide_hess.indno, ide_hess.ino, ide_hess.Mco, ide_hess.Pco] = identification_checks(AHess, 1);
            else
                [ide_hess.cond, ide_hess.ind0, ide_hess.indno, ide_hess.ino, ide_hess.Mco, ide_hess.Pco] = identification_checks(tildaM, 1);
            end
            indok = find(max(ide_hess.indno,[],1)==0);
            cparam(indok,indok) = inv(AHess(indok,indok));
182
            ide_uncert_unnormaliz(indok) = sqrt(diag(cparam(indok,indok)))';
183
184
185
            cmm = NaN(size(siJ,1),size(siJ,1));
            ind1=find(ide_hess.ind0);
            cmm = siJ(:,ind1)*((AHess(ind1,ind1))\siJ(:,ind1)');
186
187
            temp1=((AHess(ind1,ind1))\siH(:,ind1)');
            diag_chh=sum(siH(:,ind1)'.*temp1)';
188
            %             chh = siH(:,ind1)*((AHess(ind1,ind1))\siH(:,ind1)');
189
            ind1=ind1(ind1>offset);
190
            clre = siLRE(:,ind1-offset)*((AHess(ind1,ind1))\siLRE(:,ind1-offset)');
191
            rhoM=sqrt(1./diag(inv(tildaM(indok,indok))));
192
            %             deltaM = deltaM.*abs(params');
193
            flag_score=1;
194
        catch
195
            %             replic = max([replic, nparam*(nparam+1)/2*10]);
196
            replic = max([replic, length(indJJ)*3]);
197
            cmm = simulated_moment_uncertainty(indJJ, periods, replic,options_,M_,oo_);
198
            %             [V,D,W]=eig(cmm);
199
200
            sd=sqrt(diag(cmm));
            cc=cmm./(sd*sd');
201
            if isoctave || matlab_ver_less_than('8.3')
202
203
                [V,D]=eig(cc);
                %fix for older Matlab versions that do not support computing left eigenvalues, see http://mathworks.com/help/releases/R2012b/matlab/ref/eig.html
204
                [W,~] = eig(cc.');
205
206
207
208
                W = conj(W);
            else
                [V,D,W]=eig(cc);
            end
209
210
211
212
            id=find(diag(D)>1.e-8);
            siTMP=siJ./repmat(sd,[1 nparam]);
            MIM=(siTMP'*V(:,id))*(D(id,id)\(W(:,id)'*siTMP));
            clear siTMP;
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
            %           MIM=siJ(:,indok)'*(cmm\siJ(:,indok));
            %           look for independent moments!
            % % %             sd=sqrt(diag(cmm));
            % % %             cc=cmm./(sd*sd');
            % % %             ix=[];
            % % %             for jc=1:length(cmm),
            % % %                 jcheck=find(abs(cc(:,jc))>(1-1.e-6));
            % % %                 ix=[ix; jcheck(jcheck>jc)];
            % % %             end
            % % %             iy=find(~ismember([1:length(cmm)],ix));
            % % %             indJJ=indJJ(iy);
            % % %             GAM=GAM(iy);
            % % %             cmm=cmm(iy,iy);
            % % %             siJ = (JJ(indJJ,:));
            % % %             MIM=siJ'*(cmm\siJ);
228
229
230
231
            ide_hess.AHess= MIM;
            deltaM = sqrt(diag(MIM));
            iflag=any((deltaM.*deltaM)==0);
            tildaM = MIM./((deltaM)*(deltaM'));
232
            if iflag || rank(MIM)>rank(tildaM)
233
234
235
236
237
                [ide_hess.cond, ide_hess.ind0, ide_hess.indno, ide_hess.ino, ide_hess.Mco, ide_hess.Pco] = identification_checks(MIM, 1);
            else
                [ide_hess.cond, ide_hess.ind0, ide_hess.indno, ide_hess.ino, ide_hess.Mco, ide_hess.Pco] = identification_checks(tildaM, 1);
            end
            indok = find(max(ide_hess.indno,[],1)==0);
238
239
            %             rhoM=sqrt(1-1./diag(inv(tildaM)));
            %             rhoM=(1-1./diag(inv(tildaM)));
240
            ind1=find(ide_hess.ind0);
241
242
            temp1=((MIM(ind1,ind1))\siH(:,ind1)');
            diag_chh=sum(siH(:,ind1)'.*temp1)';
243
            %             chh = siH(:,ind1)*((MIM(ind1,ind1))\siH(:,ind1)');
244
            ind1=ind1(ind1>offset);
245
            clre = siLRE(:,ind1-offset)*((MIM(ind1,ind1))\siLRE(:,ind1-offset)');
246
            if ~isempty(indok)
247
                rhoM(indok)=sqrt(1./diag(inv(tildaM(indok,indok))));
248
                ide_uncert_unnormaliz(indok) = (sqrt(diag(inv(tildaM(indok,indok))))./deltaM(indok))'; %sqrt(diag(inv(MIM(indok,indok))))';
249
250
            end
            %             deltaM = deltaM.*abs(params')
251
            flag_score=0;
252
253
254
255
256
257
        end        
        ide_strength_J(indok) = (1./(ide_uncert_unnormaliz(indok)'./abs(params(indok)')));
        ide_strength_J_prior(indok) = (1./(ide_uncert_unnormaliz(indok)'./normaliz_prior_std(indok)'));
        %ide_strength_J(params==0)=1./ide_uncert_unnormaliz(params==0)';
        sensitivity_zero_pos=find(isinf(deltaM));
        deltaM_prior = deltaM.*abs(normaliz_prior_std');
258
        deltaM = deltaM.*abs(params');
259
        %deltaM(params==0)=deltaM_prior(params==0);
260
        quant = siJ./repmat(sqrt(diag(cmm)),1,nparam);
261
        if size(quant,1)==1
262
            siJnorm = abs(quant).*normaliz_prior_std;
263
        else
264
            siJnorm = vnorm(quant).*normaliz_prior_std;
265
        end
266
267
        %                 siJnorm = vnorm(siJ(inok,:)).*normaliz;
        quant=[];
268
269
270
271
272
        %         inok = find((abs(TAU)<1.e-8));
        %         isok = find((abs(TAU)>=1.e-8));
        %         quant(isok,:) = siH(isok,:)./repmat(TAU(isok,1),1,nparam);
        %         quant(inok,:) = siH(inok,:)./repmat(mean(abs(TAU)),length(inok),nparam);
        %         quant = siH./repmat(sqrt(diag(chh)),1,nparam);
273
274
275
        iy = find(diag_chh);
        indH=indH(iy);
        siH=siH(iy,:);
276
        if ~isempty(iy)
277
            quant = siH./repmat(sqrt(diag_chh(iy)),1,nparam);
278
            if size(quant,1)==1
279
                siHnorm = abs(quant).*normaliz_prior_std;
280
            else
281
                siHnorm = vnorm(quant).*normaliz_prior_std;
282
            end
283
284
285
        else
            siHnorm = [];
        end
286
287
        %                 siHnorm = vnorm(siH./repmat(TAU,1,nparam)).*normaliz;
        quant=[];
288
289
290
291
        %         inok = find((abs(LRE)<1.e-8));
        %         isok = find((abs(LRE)>=1.e-8));
        %         quant(isok,:) = siLRE(isok,:)./repmat(LRE(isok,1),1,np);
        %         quant(inok,:) = siLRE(inok,:)./repmat(mean(abs(LRE)),length(inok),np);
292
293
294
295
        diag_clre = diag(clre);
        iy = find(diag_clre);
        indLRE=indLRE(iy);
        siLRE=siLRE(iy,:);
296
        if ~isempty(iy)
297
            quant = siLRE./repmat(sqrt(diag_clre(iy)),1,np);
298
            if size(quant,1)==1
299
                siLREnorm = abs(quant).*normaliz_prior_std(offset+1:end);
300
            else
301
                siLREnorm = vnorm(quant).*normaliz_prior_std(offset+1:end);
302
            end
303
304
305
        else
            siLREnorm=[];
        end
306
        %                 siLREnorm = vnorm(siLRE./repmat(LRE,1,nparam-offset)).*normaliz(offset+1:end);
307
308
309
310
        ide_hess.ide_strength_J=ide_strength_J;
        ide_hess.ide_strength_J_prior=ide_strength_J_prior;
        ide_hess.deltaM=deltaM;
        ide_hess.deltaM_prior=deltaM_prior;
311
312
        ide_hess.sensitivity_zero_pos=sensitivity_zero_pos;
        ide_hess.identified_parameter_indices=indok;
313
314
315
316
        ide_moments.siJnorm=siJnorm;
        ide_model.siHnorm=siHnorm;
        ide_lre.siLREnorm=siLREnorm;
        ide_hess.flag_score=flag_score;
317
318
    end
    if normalize_jacobians
319
320
321
322
323
324
325
326
327
328
329
        normH = max(abs(siH)')';
        normH = normH(:,ones(nparam,1));
        normJ = max(abs(siJ)')';
        normJ = normJ(:,ones(nparam,1));
        normLRE = max(abs(siLRE)')';
        normLRE = normLRE(:,ones(size(gp,2),1));
    else
        normH = 1;
        normJ = 1;
        normLRE = 1;
    end
330
331
332
333
334
335
336
337
338
    ide_moments.indJJ=indJJ;
    ide_model.indH=indH;
    ide_lre.indLRE=indLRE;
    ide_moments.siJ=siJ;
    ide_model.siH=siH;
    ide_lre.siLRE=siLRE;
    ide_moments.GAM=GAM;
    ide_model.TAU=TAU;
    ide_lre.LRE=LRE;
339
340
341
342
343
344
345
346
    %     [ide_checks.idemodel_Mco, ide_checks.idemoments_Mco, ide_checks.idelre_Mco, ...
    %         ide_checks.idemodel_Pco, ide_checks.idemoments_Pco, ide_checks.idelre_Pco, ...
    %         ide_checks.idemodel_cond, ide_checks.idemoments_cond, ide_checks.idelre_cond, ...
    %         ide_checks.idemodel_ee, ide_checks.idemoments_ee, ide_checks.idelre_ee, ...
    %         ide_checks.idemodel_ind, ide_checks.idemoments_ind, ...
    %         ide_checks.idemodel_indno, ide_checks.idemoments_indno, ...
    %         ide_checks.idemodel_ino, ide_checks.idemoments_ino] = ...
    %         identification_checks(H(indH,:)./normH(:,ones(nparam,1)),JJ(indJJ,:)./normJ(:,ones(nparam,1)), gp(indLRE,:)./normLRE(:,ones(size(gp,2),1)));
347
    [ide_moments.cond, ide_moments.ind0, ide_moments.indno, ide_moments.ino, ide_moments.Mco, ide_moments.Pco, ide_moments.jweak, ide_moments.jweak_pair] = ...
348
        identification_checks(JJ(indJJ,:)./normJ, 0);
349
    [ide_model.cond, ide_model.ind0, ide_model.indno, ide_model.ino, ide_model.Mco, ide_model.Pco, ide_model.jweak, ide_model.jweak_pair] = ...
350
        identification_checks(H(indH,:)./normH, 0);
351
    [ide_lre.cond, ide_lre.ind0, ide_lre.indno, ide_lre.ino, ide_lre.Mco, ide_lre.Pco, ide_lre.jweak, ide_lre.jweak_pair] = ...
352
353
354
        identification_checks(gp(indLRE,:)./normLRE, 0);
    normJ=1;
    [U, S, V]=svd(JJ(indJJ,:)./normJ,0);
355
    S=diag(S);
356
    S=[S;zeros(size(JJ,2)-length(indJJ),1)];
357
358
359
360
361
362
363
364
    if nparam>8
        ide_moments.S = S([1:4, end-3:end]);
        ide_moments.V = V(:,[1:4, end-3:end]);
    else
        ide_moments.S = S;
        ide_moments.V = V;
    end
    indok = find(max(ide_moments.indno,[],1)==0);
365
    if advanced
366
        [ide_moments.pars, ide_moments.cosnJ] = ident_bruteforce(JJ(indJJ,:)./normJ,max_dim_cova_group,options_.TeX,name_tex,tittxt);
367
    end
368
end