identification_analysis

 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

0001 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)
0002 % 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)
0003 % given the parameter vector params, wraps all identification analyses
0004 %
0005 % INPUTS
0006 %    o params             [array] parameter values for identification checks
0007 %    o indx               [array] index of estimated parameters
0008 %    o indexo             [array] index of estimated shocks
0009 %    o options_ident      [structure] identification options
0010 %    o data_info          [structure] data info for Kalmna Filter
0011 %    o prior_exist        [integer]
0012 %                           =1 when prior exists and indentification is checked only for estimated params and shocks
0013 %                           =0 when prior is not defined and indentification is checked for all params and shocks
0014 %    o nem_tex            [char] list of tex names
0015 %    o init               [integer] flag  for initialization of persistent vars
0016 %
0017 % OUTPUTS
0018 %    o ide_hess           [structure] identification results using Asymptotic Hessian
0019 %    o ide_moments        [structure] identification results using theoretical moments
0020 %    o ide_model          [structure] identification results using reduced form solution
0021 %    o ide_lre            [structure] identification results using LRE model
0022 %    o derivatives_info   [structure] info about analytic derivs
0023 %    o info               output from dynare resolve
0024 %
0025 % SPECIAL REQUIREMENTS
0026 %    None
0027 
0028 % Copyright (C) 2008-2011 Dynare Team
0029 %
0030 % This file is part of Dynare.
0031 %
0032 % Dynare is free software: you can redistribute it and/or modify
0033 % it under the terms of the GNU General Public License as published by
0034 % the Free Software Foundation, either version 3 of the License, or
0035 % (at your option) any later version.
0036 %
0037 % Dynare is distributed in the hope that it will be useful,
0038 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0039 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0040 % GNU General Public License for more details.
0041 %
0042 % You should have received a copy of the GNU General Public License
0043 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0044 
0045 global oo_ M_ options_ bayestopt_ estim_params_
0046 persistent indH indJJ indLRE
0047 
0048 nparam=length(params);
0049 np=length(indx);
0050 offset=nparam-np;
0051 set_all_parameters(params);
0052 
0053 nlags = options_ident.ar;
0054 useautocorr = options_ident.useautocorr;
0055 advanced = options_ident.advanced;
0056 replic = options_ident.replic;
0057 periods = options_ident.periods;
0058 max_dim_cova_group = options_ident.max_dim_cova_group;
0059 normalize_jacobians = options_ident.normalize_jacobians;
0060 [I,J]=find(M_.lead_lag_incidence');
0061 
0062 ide_hess = struct();
0063 ide_moments = struct();
0064 ide_model = struct();
0065 ide_lre = struct();
0066 derivatives_info = struct();
0067 
0068 [A,B,ys,info,M_,options_,oo_] = dynare_resolve(M_,options_,oo_);
0069 if info(1)==0,
0070     oo0=oo_;
0071     tau=[oo_.dr.ys(oo_.dr.order_var); vec(A); dyn_vech(B*M_.Sigma_e*B')];
0072     yy0=oo_.dr.ys(I);
0073     [residual, g1 ] = feval([M_.fname,'_dynamic'],yy0, ...
0074         oo_.exo_steady_state', M_.params, ...
0075         oo_.dr.ys, 1);
0076     vg1 = [oo_.dr.ys(oo_.dr.order_var); vec(g1)];
0077 
0078     [JJ, H, gam, gp, dA, dOm, dYss] = getJJ(A, B, M_,oo0,options_,0,indx,indexo,bayestopt_.mf2,nlags,useautocorr);
0079     derivatives_info.DT=dA;
0080     derivatives_info.DOm=dOm;
0081     derivatives_info.DYss=dYss;
0082     if init,
0083         indJJ = (find(max(abs(JJ'))>1.e-8));
0084         while length(indJJ)<nparam && nlags<10,
0085             disp('The number of moments with non-zero derivative is smaller than the number of parameters')
0086             disp(['Try increasing ar = ', int2str(nlags+1)])           
0087             nlags=nlags+1;
0088             [JJ, H, gam, gp, dA, dOm, dYss] = getJJ(A, B, M_,oo0,options_,0,indx,indexo,bayestopt_.mf2,nlags,useautocorr);
0089             derivatives_info.DT=dA;
0090             derivatives_info.DOm=dOm;
0091             derivatives_info.DYss=dYss;
0092             evalin('caller',['options_ident.ar=',int2str(nlags),';']);
0093         end
0094         if length(indJJ)<nparam,
0095             disp('The number of moments with non-zero derivative is smaller than the number of parameters')
0096             disp('up to 10 lags: check your model')           
0097             disp('Either further increase ar or reduce the list of estimated parameters')           
0098             error('IDETooManyParams',''),
0099         end
0100         indH = (find(max(abs(H'))>1.e-8));
0101         indLRE = (find(max(abs(gp'))>1.e-8));
0102     end
0103     TAU(:,1)=tau(indH);
0104     LRE(:,1)=vg1(indLRE);
0105     GAM(:,1)=gam(indJJ);
0106     siJ = (JJ(indJJ,:));
0107     siH = (H(indH,:));   
0108     siLRE = (gp(indLRE,:));
0109     ide_strength_J=NaN(1,nparam);
0110     ide_strength_J_prior=NaN(1,nparam);
0111     if init, %~isempty(indok),
0112         normaliz = abs(params);
0113         if prior_exist,
0114             if ~isempty(estim_params_.var_exo),
0115                 normaliz1 = estim_params_.var_exo(:,7); % normalize with prior standard deviation
0116             else
0117                 normaliz1=[];
0118             end
0119             if ~isempty(estim_params_.param_vals),
0120                 normaliz1 = [normaliz1; estim_params_.param_vals(:,7)]'; % normalize with prior standard deviation
0121             end
0122             %                         normaliz = max([normaliz; normaliz1]);
0123             normaliz1(isinf(normaliz1)) = 1;
0124             
0125         else
0126             normaliz1 = NaN(1,nparam);
0127         end
0128         try,
0129             options_.irf = 0;
0130             options_.noprint = 1;
0131             options_.order = 1;
0132             options_.periods = data_info.info.ntobs+100;
0133             options_.kalman_algo = 1;
0134             options_.analytic_derivation = -2;
0135             info = stoch_simul(options_.varobs);
0136             data_info.data=oo_.endo_simul(options_.varobs_id,100+1:end);
0137             %                         datax=data;
0138             derivatives_info.no_DLIK=1;
0139             [fval,DLIK,AHess,cost_flag,ys,trend_coeff,info,M_,options_,bayestopt_,oo_] = dsge_likelihood(params',data_info,options_,M_,estim_params_,bayestopt_,oo_,derivatives_info);
0140 %                 fval = DsgeLikelihood(xparam1,data_info,options_,M_,estim_params_,bayestopt_,oo_);
0141             AHess=-AHess;
0142             if min(eig(AHess))<0,
0143                 error('Analytic Hessian is not positive semi-definite!')
0144             end
0145 %             chol(AHess);
0146             ide_hess.AHess= AHess;
0147             deltaM = sqrt(diag(AHess));
0148             iflag=any((deltaM.*deltaM)==0);
0149             tildaM = AHess./((deltaM)*(deltaM'));
0150             if iflag || rank(AHess)>rank(tildaM),
0151                 [ide_hess.cond, ide_hess.ind0, ide_hess.indno, ide_hess.ino, ide_hess.Mco, ide_hess.Pco] = identification_checks(AHess, 1);
0152             else
0153                 [ide_hess.cond, ide_hess.ind0, ide_hess.indno, ide_hess.ino, ide_hess.Mco, ide_hess.Pco] = identification_checks(tildaM, 1);
0154             end
0155             indok = find(max(ide_hess.indno,[],1)==0);
0156             cparam(indok,indok) = inv(AHess(indok,indok));
0157             normaliz(indok) = sqrt(diag(cparam(indok,indok)))';
0158             cmm = NaN(size(siJ,1),size(siJ,1));
0159             ind1=find(ide_hess.ind0);
0160             cmm = siJ(:,ind1)*((AHess(ind1,ind1))\siJ(:,ind1)');
0161             temp1=((AHess(ind1,ind1))\siH(:,ind1)');
0162             diag_chh=sum(siH(:,ind1)'.*temp1)';
0163 %             chh = siH(:,ind1)*((AHess(ind1,ind1))\siH(:,ind1)');
0164             ind1=ind1(ind1>offset);
0165             clre = siLRE(:,ind1-offset)*((AHess(ind1,ind1))\siLRE(:,ind1-offset)');
0166             rhoM=sqrt(1./diag(inv(tildaM(indok,indok))));
0167 %             deltaM = deltaM.*abs(params');
0168             flag_score=1;
0169         catch,
0170             replic = max([replic, length(indJJ)*3]);
0171             cmm = simulated_moment_uncertainty(indJJ, periods, replic);
0172 %             MIM=siJ(:,indok)'*(cmm\siJ(:,indok));
0173 %           look for independent moments!
0174             sd=sqrt(diag(cmm));
0175             cc=cmm./(sd*sd');
0176             ix=[];
0177             for jc=1:length(cmm),
0178                 jcheck=find(abs(cc(:,jc))>(1-1.e-6));
0179                 ix=[ix; jcheck(jcheck>jc)];
0180             end
0181             iy=find(~ismember([1:length(cmm)],ix));
0182             indJJ=indJJ(iy);
0183             GAM=GAM(iy);
0184             cmm=cmm(iy,iy);
0185             siJ = (JJ(indJJ,:));
0186             MIM=siJ'*(cmm\siJ);
0187             ide_hess.AHess= MIM;
0188             deltaM = sqrt(diag(MIM));
0189             iflag=any((deltaM.*deltaM)==0);
0190             tildaM = MIM./((deltaM)*(deltaM'));
0191             if iflag || rank(MIM)>rank(tildaM),
0192                 [ide_hess.cond, ide_hess.ind0, ide_hess.indno, ide_hess.ino, ide_hess.Mco, ide_hess.Pco] = identification_checks(MIM, 1);
0193             else
0194                 [ide_hess.cond, ide_hess.ind0, ide_hess.indno, ide_hess.ino, ide_hess.Mco, ide_hess.Pco] = identification_checks(tildaM, 1);
0195             end
0196             indok = find(max(ide_hess.indno,[],1)==0);
0197 %             rhoM=sqrt(1-1./diag(inv(tildaM)));
0198 %             rhoM=(1-1./diag(inv(tildaM)));
0199             ind1=find(ide_hess.ind0);
0200             temp1=((MIM(ind1,ind1))\siH(:,ind1)');
0201             diag_chh=sum(siH(:,ind1)'.*temp1)';
0202 %             chh = siH(:,ind1)*((MIM(ind1,ind1))\siH(:,ind1)');
0203             ind1=ind1(ind1>offset);
0204             clre = siLRE(:,ind1-offset)*((MIM(ind1,ind1))\siLRE(:,ind1-offset)');
0205             if ~isempty(indok),
0206                 rhoM(indok)=sqrt(1./diag(inv(tildaM(indok,indok))));
0207                 normaliz(indok) = (sqrt(diag(inv(tildaM(indok,indok))))./deltaM(indok))'; %sqrt(diag(inv(MIM(indok,indok))))';
0208             end
0209             %             deltaM = deltaM.*abs(params')
0210             flag_score=0;
0211         end
0212         ide_strength_J(indok) = (1./(normaliz(indok)'./abs(params(indok)')));
0213         ide_strength_J_prior(indok) = (1./(normaliz(indok)'./normaliz1(indok)'));
0214         ide_strength_J(params==0)=ide_strength_J_prior(params==0);
0215         deltaM_prior = deltaM.*abs(normaliz1');
0216         deltaM = deltaM.*abs(params');
0217         deltaM(params==0)=deltaM_prior(params==0);
0218         quant = siJ./repmat(sqrt(diag(cmm)),1,nparam);
0219         siJnorm = vnorm(quant).*normaliz1;
0220         %                 siJnorm = vnorm(siJ(inok,:)).*normaliz;
0221         quant=[];
0222 %         inok = find((abs(TAU)<1.e-8));
0223 %         isok = find((abs(TAU)>=1.e-8));
0224 %         quant(isok,:) = siH(isok,:)./repmat(TAU(isok,1),1,nparam);
0225 %         quant(inok,:) = siH(inok,:)./repmat(mean(abs(TAU)),length(inok),nparam);
0226 %         quant = siH./repmat(sqrt(diag(chh)),1,nparam);
0227         quant = siH./repmat(sqrt(diag_chh),1,nparam);
0228         siHnorm = vnorm(quant).*normaliz1;
0229         %                 siHnorm = vnorm(siH./repmat(TAU,1,nparam)).*normaliz;
0230         quant=[];
0231 %         inok = find((abs(LRE)<1.e-8));
0232 %         isok = find((abs(LRE)>=1.e-8));
0233 %         quant(isok,:) = siLRE(isok,:)./repmat(LRE(isok,1),1,np);
0234 %         quant(inok,:) = siLRE(inok,:)./repmat(mean(abs(LRE)),length(inok),np);
0235         quant = siLRE./repmat(sqrt(diag(clre)),1,np);
0236         siLREnorm = vnorm(quant).*normaliz1(offset+1:end);
0237         %                 siLREnorm = vnorm(siLRE./repmat(LRE,1,nparam-offset)).*normaliz(offset+1:end);
0238         ide_hess.ide_strength_J=ide_strength_J; 
0239         ide_hess.ide_strength_J_prior=ide_strength_J_prior; 
0240         ide_hess.deltaM=deltaM; 
0241         ide_hess.deltaM_prior=deltaM_prior; 
0242         ide_moments.siJnorm=siJnorm; 
0243         ide_model.siHnorm=siHnorm; 
0244         ide_lre.siLREnorm=siLREnorm; 
0245         ide_hess.flag_score=flag_score; 
0246     end,
0247     if normalize_jacobians,
0248         normH = max(abs(siH)')';
0249         normH = normH(:,ones(nparam,1));
0250         normJ = max(abs(siJ)')';
0251         normJ = normJ(:,ones(nparam,1));
0252         normLRE = max(abs(siLRE)')';
0253         normLRE = normLRE(:,ones(size(gp,2),1));
0254     else
0255         normH = 1;
0256         normJ = 1;
0257         normLRE = 1;
0258     end
0259     ide_moments.indJJ=indJJ;
0260     ide_model.indH=indH;
0261     ide_lre.indLRE=indLRE;
0262     ide_moments.siJ=siJ;
0263     ide_model.siH=siH;
0264     ide_lre.siLRE=siLRE;
0265     ide_moments.GAM=GAM;
0266     ide_model.TAU=TAU;
0267     ide_lre.LRE=LRE;
0268 %     [ide_checks.idemodel_Mco, ide_checks.idemoments_Mco, ide_checks.idelre_Mco, ...
0269 %         ide_checks.idemodel_Pco, ide_checks.idemoments_Pco, ide_checks.idelre_Pco, ...
0270 %         ide_checks.idemodel_cond, ide_checks.idemoments_cond, ide_checks.idelre_cond, ...
0271 %         ide_checks.idemodel_ee, ide_checks.idemoments_ee, ide_checks.idelre_ee, ...
0272 %         ide_checks.idemodel_ind, ide_checks.idemoments_ind, ...
0273 %         ide_checks.idemodel_indno, ide_checks.idemoments_indno, ...
0274 %         ide_checks.idemodel_ino, ide_checks.idemoments_ino] = ...
0275 %         identification_checks(H(indH,:)./normH(:,ones(nparam,1)),JJ(indJJ,:)./normJ(:,ones(nparam,1)), gp(indLRE,:)./normLRE(:,ones(size(gp,2),1)));
0276     [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] = ...
0277         identification_checks(JJ(indJJ,:)./normJ, 0);
0278     [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] = ...
0279         identification_checks(H(indH,:)./normH, 0);
0280     [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] = ...
0281         identification_checks(gp(indLRE,:)./normLRE, 0);
0282     normJ=1;
0283     [U, S, V]=svd(JJ(indJJ,:)./normJ,0);
0284     S=diag(S);
0285     if nparam>8
0286         ide_moments.S = S([1:4, end-3:end]);
0287         ide_moments.V = V(:,[1:4, end-3:end]);
0288     else
0289         ide_moments.S = S;
0290         ide_moments.V = V;
0291     end
0292     
0293     indok = find(max(ide_moments.indno,[],1)==0);
0294     if advanced,
0295         [ide_moments.pars, ide_moments.cosnJ] = ident_bruteforce(JJ(indJJ,:)./normJ,max_dim_cova_group,options_.TeX,name_tex);
0296     end
0297 end