diff --git a/matlab/AHessian.m b/matlab/AHessian.m
index f1046c8dc..02cd7f0a3 100644
--- a/matlab/AHessian.m
+++ b/matlab/AHessian.m
@@ -50,7 +50,7 @@ end
 %         DOm = DR(:,:,ii)*Q*transpose(R) + R*DQ(:,:,ii)*transpose(R) + R*Q*transpose(DR(:,:,ii)); 
 %     end
     
-    while notsteady & t<smpl
+    while notsteady && t<smpl
         t  = t+1;
         v  = Y(:,t)-a(mf);
         F  = P(mf,mf) + H;
@@ -112,7 +112,7 @@ end
             a = T*(a+K*v);
         lik(t) = transpose(v)*iF*v;
         end
-        AHess = AHess + .5*(smpl+t0-1)*(vecDPmf' * kron(iF,iF) * vecDPmf);
+        AHess = AHess + .5*(smpl-t0+1)*(vecDPmf' * kron(iF,iF) * vecDPmf);
         if nargout > 1
         for ii = 1:k
 %             DLIK(ii,1)  = DLIK(ii,1) + (smpl-t0+1)*trace( iF*DF(:,:,ii) );
diff --git a/matlab/dsge_likelihood.m b/matlab/dsge_likelihood.m
index 8d3962832..b4bf85e4e 100644
--- a/matlab/dsge_likelihood.m
+++ b/matlab/dsge_likelihood.m
@@ -332,7 +332,7 @@ if (kalman_algo == 2) || (kalman_algo == 4)
     else
         if all(all(abs(H-diag(diag(H)))<1e-14))% ie, the covariance matrix is diagonal...
             H = diag(H);
-            mmm = mm;
+            mmm = mm; 
         else
             Z = [Z, eye(pp)];
             T = blkdiag(T,zeros(pp));
@@ -381,7 +381,7 @@ switch DynareOptions.lik_init
         % Use standard kalman filter except if the univariate filter is explicitely choosen.
     if kalman_algo == 0
         kalman_algo = 3;
-    elseif ~((kalman_algo == 3) || (kalman_algo == 4))
+    elseif ~((kalman_algo == 3) || (kalman_algo == 4)) 
             error(['diffuse filter: options_.kalman_algo can only be equal ' ...
                    'to 0 (default), 3 or 4'])
     end
@@ -455,7 +455,7 @@ switch DynareOptions.lik_init
         DynareOptions.lik_init = 1;
         Pstar = lyapunov_symm(T,R*Q*R',DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold);
     end
-    Pinf = [];
+    Pinf  = [];
     a = zeros(mm,1);
     Zflag = 0;
   otherwise
@@ -471,9 +471,13 @@ if analytic_derivation,
     no_DLIK = 0;
     full_Hess = analytic_derivation==2;
     asy_Hess = analytic_derivation==-2;
+    outer_product_gradient = analytic_derivation==-1;
     if asy_Hess,
         analytic_derivation=1;
     end
+    if outer_product_gradient,
+        analytic_derivation=1;
+    end
     DLIK = [];
     AHess = [];
     if nargin<8 || isempty(derivatives_info)
@@ -578,7 +582,7 @@ if analytic_derivation,
         end
     end
     if analytic_derivation==1,
-        analytic_deriv_info={analytic_derivation,DT,DYss,DOm,DH,DP};
+        analytic_deriv_info={analytic_derivation,DT,DYss,DOm,DH,DP,asy_Hess};
     else
         analytic_deriv_info={analytic_derivation,DT,DYss,DOm,DH,DP,D2T,D2Yss,D2Om,D2H,D2P};
     end
@@ -614,6 +618,8 @@ if ((kalman_algo==1) || (kalman_algo==3))% Multivariate Kalman Filter
     if analytic_derivation,
         LIK1=LIK;
         LIK=LIK1{1};
+        lik1=lik;
+        lik=lik1{1};
     end
     if isinf(LIK)
         if kalman_algo == 1
@@ -676,6 +682,8 @@ if (kalman_algo==2) || (kalman_algo==4)
     if analytic_derivation,
         LIK1=LIK;
         LIK=LIK1{1};
+        lik1=lik;
+        lik=lik1{1};
     end
     if DynareOptions.lik_init==3
         LIK = LIK+dLIK;
@@ -690,13 +698,17 @@ if analytic_derivation
         DLIK = LIK1{2};
         %                 [DLIK] = score(T,R,Q,H,Pstar,Y,DT,DYss,DOm,DH,DP,start,Z,kalman_tol,riccati_tol);
     end
-    if full_Hess,
+    if full_Hess ,
         Hess = -LIK1{3};
         %                     [Hess, DLL] = get_Hessian(T,R,Q,H,Pstar,Y,DT,DYss,DOm,DH,DP,D2T,D2Yss,D2Om,D2H,D2P,start,Z,kalman_tol,riccati_tol);
         %                     Hess0 = getHessian(Y,T,DT,D2T, R*Q*transpose(R),DOm,D2Om,Z,DYss,D2Yss);
     end
     if asy_Hess,
-        [Hess] = AHessian(T,R,Q,H,Pstar,Y,DT,DYss,DOm,DH,DP,start,Z,kalman_tol,riccati_tol);
+        if ~((kalman_algo==2) || (kalman_algo==4)),
+            [Hess] = AHessian(T,R,Q,H,Pstar,Y,DT,DYss,DOm,DH,DP,start,Z,kalman_tol,riccati_tol);
+        else
+        Hess = LIK1{3};
+        end
     end
 end
 
@@ -725,6 +737,11 @@ if analytic_derivation
     if no_DLIK==0
         DLIK = DLIK - dlnprior';
     end
+    if outer_product_gradient,
+        dlik = lik1{2};
+        dlik=[- dlnprior; dlik(start:end,:)];
+        Hess = dlik'*dlik;
+    end
 else
     lnprior = priordens(xparam1,BayesInfo.pshape,BayesInfo.p6,BayesInfo.p7,BayesInfo.p3,BayesInfo.p4);
 end
diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m
index 27c06afd7..5fdec1d78 100644
--- a/matlab/dynare_estimation_1.m
+++ b/matlab/dynare_estimation_1.m
@@ -213,22 +213,29 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
         else
             flag = 1;
         end
-        if ~exist('igg','var'),  % by M. Ratto
-            hh=[];
-            gg=[];
-            igg=[];
-        end   % by M. Ratto
         if isfield(options_,'ftol')
             crit = options_.ftol;
         else
             crit = 1.e-5;
         end
+        if options_.analytic_derivation,
+            analytic_grad=1;
+            ana_deriv = options_.analytic_derivation;
+            options_.analytic_derivation = -1;
+            crit = 1.e-7;
+            flag = 0;
+        else
+            analytic_grad=0;
+        end
         if isfield(options_,'nit')
             nit = options_.nit;
         else
             nit=1000;
         end
-        [xparam1,hh,gg,fval,invhess] = newrat(objective_function,xparam1,hh,gg,igg,crit,nit,flag,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
+        [xparam1,hh,gg,fval,invhess] = newrat(objective_function,xparam1,analytic_grad,crit,nit,flag,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
+        if options_.analytic_derivation,
+            options_.analytic_derivation = ana_deriv;
+        end
         parameter_names = bayestopt_.name;
         save([M_.fname '_mode.mat'],'xparam1','hh','gg','fval','invhess','parameter_names');
       case 6
diff --git a/matlab/dynare_identification.m b/matlab/dynare_identification.m
index d75299382..cc73b07d8 100644
--- a/matlab/dynare_identification.m
+++ b/matlab/dynare_identification.m
@@ -123,10 +123,13 @@ if isempty(dataset_),
     dataset_.info.varobs = options_.varobs;
     dataset_.rawdata = [];
     dataset_.missing.state = 0;
-    dataset_.missing.aindex = [];
+    for jdata=1:periods,
+        temp1{jdata}=[1:dataset_.info.nvobs]';
+    end
+    dataset_.missing.aindex = temp1;
     dataset_.missing.vindex = [];
     dataset_.missing.number_of_observations = [];
-    dataset_.missing.no_more_missing_observations = [];
+    dataset_.missing.no_more_missing_observations = 1;
     dataset_.descriptive.mean = [];
     dataset_.data = [];
 
diff --git a/matlab/identification_analysis.m b/matlab/identification_analysis.m
index baf3550c6..59846c1a8 100644
--- a/matlab/identification_analysis.m
+++ b/matlab/identification_analysis.m
@@ -7,7 +7,7 @@ function [ide_hess, ide_moments, ide_model, ide_lre, derivatives_info, info] = i
 %    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 data_info          [structure] data info for Kalman 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
diff --git a/matlab/kalman/likelihood/computeDLIK.m b/matlab/kalman/likelihood/computeDLIK.m
index bb2ea4c98..d08901356 100644
--- a/matlab/kalman/likelihood/computeDLIK.m
+++ b/matlab/kalman/likelihood/computeDLIK.m
@@ -24,18 +24,18 @@ persistent DK DF D2K D2F
 if notsteady
 if Zflag
     [DK,DF,DP1] = computeDKalmanZ(T,DT,DOm,P,DP,DH,Z,iF,K);
-    if nargout>3,
+    if nargout>4,
         [D2K,D2F,D2P1] = computeD2KalmanZ(T,DT,D2T,D2Om,P,DP,D2P,DH,Z,iF,K,DK);
     end
 else
     [DK,DF,DP1] = computeDKalman(T,DT,DOm,P,DP,DH,Z,iF,K);
-    if nargout>3,
+    if nargout>4,
         [D2K,D2F,D2P1] = computeD2Kalman(T,DT,D2T,D2Om,P,DP,D2P,DH,Z,iF,K,DK);
     end
 end
 else
     DP1=DP;
-    if nargout>3,
+    if nargout>4,
         D2P1=D2P;
     end
 end
@@ -95,10 +95,27 @@ for ii = 1:k
             end
         end
     end
+        
     Da(:,ii)   = DT(:,:,ii)*tmp + T*dtmp(:,ii);
     DLIK(ii,1)  = trace( iF*DF(:,:,ii) ) + 2*Dv(:,ii)'*iF*v - v'*(iF*DF(:,:,ii)*iF)*v;
 end
 
+if nargout==4,
+    %         Hesst(ii,jj) = getHesst_ij(v,Dv(:,ii),Dv(:,jj),0,iF,diFi,diFj,0,dFj,0);
+    vecDPmf = reshape(DF,[],k);
+    D2a = 2*Dv'*iF*Dv + (vecDPmf' * kron(iF,iF) * vecDPmf);
+%     for ii = 1:k
+%         
+%         diFi = -iF*DF(:,:,ii)*iF;
+%         for jj = 1:ii
+%             dFj    = DF(:,:,jj);
+%             diFj   = -iF*DF(:,:,jj)*iF;
+%             
+%             Hesst(ii,jj) = getHesst_ij(v*0,Dv(:,ii),Dv(:,jj),v*0,iF,diFi,diFj,0,-dFj,0);
+%         end
+%     end
+end
+
 % end of computeDLIK
 
 function Hesst_ij = getHesst_ij(e,dei,dej,d2eij,iS,diSi,diSj,d2iSij,dSj,d2Sij);
diff --git a/matlab/kalman/likelihood/kalman_filter.m b/matlab/kalman/likelihood/kalman_filter.m
index e6ab61615..46e56abc1 100644
--- a/matlab/kalman/likelihood/kalman_filter.m
+++ b/matlab/kalman/likelihood/kalman_filter.m
@@ -1,4 +1,4 @@
-function [LIK, likk, a, P] = kalman_filter(Y,start,last,a,P,kalman_tol,riccati_tol,presample,T,Q,R,H,Z,mm,pp,rr,Zflag,diffuse_periods,analytic_derivation,DT,DYss,DOm,DH,DP,D2T,D2Yss,D2Om,D2H,D2P)
+function [LIK, LIKK, a, P] = kalman_filter(Y,start,last,a,P,kalman_tol,riccati_tol,presample,T,Q,R,H,Z,mm,pp,rr,Zflag,diffuse_periods,analytic_derivation,DT,DYss,DOm,DH,DP,D2T,D2Yss,D2Om,D2H,D2P)
 % Computes the likelihood of a stationnary state space model.
 
 %@info:
@@ -123,6 +123,7 @@ LIK  = Inf;                % Default value of the log likelihood.
 oldK = Inf;
 notsteady   = 1;
 F_singular  = 1;
+asy_hess=0;
 
 if  analytic_derivation == 0,
     DLIK=[];
@@ -131,6 +132,7 @@ else
     k = size(DT,3);                                 % number of structural parameters
     DLIK  = zeros(k,1);                             % Initialization of the score.
     Da    = zeros(mm,k);                            % Derivative State vector.
+    dlikk = zeros(smpl,k);
     
     if Zflag==0,
         C = zeros(pp,mm);
@@ -144,12 +146,17 @@ else
         D2a    = zeros(mm,k,k);                             % State vector.
         d2C = zeros(pp,mm,k,k);
     else
+        asy_hess=D2T;
         Hess=[];
         D2a=[];
         D2T=[];
         D2Yss=[];
     end
+    if asy_hess,
+        Hess  = zeros(k,k);                             % Initialization of the Hessian
+    end
     LIK={inf,DLIK,Hess};
+    LIKK={likk,dlikk};
 end
 
 while notsteady && t<=last
@@ -185,14 +192,15 @@ while notsteady && t<=last
             if analytic_derivation==2,
                 [Da,DP,DLIKt,D2a,D2P, Hesst] = computeDLIK(k,tmp,Z,Zflag,v,T,K,P,iF,Da,DYss,DT,DOm,DP,DH,notsteady,D2a,D2Yss,D2T,D2Om,D2P);
             else
-                [Da,DP,DLIKt] = computeDLIK(k,tmp,Z,Zflag,v,T,K,P,iF,Da,DYss,DT,DOm,DP,DH,notsteady);
+                [Da,DP,DLIKt,Hesst] = computeDLIK(k,tmp,Z,Zflag,v,T,K,P,iF,Da,DYss,DT,DOm,DP,DH,notsteady);
             end
             if t>presample
                 DLIK = DLIK + DLIKt;
-                if analytic_derivation==2,
+                if analytic_derivation==2 || asy_hess,
                     Hess = Hess + Hesst;
                 end
             end
+            dlikk(s,:)=DLIKt;
         end
         a = T*tmp;
         P=Ptmp;
@@ -208,19 +216,31 @@ end
 
 % Add observation's densities constants and divide by two.
 likk(1:s) = .5*(likk(1:s) + pp*log(2*pi));
+if analytic_derivation,
+    DLIK = DLIK/2;
+    dlikk = dlikk/2;
+    if analytic_derivation==2 || asy_hess,
+        if asy_hess==0,
+        Hess = Hess + tril(Hess,-1)';
+        end
+        Hess = -Hess/2;
+    end
+end
 
 % Call steady state Kalman filter if needed.
 if t <= last
     if analytic_derivation,
         if analytic_derivation==2,
-            [tmp, likk(s+1:end)] = kalman_filter_ss(Y,t,last,a,T,K,iF,dF,Z,pp,Zflag, ...
+            [tmp, tmp2] = kalman_filter_ss(Y,t,last,a,T,K,iF,dF,Z,pp,Zflag, ...
                 analytic_derivation,Da,DT,DYss,D2a,D2T,D2Yss);
         else
-            [tmp, likk(s+1:end)] = kalman_filter_ss(Y,t,last,a,T,K,iF,dF,Z,pp,Zflag, ...
-                analytic_derivation,Da,DT,DYss);
+            [tmp, tmp2] = kalman_filter_ss(Y,t,last,a,T,K,iF,dF,Z,pp,Zflag, ...
+                analytic_derivation,Da,DT,DYss,asy_hess);
         end
+        likk(s+1:end)=tmp2{1};
+        dlikk(s+1:end,:)=tmp2{2};
         DLIK = DLIK + tmp{2};
-        if analytic_derivation==2,
+        if analytic_derivation==2 || asy_hess,
             Hess = Hess + tmp{3};
         end
     else
@@ -229,20 +249,19 @@ if t <= last
 end
 
 % Compute minus the log-likelihood.
-if presample
-    if presample>=diffuse_periods
-        likk = likk(1+(presample-diffuse_periods):end);
-    end
+if presample>diffuse_periods,
+    LIK = sum(likk(1+(presample-diffuse_periods):end));
+else
+    LIK = sum(likk);
 end
-LIK = sum(likk);
 
 if analytic_derivation,
-    DLIK = DLIK/2;
-    if analytic_derivation==2,
-        Hess = Hess + tril(Hess,-1)';
-        Hess = -Hess/2;
+    if analytic_derivation==2 || asy_hess,
         LIK={LIK, DLIK, Hess};
     else
         LIK={LIK, DLIK};
     end
+    LIKK={likk, dlikk};
+else
+    LIKK=likk;
 end
diff --git a/matlab/kalman/likelihood/kalman_filter_ss.m b/matlab/kalman/likelihood/kalman_filter_ss.m
index a99e51b3e..ccd8a87cc 100644
--- a/matlab/kalman/likelihood/kalman_filter_ss.m
+++ b/matlab/kalman/likelihood/kalman_filter_ss.m
@@ -81,6 +81,7 @@ t    = start;              % Initialization of the time index.
 likk = zeros(smpl,1);      % Initialization of the vector gathering the densities.
 LIK  = Inf;                % Default value of the log likelihood.
 notsteady = 0;
+asy_hess=0;
 if nargin<12
     analytic_derivation = 0;
 end
@@ -91,10 +92,16 @@ if  analytic_derivation == 0,
 else
     k = size(DT,3);                                 % number of structural parameters
     DLIK  = zeros(k,1);                             % Initialization of the score.
+    dlikk = zeros(smpl,k);
     if analytic_derivation==2,
         Hess  = zeros(k,k);                             % Initialization of the Hessian
     else
-        Hess=[];
+        asy_hess=D2a;
+        if asy_hess,
+            Hess  = zeros(k,k);                             % Initialization of the Hessian
+        else
+            Hess=[];
+        end
     end
 end
 
@@ -109,12 +116,13 @@ while t <= last
         if analytic_derivation==2,
             [Da,junk,DLIKt,D2a,junk2, Hesst] = computeDLIK(k,tmp,Z,Zflag,v,T,K,[],iF,Da,DYss,DT,[],[],[],notsteady,D2a,D2Yss,D2T,[],[]);
         else
-            [Da,junk,DLIKt] = computeDLIK(k,tmp,Z,Zflag,v,T,K,[],iF,Da,DYss,DT,[],[],[],notsteady);
+            [Da,junk,DLIKt,Hesst] = computeDLIK(k,tmp,Z,Zflag,v,T,K,[],iF,Da,DYss,DT,[],[],[],notsteady);
         end
         DLIK = DLIK + DLIKt;
-        if analytic_derivation==2,
+        if analytic_derivation==2 || asy_hess,
             Hess = Hess + Hesst;
         end
+        dlikk(t-start+1,:)=DLIKt;
     end
     a = T*tmp;
     likk(t-start+1) = transpose(v)*iF*v;
@@ -129,9 +137,17 @@ likk = .5*(likk + pp*log(2*pi));
 
 % Sum the observation's densities (minus the likelihood)
 LIK = sum(likk);
-if analytic_derivation==2,
-    LIK={LIK,DLIK,Hess};
+if analytic_derivation,
+    dlikk = dlikk/2;
+    DLIK = DLIK/2;
+    likk = {likk, dlikk};
 end
-if analytic_derivation==1,
+if analytic_derivation==2 || asy_hess,
+    if asy_hess==0,
+        Hess = Hess + tril(Hess,-1)';
+    end
+    Hess = -Hess/2;
+    LIK={LIK,DLIK,Hess};
+elseif analytic_derivation==1,
     LIK={LIK,DLIK};
 end
diff --git a/matlab/kalman/likelihood/univariate_computeDLIK.m b/matlab/kalman/likelihood/univariate_computeDLIK.m
index c2662ea63..a9eb34459 100644
--- a/matlab/kalman/likelihood/univariate_computeDLIK.m
+++ b/matlab/kalman/likelihood/univariate_computeDLIK.m
@@ -30,7 +30,7 @@ if notsteady,
             DF(j)=Z*DP(:,:,j)*Z'+DH;
             DK(:,j) = (DP(:,:,j)*Z')/F-PZ*DF(j)/F^2;
         end
-        if nargout>3
+        if nargout>4
             D2F = zeros(k,k);
             D2v = zeros(k,k);
             D2K = zeros(rows(K),k,k);
@@ -50,7 +50,7 @@ if notsteady,
         Dv   = -Da(Z,:) - DYss(Z,:);
         DF = squeeze(DP(Z,Z,:))+DH';
         DK = squeeze(DP(:,Z,:))/F-PZ*transpose(DF)/F^2;
-        if nargout>3
+        if nargout>4
             D2v   = squeeze(-D2a(Z,:,:) - D2Yss(Z,:,:));
             D2F = squeeze(D2P(Z,Z,:,:));
             D2K = squeeze(D2P(:,Z,:,:))/F;
@@ -60,7 +60,7 @@ if notsteady,
             end
         end
     end
-    if nargout>3
+    if nargout>4
         DD2K(:,indx,:,:)=D2K;
         DD2F(indx,:,:)=D2F;
     end
@@ -69,13 +69,13 @@ if notsteady,
 else
     DK = squeeze(DDK(:,indx,:));
     DF = DDF(indx,:)';
-    if nargout>3
+    if nargout>4
         D2K = squeeze(DD2K(:,indx,:,:));
         D2F = squeeze(DD2F(indx,:,:));
     end
     if Zflag
         Dv   = -Z*Da(:,:) - Z*DYss(:,:);
-        if nargout>3
+        if nargout>4
             D2v = zeros(k,k);
             for j=1:k,
                 D2v(:,j)   = -Z*D2a(:,:,j) - Z*D2Yss(:,:,j);
@@ -83,7 +83,7 @@ else
         end
     else
         Dv   = -Da(Z,:) - DYss(Z,:);
-        if nargout>3
+        if nargout>4
             D2v   = squeeze(-D2a(Z,:,:) - D2Yss(Z,:,:));
         end
     end
@@ -93,10 +93,15 @@ DLIK = DF/F + 2*Dv'/F*v - v^2/F^2*DF;
 if nargout==6
     Hesst = D2F/F-1/F^2*(DF*DF') + 2*D2v/F*v + 2*(Dv'*Dv)/F - 2*(DF*Dv)*v/F^2 ...
         - v^2/F^2*D2F - 2*v/F^2*(Dv'*DF') + 2*v^2/F^3*(DF*DF');
+elseif nargout==4,
+    D2a = 1/F^2*(DF*DF') + 2*(Dv'*Dv)/F ;
+%     D2a = -1/F^2*(DF*DF') + 2*(Dv'*Dv)/F  + 2*v^2/F^3*(DF*DF') ...
+%         - 2*(DF*Dv)*v/F^2 - 2*v/F^2*(Dv'*DF');
+%     D2a = +2*(Dv'*Dv)/F + (DF' * DF)/F^2;
 end
 
 Da = Da + DK*v+K*Dv;
-if nargout>3
+if nargout>4
     
     D2a = D2a + D2K*v;
     for j=1:k,
@@ -121,7 +126,7 @@ if notsteady,
             DP1(:,:,j)=DP(:,:,j) - (DP(:,Z,j))*K'-PZ*DK(:,j)';
         end
     end
-    if nargout>3,
+    if nargout>4,
         if Zflag,
             for j=1:k,
                 D2P = D2P;
diff --git a/matlab/kalman/likelihood/univariate_kalman_filter.m b/matlab/kalman/likelihood/univariate_kalman_filter.m
index d5efcb0a7..91b37ecfc 100644
--- a/matlab/kalman/likelihood/univariate_kalman_filter.m
+++ b/matlab/kalman/likelihood/univariate_kalman_filter.m
@@ -119,6 +119,7 @@ notsteady = 1;
 
 oldK = Inf;
 K = NaN(mm,pp);
+asy_hess=0;
 
 if  analytic_derivation == 0,
     DLIK=[];
@@ -127,6 +128,7 @@ else
     k = size(DT,3);                                 % number of structural parameters
     DLIK  = zeros(k,1);                             % Initialization of the score.
     Da    = zeros(mm,k);                            % Derivative State vector.
+    dlik  = zeros(smpl,k);
     
     if Zflag==0,
         C = zeros(pp,mm);
@@ -140,11 +142,15 @@ else
         D2a    = zeros(mm,k,k);                             % State vector.
         d2C = zeros(pp,mm,k,k);
     else
+        asy_hess=D2T;
         Hess=[];
         D2a=[];
         D2T=[];
         D2Yss=[];
     end
+    if asy_hess,
+        Hess  = zeros(k,k);                             % Initialization of the Hessian
+    end
     LIK={inf,DLIK,Hess};
 end
 
@@ -177,14 +183,15 @@ while notsteady && t<=last
                 if analytic_derivation==2,
                     [Da,DP,DLIKt,D2a,D2P, Hesst] = univariate_computeDLIK(k,i,z(i,:),Zflag,prediction_error,Ki,PZ,Fi,Da,DYss,DP,DH(d_index(i),:),notsteady,D2a,D2Yss,D2P);
                 else
-                    [Da,DP,DLIKt] = univariate_computeDLIK(k,i,z(i,:),Zflag,prediction_error,Ki,PZ,Fi,Da,DYss,DP,DH(d_index(i),:),notsteady);
+                    [Da,DP,DLIKt,Hesst] = univariate_computeDLIK(k,i,z(i,:),Zflag,prediction_error,Ki,PZ,Fi,Da,DYss,DP,DH(d_index(i),:),notsteady);
                 end
                 if t>presample
                     DLIK = DLIK + DLIKt;
-                    if analytic_derivation==2,
+                    if analytic_derivation==2 || asy_hess,
                         Hess = Hess + Hesst;
                     end
                 end
+                dlik(s,:)=DLIKt;
             end
             a = a + Ki*prediction_error;
             P = P - PZ*Ki';
@@ -208,19 +215,29 @@ end
 
 % Divide by two.
 lik(1:s) = .5*lik(1:s);
+if analytic_derivation,
+    DLIK = DLIK/2;
+    dlik = dlik/2;
+    if analytic_derivation==2 || asy_hess,
+%         Hess = (Hess + Hess')/2;
+        Hess = -Hess/2;
+    end
+end
 
 % Call steady state univariate kalman filter if needed.
 if t <= last
     if analytic_derivation,
         if analytic_derivation==2,
-            [tmp, lik(s+1:end)] = univariate_kalman_filter_ss(Y,t,last,a,P,kalman_tol,T,H,Z,pp,Zflag, ...
+            [tmp, tmp2] = univariate_kalman_filter_ss(Y,t,last,a,P,kalman_tol,T,H,Z,pp,Zflag, ...
                 analytic_derivation,Da,DT,DYss,DP,DH,D2a,D2T,D2Yss,D2P);
         else
-            [tmp, lik(s+1:end)] = univariate_kalman_filter_ss(Y,t,last,a,P,kalman_tol,T,H,Z,pp,Zflag, ...
-                analytic_derivation,Da,DT,DYss,DP,DH);
+            [tmp, tmp2] = univariate_kalman_filter_ss(Y,t,last,a,P,kalman_tol,T,H,Z,pp,Zflag, ...
+                analytic_derivation,Da,DT,DYss,DP,DH,asy_hess);
         end
+        lik(s+1:end)=tmp2{1};
+        dlik(s+1:end,:)=tmp2{2};
         DLIK = DLIK + tmp{2};
-        if analytic_derivation==2,
+        if analytic_derivation==2 || asy_hess,
             Hess = Hess + tmp{3};
         end
     else
@@ -236,11 +253,10 @@ else
 end
 
 if analytic_derivation,
-    DLIK = DLIK/2;
-    if analytic_derivation==2,
-        Hess = -Hess/2;
+    if analytic_derivation==2 || asy_hess,
         LIK={LIK, DLIK, Hess};
     else
         LIK={LIK, DLIK};
     end
+    lik={lik, dlik};
 end
diff --git a/matlab/kalman/likelihood/univariate_kalman_filter_ss.m b/matlab/kalman/likelihood/univariate_kalman_filter_ss.m
index 5747b0870..321451862 100644
--- a/matlab/kalman/likelihood/univariate_kalman_filter_ss.m
+++ b/matlab/kalman/likelihood/univariate_kalman_filter_ss.m
@@ -82,6 +82,7 @@ t    = start;              % Initialization of the time index.
 likk = zeros(smpl,1);      % Initialization of the vector gathering the densities.
 LIK  = Inf;                % Default value of the log likelihood.
 l2pi = log(2*pi);
+asy_hess=0;
 
 if nargin<12
     analytic_derivation = 0;
@@ -93,10 +94,16 @@ if  analytic_derivation == 0,
 else
     k = size(DT,3);                                 % number of structural parameters
     DLIK  = zeros(k,1);                             % Initialization of the score.
+    dlikk = zeros(smpl,k);
     if analytic_derivation==2,
         Hess  = zeros(k,k);                             % Initialization of the Hessian
     else
-        Hess=[];
+        asy_hess=D2a;
+        if asy_hess,
+            Hess  = zeros(k,k);                             % Initialization of the Hessian
+        else
+            Hess=[];
+        end
     end
 end
 
@@ -129,12 +136,13 @@ while t<=last
                 if analytic_derivation==2,
                     [Da,DPP,DLIKt,D2a,D2PP, Hesst] = univariate_computeDLIK(k,i,Z(i,:),Zflag,prediction_error,Ki,PPZ,Fi,Da,DYss,DPP,DH(i,:),0,D2a,D2Yss,D2PP);
                 else
-                    [Da,DPP,DLIKt] = univariate_computeDLIK(k,i,Z(i,:),Zflag,prediction_error,Ki,PPZ,Fi,Da,DYss,DPP,DH(i,:),0);
+                    [Da,DPP,DLIKt,Hesst] = univariate_computeDLIK(k,i,Z(i,:),Zflag,prediction_error,Ki,PPZ,Fi,Da,DYss,DPP,DH(i,:),0);
                 end
                 DLIK = DLIK + DLIKt;
-                if analytic_derivation==2,
+                if analytic_derivation==2 || asy_hess,
                     Hess = Hess + Hesst;
                 end
+                dlikk(s,:)=DLIKt;
             end
         end
     end
@@ -152,9 +160,15 @@ end
 likk = .5*likk;
 
 LIK = sum(likk);
-if analytic_derivation==2,
-    LIK={LIK,DLIK,Hess};
+if analytic_derivation,
+    dlikk = dlikk/2;
+    DLIK = DLIK/2;
+    likk = {likk, dlikk};
 end
-if analytic_derivation==1,
+if analytic_derivation==2 || asy_hess,
+%     Hess = (Hess + Hess')/2;
+    Hess = -Hess/2;
+    LIK={LIK,DLIK,Hess};
+elseif analytic_derivation==1,
     LIK={LIK,DLIK};
 end
\ No newline at end of file
diff --git a/matlab/mr_gstep.m b/matlab/mr_gstep.m
index 381b21fa0..a91d60ade 100644
--- a/matlab/mr_gstep.m
+++ b/matlab/mr_gstep.m
@@ -21,6 +21,10 @@ function [f0, x, ig] = mr_gstep(h1,x,func0,htol0,DynareDataset,DynareOptions,Mod
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
 n=size(x,1);
+if isempty(h1),
+    h1=DynareOptions.gradient_epsilon*ones(n,1);
+end
+
 
 if isempty(htol0)
     htol = 1.e-6;
diff --git a/matlab/newrat.m b/matlab/newrat.m
index 60c53f3b7..823de2618 100644
--- a/matlab/newrat.m
+++ b/matlab/newrat.m
@@ -1,4 +1,4 @@
-function [xparam1, hh, gg, fval, igg] = newrat(func0, x, hh, gg, igg, ftol0, nit, flagg, DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults)
+function [xparam1, hh, gg, fval, igg] = newrat(func0, x, analytic_derivation, ftol0, nit, flagg, DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults)
 %  [xparam1, hh, gg, fval, igg] = newrat(func0, x, hh, gg, igg, ftol0, nit, flagg, varargin)
 %
 %  Optimiser with outer product gradient and with sequences of univariate steps
@@ -10,9 +10,7 @@ function [xparam1, hh, gg, fval, igg] = newrat(func0, x, hh, gg, igg, ftol0, nit
 %    of the log-likelihood to compute outer product gradient
 %
 %  x = starting guess
-%  hh = initial Hessian [OPTIONAL]
-%  gg = initial gradient [OPTIONAL]
-%  igg = initial inverse Hessian [OPTIONAL]
+%  analytic_derivation = 1 if analytic derivs
 %  ftol0 = ending criterion for function change
 %  nit = maximum number of iterations
 %
@@ -56,14 +54,14 @@ gibbstol=length(BayesInfo.pshape)/50; %25;
 
 % func0 = str2func([func2str(func0),'_hh']);
 % func0 = func0;
-fval0=feval(func0,x,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
+[fval0,gg,hh]=feval(func0,x,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
 fval=fval0;
 
 % initialize mr_gstep and mr_hessian
-mr_hessian(1,x,[],[],[],DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
-outer_product_gradient=1;
 
+outer_product_gradient=1;
 if isempty(hh)
+    mr_hessian(1,x,[],[],[],DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
     [dum, gg, htol0, igg, hhg, h1]=mr_hessian(0,x,func0,flagit,htol,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
     if isempty(dum),
         outer_product_gradient=0;
@@ -85,6 +83,7 @@ else
     hh0=hh;
     hhg=hh;
     igg=inv(hh);
+    h1=[];
 end
 H = igg;
 disp(['Gradient norm ',num2str(norm(gg))])
@@ -125,7 +124,11 @@ while norm(gg)>gtol && check==0 && jit<nit
         if length(find(ig))<nx
             ggx=ggx*0;
             ggx(find(ig))=gg(find(ig));
-            hhx = reshape(dum,nx,nx);
+            if analytic_derivation,
+                hhx=hh;
+            else
+                hhx = reshape(dum,nx,nx);
+            end
             iggx=eye(length(gg));
             iggx(find(ig),find(ig)) = inv( hhx(find(ig),find(ig)) );
             [fvala,x0,fc,retcode] = csminit1(func0,x0,fval,ggx,0,iggx,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
@@ -159,6 +162,11 @@ while norm(gg)>gtol && check==0 && jit<nit
     if (fval0(icount)-fval)<ftol
         disp('No further improvement is possible!')
         check=1;
+        if analytic_derivation,
+            [fvalx,gg,hh]=feval(func0,xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
+            hhg=hh;
+            H = inv(hh);            
+        else
         if flagit==2
             hh=hh0;
         elseif flagg>0
@@ -173,6 +181,7 @@ while norm(gg)>gtol && check==0 && jit<nit
                 hh=hhg;
             end
         end
+        end
         disp(['Actual dxnorm ',num2str(norm(x(:,end)-x(:,end-1)))])
         disp(['FVAL          ',num2str(fval)])
         disp(['Improvement   ',num2str(fval0(icount)-fval)])
@@ -191,7 +200,7 @@ while norm(gg)>gtol && check==0 && jit<nit
         disp(['Ftol          ',num2str(ftol)])
         disp(['Htol          ',num2str(htol0)])
         htol=htol_base;
-        if norm(x(:,icount)-xparam1)>1.e-12
+        if norm(x(:,icount)-xparam1)>1.e-12 && analytic_derivation==0,
             try
                 save m1.mat x fval0 nig -append
             catch
@@ -225,6 +234,10 @@ while norm(gg)>gtol && check==0 && jit<nit
                 end
                 H = igg;
             end
+        else
+            [fvalx,gg,hh]=feval(func0,xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
+            hhg=hh;
+            H = inv(hh);
         end
         disp(['Gradient norm  ',num2str(norm(gg))])
         ee=eig(hh);