New wrapper for identification analysis, given the parameter set provides results.

Allows more flexible coding and checks for particular parameter sets.

matlab/identification_analysis.m

+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)
+% 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)
+% 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
+%    o data_info          [structure] data info for Kalmna Filter
+%    o prior_exist        [integer] 
+%                           =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
+%    o nem_tex            [char] list of tex names
+%    o init               [integer] flag  for initialization of persistent vars
+%    
+% 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
+%    
+% SPECIAL REQUIREMENTS
+%    None
+
+% Copyright (C) 2008-2011 Dynare Team
+%
+% 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 oo_ M_ options_ bayestopt_ estim_params_
+persistent indH indJJ indLRE
+
+nparam=length(params);
+np=length(indx);
+offset=nparam-np;
+set_all_parameters(params);
+
+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;
+[I,J]=find(M_.lead_lag_incidence');
+
+ide_hess = struct();
+ide_moments = struct();
+ide_model = struct();
+ide_lre = struct();
+derivatives_info = struct();
+
+[A,B,ys,info]=dynare_resolve;
+if info(1)==0,
+    oo0=oo_;
+    tau=[oo_.dr.ys(oo_.dr.order_var); vec(A); dyn_vech(B*M_.Sigma_e*B')];
+    yy0=oo_.dr.ys(I);
+    [residual, g1 ] = feval([M_.fname,'_dynamic'],yy0, ...
+        oo_.exo_steady_state', M_.params, ...
+        oo_.dr.ys, 1);
+    vg1 = [oo_.dr.ys(oo_.dr.order_var); vec(g1)];
+
+    [JJ, H, gam, gp, dA, dOm, dYss] = getJJ(A, B, M_,oo0,options_,0,indx,indexo,bayestopt_.mf2,nlags,useautocorr);
+    derivatives_info.DT=dA;
+    derivatives_info.DOm=dOm;
+    derivatives_info.DYss=dYss;
+    if init,
+        indJJ = (find(max(abs(JJ'))>1.e-8));
+        indH = (find(max(abs(H'))>1.e-8));
+        indLRE = (find(max(abs(gp'))>1.e-8));
+    end
+    ide_moments.indJJ=indJJ;
+    ide_model.indH=indH;
+    ide_lre.indLRE=indLRE;
+    TAU(:,1)=tau(indH);
+    LRE(:,1)=vg1(indLRE);
+    GAM(:,1)=gam(indJJ);
+    siJ = (JJ(indJJ,:));
+    siH = (H(indH,:));   
+    siLRE = (gp(indLRE,:));
+    normH = max(abs(siH)')';
+    normJ = max(abs(siJ)')';
+    normLRE = max(abs(siLRE)')';
+    ide_moments.siJ=siJ;
+    ide_model.siH=siH;
+    ide_lre.siLRE=siLRE;
+    ide_moments.GAM=GAM;
+    ide_model.TAU=TAU;
+    ide_lre.LRE=LRE;
+%     [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)));
+    [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] = ...
+        identification_checks(JJ(indJJ,:)./normJ(:,ones(nparam,1)), 0);
+    [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] = ...
+        identification_checks(H(indH,:)./normH(:,ones(nparam,1)), 0);
+    [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] = ...
+        identification_checks(gp(indLRE,:)./normLRE(:,ones(size(gp,2),1)), 0);
+    
+    indok = find(max(ide_moments.indno,[],1)==0);
+    ide_strength_J=NaN(1,nparam);
+    ide_strength_J_prior=NaN(1,nparam);
+    if advanced,
+        [ide_moments.pars, ide_moments.cosnJ] = ident_bruteforce(JJ(indJJ,:)./normJ(:,ones(nparam,1)),max_dim_cova_group,options_.TeX,name_tex);
+    end
+    if init, %~isempty(indok),
+        normaliz = abs(params);
+        if prior_exist,
+            if ~isempty(estim_params_.var_exo),
+                normaliz1 = estim_params_.var_exo(:,7); % normalize with prior standard deviation
+            else
+                normaliz1=[];
+            end
+            if ~isempty(estim_params_.param_vals),
+                normaliz1 = [normaliz1; estim_params_.param_vals(:,7)]'; % normalize with prior standard deviation
+            end
+            %                         normaliz = max([normaliz; normaliz1]);
+            normaliz1(isinf(normaliz1)) = 1;
+            
+        else
+            normaliz1 = NaN(1,nparam);
+        end
+        try,
+            options_.irf = 0;
+            options_.noprint = 1;
+            options_.order = 1;
+            options_.periods = data_info.gend+100;
+            options_.kalman_algo = 1;
+            info = stoch_simul(options_.varobs);
+            datax=oo_.endo_simul(options_.varobs_id,100+1:end);
+            %                         datax=data;
+            derivatives_info.no_DLIK=1;
+            [fval,cost_flag,ys,trend_coeff,info,DLIK,AHess] = DsgeLikelihood(params',data_info.gend,datax,data_info.data_index,data_info.number_of_observations,data_info.no_more_missing_observations,derivatives_info);
+            AHess=-AHess;
+            ide_hess.AHess= AHess;
+            deltaM = sqrt(diag(AHess));
+            iflag=any((deltaM.*deltaM)==0);
+            tildaM = AHess./((deltaM)*(deltaM'));
+            if iflag || rank(AHess)>rank(tildaM),
+                [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));
+            normaliz(indok) = sqrt(diag(cparam(indok,indok)))';
+            cmm = siJ*((AHess)\siJ');
+            rhoM=sqrt(1./diag(inv(tildaM(indok,indok))));
+            deltaM = deltaM.*params';
+            flag_score=1;
+        catch,
+            replic = max([replic, length(indJJ)*3]);
+            cmm = simulated_moment_uncertainty(indJJ, periods, replic);
+%             MIM=siJ(:,indok)'*(cmm\siJ(:,indok));
+            MIM=siJ'*(cmm\siJ);
+            ide_hess.AHess= MIM;
+            deltaM = sqrt(diag(MIM));
+            iflag=any((deltaM.*deltaM)==0);
+            tildaM = MIM./((deltaM)*(deltaM'));
+            if iflag || rank(MIM)>rank(tildaM),
+                [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);
+%             rhoM=sqrt(1-1./diag(inv(tildaM)));
+%             rhoM=(1-1./diag(inv(tildaM)));
+            rhoM(indok)=sqrt(1./diag(inv(tildaM(indok,indok))));
+            normaliz(indok) = (sqrt(diag(inv(tildaM(indok,indok))))./deltaM(indok))'; %sqrt(diag(inv(MIM(indok,indok))))';
+            deltaM = deltaM.*params';
+            flag_score=0;
+        end
+        ide_strength_J(indok) = (1./(normaliz(indok)'./abs(params(indok)')));
+        ide_strength_J_prior(indok) = (1./(normaliz(indok)'./normaliz1(indok)'));
+        ide_strength_J(params==0)=ide_strength_J_prior(params==0);
+        quant = siJ./repmat(sqrt(diag(cmm)),1,nparam);
+        siJnorm = vnorm(quant).*normaliz1;
+        %                 siJnorm = vnorm(siJ(inok,:)).*normaliz;
+        quant=[];
+        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);
+        siHnorm = vnorm(quant).*normaliz1;
+        %                 siHnorm = vnorm(siH./repmat(TAU,1,nparam)).*normaliz;
+        quant=[];
+        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);
+        siLREnorm = vnorm(quant).*normaliz1(offset+1:end);
+        %                 siLREnorm = vnorm(siLRE./repmat(LRE,1,nparam-offset)).*normaliz(offset+1:end);
+        ide_hess.ide_strength_J=ide_strength_J; 
+        ide_hess.ide_strength_J_prior=ide_strength_J_prior; 
+        ide_moments.siJnorm=siJnorm; 
+        ide_model.siHnorm=siHnorm; 
+        ide_lre.siLREnorm=siLREnorm; 
+        ide_hess.flag_score=flag_score; 
+    end,
+end