diff --git a/matlab/DsgeLikelihood.m b/matlab/DsgeLikelihood.m
index cf50dfadb..9a66a4632 100644
--- a/matlab/DsgeLikelihood.m
+++ b/matlab/DsgeLikelihood.m
@@ -200,8 +200,19 @@ Y   = data-trend;
 if analytic_derivation,
     if nargin<7 || isempty(derivatives_info)
         [A,B] = dynare_resolve;
+        if ~isempty(estim_params_.var_exo),
+            indexo=estim_params_.var_exo(:,1);
+        else
+            indexo=[];
+        end
+        if ~isempty(estim_params_.param_vals),
+            indparam=estim_params_.param_vals(:,1);;
+        else
+            indparam=[];
+        end
+            
         [dum, DT, DOm, DYss] = getH(A, B, M_,oo_,0, ...
-            estim_params_.param_vals(:,1),estim_params_.var_exo(:,1));
+            indparam,indexo);
     else
         DT = derivatives_info.DT;
         DOm = derivatives_info.DOm;
@@ -209,7 +220,7 @@ if analytic_derivation,
         clear derivatives_info,
     end
     iv = oo_.dr.restrict_var_list;
-    DYss = [zeros(length(DYss),offset) DYss];
+    DYss = [zeros(size(DYss,1),offset) DYss];
     DT = DT(iv,iv,:);
     DOm = DOm(iv,iv,:);
     DYss = DYss(iv,:);
diff --git a/matlab/dynare_identification.m b/matlab/dynare_identification.m
index d1125ac76..4521d938d 100644
--- a/matlab/dynare_identification.m
+++ b/matlab/dynare_identification.m
@@ -257,9 +257,9 @@ if iload <=0,
 
             if iteration,
                 [JJ, H, gam, gp, dA, dOm, dYss] = getJJ(A, B, M_,oo0,options_,0,indx,indexo,bayestopt_.mf2,nlags,useautocorr);
-                derivatives_info.dA=dA;
-                derivatives_info.dOm=dOm;
-                derivatives_info.dYss=dYss;
+                derivatives_info.DT=dA;
+                derivatives_info.DOm=dOm;
+                derivatives_info.DYss=dYss;
                 if BurninSampleSize == 0,
                     indJJ = (find(max(abs(JJ'))>1.e-8));
                     indH = (find(max(abs(H'))>1.e-8));
@@ -361,10 +361,11 @@ if iload <=0,
                         info = stoch_simul(options_.varobs);
                         datax=oo_.endo_simul(options_.varobs_id,100+1:end);
 %                         datax=data;
-                        [fval,cost_flag,ys,trend_coeff,info,DLIK,AHess] = DsgeLikelihood(params',data_info.gend,datax,data_info.data_index,data_info.number_of_observations,derivatives_info);
+                        [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);
                         cparam = inv(-AHess);
                         normaliz(indok) = sqrt(diag(cparam))';
                         cmm = siJ*((-AHess)\siJ');
+                        flag_score=1;
                     catch,
                         cmm = simulated_moment_uncertainty(indJJ, periods, replic);
                         %                 Jinv=(siJ(:,indok)'*siJ(:,indok))\siJ(:,indok)';
@@ -375,6 +376,7 @@ if iload <=0,
                         rhoM=sqrt(1-1./diag(inv(tildaM)));
                         deltaM = deltaM.*normaliz(indok)';
                         normaliz(indok) = sqrt(diag(inv(MIM)))';
+                        flag_score=0;
                     end
                     ide_strength_J(indok) = (1./(normaliz(indok)'./abs(params(indok)')));
                     ide_strength_J_prior(indok) = (1./(normaliz(indok)'./normaliz1(indok)'));
@@ -589,13 +591,15 @@ subplot(211)
 set(gca,'xlim',[0 nparam+1])
 set(gca,'xticklabel','')
 dy = get(gca,'ylim');
-% dy=dy(2)-dy(1);
 for ip=1:nparam,
     text(ip,dy(1),name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none')
 end
 legend('relative to param value','relative to prior std','Location','Best')
-title('Identification strength in the moments (log-scale)')
-
+if  flag_score,
+    title('Identification strength using likelihood scores (log-scale)')
+else
+    title('Identification strength in the moments (log-scale)')
+end
 subplot(212)
 % mmm = mean(siJmean);
 % [ss, is] = sort(mmm);