diff --git a/matlab/AHessian.m b/matlab/AHessian.m
index d03a018896cbce9189f9419d63b2cdf9bc7bb386..f1046c8dc08a9cb6e9240a8b63de7d8970131ed8 100644
--- a/matlab/AHessian.m
+++ b/matlab/AHessian.m
@@ -1,5 +1,5 @@
-function [AHess] = AHessian(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,start,mf,kalman_tol,riccati_tol)
-% function [AHess] = AHessian(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,start,mf,kalman_tol,riccati_tol)
+function [AHess, DLIK, LIK] = AHessian(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,start,mf,kalman_tol,riccati_tol)
+% function [AHess, DLIK, LIK] = AHessian(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,start,mf,kalman_tol,riccati_tol)
 %
 % computes the asymptotic hessian matrix of the log-likelihood function of
 % a state space model (notation as in kalman_filter.m in DYNARE
@@ -28,6 +28,7 @@ function [AHess] = AHessian(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,start,mf,kalman_tol,ri
 
     k = size(DT,3);                                 % number of structural parameters
     smpl = size(Y,2);                               % Sample size.
+    pp   = size(Y,1);                               % Maximum number of observed variables.
     mm   = size(T,2);                               % Number of state variables.
     a    = zeros(mm,1);                             % State vector.
     Om   = R*Q*transpose(R);                        % Variance of R times the vector of structural innovations.
@@ -36,6 +37,11 @@ function [AHess] = AHessian(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,start,mf,kalman_tol,ri
     notsteady   = 1;                                % Steady state flag.
     F_singular  = 1;
 
+lik  = zeros(smpl,1);                           % Initialization of the vector gathering the densities.
+LIK  = Inf;                                     % Default value of the log likelihood.
+if nargout > 1,
+    DLIK  = zeros(k,1);                             % Initialization of the score.
+end
     AHess  = zeros(k,k);                             % Initialization of the Hessian
     Da    = zeros(mm,k);                             % State vector.
     Dv = zeros(length(mf),k);
@@ -59,17 +65,21 @@ function [AHess] = AHessian(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,start,mf,kalman_tol,ri
             F_singular = 0;
             iF     = inv(F);
             K      = P(:,mf)*iF;
+            lik(t) = log(det(F))+transpose(v)*iF*v;
 
             [DK,DF,DP1] = computeDKalman(T,DT,DOm,P,DP,DH,mf,iF,K);
             
             		for ii = 1:k
                         Dv(:,ii)   = -Da(mf,ii) - DYss(mf,ii);
                         Da(:,ii)   = DT(:,:,ii)*(a+K*v) + T*(Da(:,ii)+DK(:,:,ii)*v + K*Dv(:,ii));
+                        if t>=start && nargout > 1
+                            DLIK(ii,1)  = DLIK(ii,1) + trace( iF*DF(:,:,ii) ) + 2*Dv(:,ii)'*iF*v - v'*(iF*DF(:,:,ii)*iF)*v;
+                        end
                     end
                     vecDPmf = reshape(DP(mf,mf,:),[],k);
-                    iPmf = inv(P(mf,mf));
+%                     iPmf = inv(P(mf,mf));
                     if t>=start
-                        AHess = AHess + Dv'*iPmf*Dv + .5*(vecDPmf' * kron(iPmf,iPmf) * vecDPmf);
+                        AHess = AHess + Dv'*iF*Dv + .5*(vecDPmf' * kron(iF,iF) * vecDPmf);
                     end
             a      = T*(a+K*v);                   
             P      = T*(P-K*P(mf,:))*transpose(T)+Om;
@@ -92,13 +102,23 @@ function [AHess] = AHessian(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,start,mf,kalman_tol,ri
                   	for ii = 1:k
                         Dv(:,ii)   = -Da(mf,ii)-DYss(mf,ii);
                         Da(:,ii)   = DT(:,:,ii)*(a+K*v) + T*(Da(:,ii)+DK(:,:,ii)*v + K*Dv(:,ii));
+                if t>=start && nargout >1
+                   DLIK(ii,1)  = DLIK(ii,1) + trace( iF*DF(:,:,ii) ) + 2*Dv(:,ii)'*iF*v - v'*(iF*DF(:,:,ii)*iF)*v;
+                end
                     end
              if t>=start
-                AHess = AHess + Dv'*iPmf*Dv; 
+                AHess = AHess + Dv'*iF*Dv; 
              end   
             a = T*(a+K*v);
+        lik(t) = transpose(v)*iF*v;
+        end
+        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) );
+        end
         end
-        AHess = AHess + .5*(smpl+t0-1)*(vecDPmf' * kron(iPmf,iPmf) * vecDPmf);
+        lik(t0:smpl) = lik(t0:smpl) + log(det(F));
 %         for ii = 1:k;
 %             for jj = 1:ii
 %              H(ii,jj) = trace(iPmf*(.5*DP(mf,mf,ii)*iPmf*DP(mf,mf,jj) + Dv(:,ii)*Dv(:,jj)'));
@@ -107,6 +127,13 @@ function [AHess] = AHessian(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,start,mf,kalman_tol,ri
     end    
     
 AHess = -AHess;  
+if nargout > 1,
+    DLIK = DLIK/2;
+end
+% adding log-likelihhod constants
+lik = (lik + pp*log(2*pi))/2;
+
+LIK = sum(lik(start:end)); % Minus the log-likelihood.
 % end of main function    
     
 function [DK,DF,DP1] = computeDKalman(T,DT,DOm,P,DP,DH,mf,iF,K)
diff --git a/matlab/DsgeLikelihood.m b/matlab/DsgeLikelihood.m
index bd09d245c34641e413c8f927baf85048dce522be..e2099fca1ca4c3ba7a29a19a528a8ebc594a3e84 100644
--- a/matlab/DsgeLikelihood.m
+++ b/matlab/DsgeLikelihood.m
@@ -439,9 +439,10 @@ end
 
 if analytic_derivation
     no_DLIK = 0;
+    full_Hess = 0;
     DLIK = [];
     AHess = [];
-    if nargin<7 || isempty(derivatives_info)
+    if nargin<8 || isempty(derivatives_info)
         [A,B,nou,nou,Model,DynareOptions,DynareResults] = dynare_resolve(Model,DynareOptions,DynareResults);
         if ~isempty(EstimatedParameters.var_exo)
             indexo=EstimatedParameters.var_exo(:,1);
@@ -453,12 +454,25 @@ if analytic_derivation
         else
             indparam=[];
         end
-        [dum, DT, DOm, DYss] = getH(A,B,Model,DynareResults,0,indparam,indexo);
+
+        if full_Hess,
+        [dum, DT, DOm, DYss, dum2, D2T, D2Om, D2Yss] = getH(A, B, Model,DynareResults,0,indparam,indexo);
+        else
+        [dum, DT, DOm, DYss] = getH(A, B, Model,DynareResults,0,indparam,indexo);
+        end
     else
         DT = derivatives_info.DT;
         DOm = derivatives_info.DOm;
         DYss = derivatives_info.DYss;
-        if isfield(derivatives_info,'no_DLIK')
+        if isfield(derivatives_info,'full_Hess'),
+            full_Hess = derivatives_info.full_Hess;
+        end
+        if full_Hess,
+        D2T = derivatives_info.D2T;
+        D2Om = derivatives_info.D2Om;
+        D2Yss = derivatives_info.D2Yss;
+        end
+        if isfield(derivatives_info,'no_DLIK'),
             no_DLIK = derivatives_info.no_DLIK;
         end
         clear('derivatives_info');
@@ -471,6 +485,17 @@ if analytic_derivation
     DH=zeros([size(H),length(xparam1)]);
     DQ=zeros([size(Q),length(xparam1)]);
     DP=zeros([size(T),length(xparam1)]);
+    if full_Hess,
+        for j=1:size(D2Yss,1),
+        tmp(j,:,:) = blkdiag(zeros(offset,offset), squeeze(D2Yss(j,:,:)));
+        end
+        D2Yss = tmp;
+        D2T = D2T(iv,iv,:,:);
+        D2Om = D2Om(iv,iv,:,:);
+        D2Yss = D2Yss(iv,:,:);
+        D2H=zeros([size(H),length(xparam1),length(xparam1)]);
+        D2P=zeros([size(T),length(xparam1),length(xparam1)]);
+    end
     for i=1:EstimatedParameters.nvx
         k =EstimatedParameters.var_exo(i,1);
         DQ(k,k,i) = 2*sqrt(Q(k,k));
@@ -478,11 +503,23 @@ if analytic_derivation
         kk = find(abs(dum) < 1e-12);
         dum(kk) = 0;
         DP(:,:,i)=dum;
+        if full_Hess
+        for j=1:i,
+            dum =  lyapunov_symm(T,D2Om(:,:,i,j),DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold);
+            kk = (abs(dum) < 1e-12);
+            dum(kk) = 0;
+            D2P(:,:,i,j)=dum;
+            D2P(:,:,j,i)=dum;
+        end
+        end
     end
     offset = EstimatedParameters.nvx;
     for i=1:EstimatedParameters.nvn
         k = EstimatedParameters.var_endo(i,1);
         DH(k,k,i+offset) = 2*sqrt(H(k,k));
+        if full_Hess
+        D2H(k,k,i+offset,i+offset) = 2;
+        end
     end
     offset = offset + EstimatedParameters.nvn;
     for j=1:EstimatedParameters.np
@@ -490,6 +527,21 @@ if analytic_derivation
         kk = find(abs(dum) < 1e-12);
         dum(kk) = 0;
         DP(:,:,j+offset)=dum;
+        if full_Hess
+        DTj = DT(:,:,j+offset);
+        DPj = dum;
+        for i=1:j,
+            DTi = DT(:,:,i+offset);
+            DPi = DP(:,:,i+offset);
+            D2Tij = D2T(:,:,i,j);
+            D2Omij = D2Om(:,:,i,j);
+            tmp = D2Tij*Pstar*T' + T*Pstar*D2Tij' + DTi*DPj*T' + DTj*DPi*T' + T*DPj*DTi' + T*DPi*DTj' + DTi*Pstar*DTj' + DTj*Pstar*DTi' + D2Omij;
+            dum = lyapunov_symm(T,tmp,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold);
+            dum(abs(dum)<1.e-12) = 0;
+            D2P(:,:,i+offset,j+offset) = dum;
+            D2P(:,:,j+offset,i+offset) = dum;
+        end
+        end
     end
 end
 
@@ -515,6 +567,10 @@ if ((kalman_algo==1) || (kalman_algo==3))% Multivariate Kalman Filter
             end
             if nargout==11
                 [AHess] = AHessian(T,R,Q,H,Pstar,Y,DT,DYss,DOm,DH,DP,start,Z,kalman_tol,riccati_tol);
+                if full_Hess,
+                    Hess = 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
             end
         end
     else
@@ -587,7 +643,19 @@ end
 % ------------------------------------------------------------------------------
 % 5. Adds prior if necessary
 % ------------------------------------------------------------------------------
-lnprior = priordens(xparam1,BayesInfo.pshape,BayesInfo.p6,BayesInfo.p7,BayesInfo.p3,BayesInfo.p4);
+if analytic_derivation
+    if full_Hess,
+        [lnprior, dlnprior, d2lnprior] = priordens(xparam1,BayesInfo.pshape,BayesInfo.p6,BayesInfo.p7,BayesInfo.p3,BayesInfo.p4);
+        AHess = Hess + d2lnprior;
+    else
+        [lnprior, dlnprior] = priordens(xparam1,BayesInfo.pshape,BayesInfo.p6,BayesInfo.p7,BayesInfo.p3,BayesInfo.p4);
+    end
+    if no_DLIK==0
+        DLIK = DLIK - dlnprior';
+    end
+else
+    lnprior = priordens(xparam1,BayesInfo.pshape,BayesInfo.p6,BayesInfo.p7,BayesInfo.p3,BayesInfo.p4);
+end
 fval    = (likelihood-lnprior);
 
 % Update DynareOptions.kalman_algo.
diff --git a/matlab/DsgeLikelihood_hh.m b/matlab/DsgeLikelihood_hh.m
index 896b474f6472cc538d6d4a552ff7afe7d2ac3861..3f0d0749ac6565b08f5a157c0dc346568dac9ba9 100644
--- a/matlab/DsgeLikelihood_hh.m
+++ b/matlab/DsgeLikelihood_hh.m
@@ -1,5 +1,5 @@
 function [fval,llik,cost_flag,ys,trend_coeff,info] = DsgeLikelihood_hh(xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults)
-% function [fval,cost_flag,ys,trend_coeff,info] = DsgeLikelihood(xparam1,gend,data,data_index,number_of_observations,no_more_missing_observations)
+% function [fval,llik,cost_flag,ys,trend_coeff,info] = DsgeLikelihood_hh(xparam1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults)
 % Evaluates the posterior kernel of a dsge model.
 %
 % INPUTS
@@ -11,6 +11,7 @@ function [fval,llik,cost_flag,ys,trend_coeff,info] = DsgeLikelihood_hh(xparam1,D
 %   no_more_missing_observations   [integer]
 % OUTPUTS
 %   fval        :     value of the posterior kernel at xparam1.
+%   llik        :     probabilities at each time point
 %   cost_flag   :     zero if the function returns a penalty, one otherwise.
 %   ys          :     steady state of original endogenous variables
 %   trend_coeff :
@@ -52,11 +53,14 @@ if nargin==1
     return
 end
 
-fval          = [];
-ys            = [];
-trend_coeff   = [];
-cost_flag     = 1;
+% Initialization of the returned variables and others...
+fval        = [];
+ys          = [];
+trend_coeff = [];
+cost_flag   = 1;
 llik=NaN;
+info        = 0;
+singularity_flag = 0;
 
 if DynareOptions.block == 1
     error('DsgeLikelihood_hh:: This routine (called if mode_compute==5) is not compatible with the block option!')
@@ -159,8 +163,6 @@ end
 Model.Sigma_e = Q;
 Model.H = H;
 
-
-
 %------------------------------------------------------------------------------
 % 2. call model setup & reduction program
 %------------------------------------------------------------------------------
@@ -223,8 +225,31 @@ Y   = DynareDataset.data-trend;
 %------------------------------------------------------------------------------
 % 3. Initial condition of the Kalman filter
 %------------------------------------------------------------------------------
-
 kalman_algo = DynareOptions.kalman_algo;
+
+% resetting measurement errors covariance matrix for univariate filters
+if (kalman_algo == 2) || (kalman_algo == 4)
+    if isequal(H,0)
+        H = zeros(nobs,1);
+        mmm = mm;
+    else
+        if all(all(abs(H-diag(diag(H)))<1e-14))% ie, the covariance matrix is diagonal...
+            H = diag(H);
+            mmm = mm; 
+        else
+            Z = [Z, eye(pp)];
+            T = blkdiag(T,zeros(pp));
+            Q = blkdiag(Q,H);
+            R = blkdiag(R,eye(pp));
+            Pstar = blkdiag(Pstar,H);
+            Pinf  = blckdiag(Pinf,zeros(pp));
+            H = zeros(nobs,1);
+            mmm   = mm+pp;
+        end
+    end
+end
+
+
 diffuse_periods = 0;
 switch DynareOptions.lik_init
   case 1% Standard initialization with the steady state of the state equation.
@@ -271,20 +296,32 @@ switch DynareOptions.lik_init
         if isinf(dLIK)
             % Go to univariate diffuse filter if singularity problem.
             kalman_algo = 4;
+            singularity_flag = 1;
         end
     end
     if (kalman_algo==4)
         % Univariate Diffuse Kalman Filter
-        if no_correlation_flag
-            mmm = mm;
-        else
-            Z = [Z, eye(pp)];
-            T = blkdiag(T,zeros(pp));
-            Q = blkdiag(Q,H);
-            R = blkdiag(R,eye(pp));
-            Pstar = blkdiag(Pstar,H);
-            Pinf  = blckdiag(Pinf,zeros(pp));
-            mmm   = mm+pp;
+        if singularity_flag
+            if isequal(H,0)
+                H = zeros(nobs,1);
+                mmm = mm;
+            else
+                if all(all(abs(H-diag(diag(H)))<1e-14))% ie, the covariance matrix is diagonal...
+                    H = diag(H);
+                    mmm = mm; 
+                else
+                    Z = [Z, eye(pp)];
+                    T = blkdiag(T,zeros(pp));
+                    Q = blkdiag(Q,H);
+                    R = blkdiag(R,eye(pp));
+                    Pstar = blkdiag(Pstar,H);
+                    Pinf  = blckdiag(Pinf,zeros(pp));
+                    H = zeros(nobs,1);
+                    mmm   = mm+pp;
+                end
+            end
+            % no need to test again for correlation elements
+            singularity_flag = 0;
         end
         [dLIK,dlik,a,Pstar] = univariate_kalman_filter_d(DynareDataset.missing.aindex,DynareDataset.missing.number_of_observations,DynareDataset.missing.no_more_missing_observations, ...
                                                               Y, 1, size(Y,2), ...
@@ -294,13 +331,13 @@ switch DynareOptions.lik_init
         diffuse_periods = length(dlik);
     end
   case 4% Start from the solution of the Riccati equation.
-        if kalman_algo ~= 2
+    if kalman_algo ~= 2
         kalman_algo = 1;
     end
     if isequal(H,0)
-        [err,Pstar] = kalman_steady_state(transpose(T),R*Q*transpose(R),transpose(build_selection_matrix(mf,np,length(mf))));
+        [err,Pstar] = kalman_steady_state(transpose(T),R*Q*transpose(R),transpose(build_selection_matrix(Z,np,length(Z))));
     else
-        [err,Pstar] = kalman_steady_state(transpose(T),R*Q*transpose(R),transpose(build_selection_matrix(mf,np,length(mf))),H);
+        [err,Pstar] = kalman_steady_state(transpose(T),R*Q*transpose(R),transpose(build_selection_matrix(Z,np,length(Z))),H);
     end
     if err
         disp(['DsgeLikelihood:: I am not able to solve the Riccati equation, so I switch to lik_init=1!']);
@@ -316,8 +353,6 @@ end
 % 4. Likelihood evaluation
 %------------------------------------------------------------------------------
 
-singularity_flag = 0;
-
 if ((kalman_algo==1) || (kalman_algo==3))% Multivariate Kalman Filter
     if no_missing_data_flag
         [LIK,lik] = kalman_filter(Y,diffuse_periods+1,size(Y,2), ...
@@ -333,6 +368,11 @@ if ((kalman_algo==1) || (kalman_algo==3))% Multivariate Kalman Filter
                                                T,Q,R,H,Z,mm,pp,rr,Zflag,diffuse_periods);
     end
     if isinf(LIK)
+        if kalman_algo == 1
+            kalman_algo = 2;
+        else
+            kalman_algo = 4;
+        end
         singularity_flag = 1;
     else
         if DynareOptions.lik_init==3
@@ -342,21 +382,30 @@ if ((kalman_algo==1) || (kalman_algo==3))% Multivariate Kalman Filter
     end
 end
 
-if ( (singularity_flag) || (kalman_algo==2) || (kalman_algo==4) )% Univariate Kalman Filter
+if ( singularity_flag || (kalman_algo==2) || (kalman_algo==4) )
+    % Univariate Kalman Filter
+    % resetting measurement error covariance matrix when necessary                                                           % 
     if singularity_flag
-        if no_correlation
+        if isequal(H,0)
+            H = zeros(nobs,1);
             mmm = mm;
         else
-            Z = [Z, eye(pp)];
-            T = blkdiag(T,zeros(pp));
-            Q = blkdiag(Q,H);
-            R = blkdiag(R,eye(pp));
-            Pstar = blkdiag(Pstar,H);
-            Pinf  = blckdiag(Pinf,zeros(pp));
-            mmm   = mm+pp;
-            a = [a; zeros(pp,1)];
+            if all(all(abs(H-diag(diag(H)))<1e-14))% ie, the covariance matrix is diagonal...
+                H = diag(H);
+                mmm = mm;
+            else
+                Z = [Z, eye(pp)];
+                T = blkdiag(T,zeros(pp));
+                Q = blkdiag(Q,H);
+                R = blkdiag(R,eye(pp));
+                Pstar = blkdiag(Pstar,H);
+                Pinf  = blckdiag(Pinf,zeros(pp));
+                H = zeros(nobs,1);
+                mmm   = mm+pp;
+            end
         end
     end
+
     [LIK,lik] = univariate_kalman_filter(DynareDataset.missing.aindex,DynareDataset.missing.number_of_observations,DynareDataset.missing.no_more_missing_observations,Y,diffuse_periods+1,size(Y,2), ...
                                        a,Pstar, ...
                                        DynareOptions.kalman_tol, ...
diff --git a/matlab/disp_identification.m b/matlab/disp_identification.m
index 71c226f8a1fe729ea3189749cd42afcd01463f66..057981de8862baaac86dea86e3f6d4e8afd5007d 100644
--- a/matlab/disp_identification.m
+++ b/matlab/disp_identification.m
@@ -188,7 +188,7 @@ if any(idemoments.ino),
         iweak = length(find(idemoments.jweak_pair(:,j)));
         if iweak,
             [jx,jy]=find(jmap_pair==j);
-            jstore=[jstore  jx(1) jy(1)]';
+            jstore=[jstore'  jx(1) jy(1)]';
             if SampleSize > 1
                 disp(['    [',name{jx(1)},',',name{jy(1)},'] are PAIRWISE collinear (with tol = 1.e-10) for ',num2str((iweak)/SampleSize*100),'% of MC runs!' ])
             else
diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m
index 0ec5dc77928cec59537cc4a9bbc6159565aab48d..8eec84276eb3ee5080e1b75c5f54c0e2d950f0d8 100644
--- a/matlab/dynare_estimation_1.m
+++ b/matlab/dynare_estimation_1.m
@@ -353,7 +353,7 @@ if ~options_.mh_posterior_mode_estimation && options_.cova_compute
 end
 
 if options_.mode_check == 1 && ~options_.mh_posterior_mode_estimation
-    mode_check('objective_function',xparam1,hh,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
+    mode_check(objective_function,xparam1,hh,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
 end
 
 if ~options_.mh_posterior_mode_estimation
diff --git a/matlab/evaluate_steady_state_file.m b/matlab/evaluate_steady_state_file.m
index 9671ff06db59913e4f83b6dfc7639783b8331d0b..5ee03c0b51eff84982f2fd86641930098eea5910 100644
--- a/matlab/evaluate_steady_state_file.m
+++ b/matlab/evaluate_steady_state_file.m
@@ -41,6 +41,8 @@ function [ys,params1,check] = evaluate_steady_state_file(ys_init,exo_ss,params,f
 
     if steadystate_flag == 1
             % old format
+            assignin('base','tmp_00_',params);
+            evalin('base','M_.params=tmp_00_; clear(''tmp_00_'')');
             h_steadystate = str2func([fname '_steadystate']);                       
             [ys,check] = h_steadystate(ys_init, exo_ss);
             params1 = evalin('base','M_.params');
diff --git a/matlab/getH.m b/matlab/getH.m
index df7b1cea35a8973d4b2b57a57120829b81ea73fe..a457969bb3659c96e152b108f858b5481b599d1b 100644
--- a/matlab/getH.m
+++ b/matlab/getH.m
@@ -1,4 +1,4 @@
-function [H, dA, dOm, Hss, gp, info] = getH(A, B, M_,oo_,kronflag,indx,indexo)
+function [H, dA, dOm, Hss, gp, d2A, d2Om, H2ss] = getH(A, B, M_,oo_,kronflag,indx,indexo)
 
 % computes derivative of reduced form linear model w.r.t. deep params
 %
@@ -31,10 +31,29 @@ yy0=oo_.dr.ys(I);
 %     yy0 = [ yy0; oo_.dr.ys(find(M_.lead_lag_incidence(j,:)))];
 % end
 dyssdtheta=zeros(length(oo_.dr.ys),M_.param_nbr);
-df = feval([M_.fname,'_params_derivs'],yy0, oo_.exo_steady_state', ...
-           M_.params, oo_.dr.ys, 1, dyssdtheta);
-[residual, gg1] = feval([M_.fname,'_static'],oo_.dr.ys, oo_.exo_steady_state', M_.params);
+if nargout>5,
+    [residual, gg1, gg2] = feval([M_.fname,'_static'],oo_.dr.ys, oo_.exo_steady_state', M_.params);
+    [df, gp, d2f] = feval([M_.fname,'_params_derivs'],yy0, oo_.exo_steady_state', ...
+        M_.params, oo_.dr.ys, 1, dyssdtheta);
+    d2f = get_all_resid_2nd_derivs(d2f,length(oo_.dr.ys),M_.param_nbr);
+    d2f = d2f(:,indx,indx);
+    if isempty(find(gg2)),
+        for j=1:length(indx),
+        d2yssdtheta(:,:,j) = -gg1\d2f(:,:,j);
+        end
+    else
+        gam = d2f*0;
+        for j=1:length(indx),
+        d2yssdtheta(:,:,j) = -gg1\(d2f(:,:,j)+gam(:,:,j));
+        end
+    end
+else
+    [residual, gg1] = feval([M_.fname,'_static'],oo_.dr.ys, oo_.exo_steady_state', M_.params);
+    df = feval([M_.fname,'_params_derivs'],yy0, oo_.exo_steady_state', ...
+        M_.params, oo_.dr.ys, 1, dyssdtheta);
+end
 dyssdtheta = -gg1\df;
+
 if any(any(isnan(dyssdtheta))),    
     [U,T] = schur(gg1);
     global options_
@@ -45,15 +64,29 @@ if any(any(isnan(dyssdtheta))),
     [U,T] = ordschur(U,T,e1);
     T = T(k+1:end,k+1:end);
     dyssdtheta = -U(:,k+1:end)*(T\U(:,k+1:end)')*df;
+    if nargout>5,
+        for j=1:length(indx),
+            d2yssdtheta(:,:,j) = -U(:,k+1:end)*(T\U(:,k+1:end)')*d2f(:,:,j);
+        end
+    end
 end
+if nargout>5,
+[df, gp, d2f, gpp] = feval([M_.fname,'_params_derivs'],yy0, oo_.exo_steady_state', ...
+           M_.params, oo_.dr.ys, 1, dyssdtheta);
+H2ss = d2yssdtheta(oo_.dr.order_var,:,:);
+[residual, g1, g2, g3] = feval([M_.fname,'_dynamic'],yy0, oo_.exo_steady_state', ...
+                            M_.params, oo_.dr.ys, 1);
+else
 [df, gp] = feval([M_.fname,'_params_derivs'],yy0, oo_.exo_steady_state', ...
            M_.params, oo_.dr.ys, 1, dyssdtheta);
 [residual, g1, g2 ] = feval([M_.fname,'_dynamic'],yy0, oo_.exo_steady_state', ...
                             M_.params, oo_.dr.ys, 1);
+end
 
 Hss = dyssdtheta(oo_.dr.order_var,indx);
 dyssdtheta = dyssdtheta(I,:);
 [nr, nc]=size(g2);
+[m, nelem]=size(g1);
 nc = sqrt(nc);
 ns = max(max(M_.lead_lag_incidence));
 gp2 = gp*0;
@@ -73,19 +106,9 @@ end
 gp = gp+gp2;
 gp = gp(:,:,indx);
 param_nbr = length(indx);
-
-% order_var = [oo_.dr.order_var; ...
-%     [size(oo_dr.ghx,2)+1:size(oo_dr.ghx,2)+size(oo_.dr.transition_auxiliary_variables,1)]' ];
-% [A(order_var,order_var),B(order_var,:)]=dynare_resolve;
-% [A,B,ys,info]=dynare_resolve;
-%   if info(1) > 0
-%     H = [];
-%     A0 = [];
-%     B0 = [];
-%     dA = [];
-%     dOm = [];
-%     return
-%   end
+if nargout>5,
+    param_nbr_2 = param_nbr*(param_nbr+1)/2;
+end
 
 m = size(A,1);
 n = size(B,2);
@@ -93,9 +116,15 @@ n = size(B,2);
 
 
 klen = M_.maximum_endo_lag + M_.maximum_endo_lead + 1;
-k1 = M_.lead_lag_incidence(find([1:klen] ~= M_.maximum_endo_lag+1),:);
-a = g1(:,nonzeros(k1'));
-da = gp(:,nonzeros(k1'),:);
+k11 = M_.lead_lag_incidence(find([1:klen] ~= M_.maximum_endo_lag+1),:);
+a = g1(:,nonzeros(k11'));
+da = gp(:,nonzeros(k11'),:);
+if nargout > 5,
+    nelem = size(g1,2);
+    g22 = get_all_2nd_derivs(gpp,m,nelem,M_.param_nbr);
+    g22 = g22(:,:,indx,indx);
+    d2a = g22(:,nonzeros(k11'),:,:);
+end
 kstate = oo_.dr.kstate;
 
 GAM1 = zeros(M_.endo_nbr,M_.endo_nbr);
@@ -103,22 +132,26 @@ Dg1 = zeros(M_.endo_nbr,M_.endo_nbr,param_nbr);
 % nf = find(M_.lead_lag_incidence(M_.maximum_endo_lag+2,:));
 % GAM1(:, nf) = -g1(:,M_.lead_lag_incidence(M_.maximum_endo_lag+2,nf));
 
-k = find(kstate(:,2) == M_.maximum_endo_lag+2 & kstate(:,3));
-GAM1(:, kstate(k,1)) = -a(:,kstate(k,3));
-Dg1(:, kstate(k,1), :) = -da(:,kstate(k,3),:);
-k = find(kstate(:,2) > M_.maximum_endo_lag+2 & kstate(:,3));
-kk = find(kstate(:,2) > M_.maximum_endo_lag+2 & ~kstate(:,3));
-nauxe = 0;
-if ~isempty(k),
-    ax(:,k)= -a(:,kstate(k,3));
-    ax(:,kk)= 0;
-    dax(:,k,:)= -da(:,kstate(k,3),:);
-    dax(:,kk,:)= 0;
-    nauxe=size(ax,2);
-    GAM1 = [ax GAM1];
-    Dg1 = cat(2,dax,Dg1);
-    clear ax
+k1 = find(kstate(:,2) == M_.maximum_endo_lag+2 & kstate(:,3));
+GAM1(:, kstate(k1,1)) = -a(:,kstate(k1,3));
+Dg1(:, kstate(k1,1), :) = -da(:,kstate(k1,3),:);
+if nargout > 5,
+    D2g1 = zeros(M_.endo_nbr,M_.endo_nbr,param_nbr,param_nbr);
+    D2g1(:, kstate(k1,1), :, :) = -d2a(:,kstate(k1,3),:,:);
 end
+% k = find(kstate(:,2) > M_.maximum_endo_lag+2 & kstate(:,3));
+% kk = find(kstate(:,2) > M_.maximum_endo_lag+2 & ~kstate(:,3));
+% nauxe = 0;
+% if ~isempty(k),
+%     ax(:,k)= -a(:,kstate(k,3));
+%     ax(:,kk)= 0;
+%     dax(:,k,:)= -da(:,kstate(k,3),:);
+%     dax(:,kk,:)= 0;
+%     nauxe=size(ax,2);
+%     GAM1 = [ax GAM1];
+%     Dg1 = cat(2,dax,Dg1);
+%     clear ax
+% end
 
 
 [junk,cols_b,cols_j] = find(M_.lead_lag_incidence(M_.maximum_endo_lag+1, ...
@@ -127,84 +160,95 @@ GAM0 = zeros(M_.endo_nbr,M_.endo_nbr);
 Dg0 = zeros(M_.endo_nbr,M_.endo_nbr,param_nbr);
 GAM0(:,cols_b) = g1(:,cols_j);
 Dg0(:,cols_b,:) = gp(:,cols_j,:);
+if nargout > 5,
+    D2g0 = zeros(M_.endo_nbr,M_.endo_nbr,param_nbr,param_nbr);
+    D2g0(:, cols_b, :, :) = g22(:,cols_j,:,:);
+end
 % GAM0 = g1(:,M_.lead_lag_incidence(M_.maximum_endo_lag+1,:));
 
 
-k = find(kstate(:,2) == M_.maximum_endo_lag+1 & kstate(:,4));
+k2 = find(kstate(:,2) == M_.maximum_endo_lag+1 & kstate(:,4));
 GAM2 = zeros(M_.endo_nbr,M_.endo_nbr);
 Dg2 = zeros(M_.endo_nbr,M_.endo_nbr,param_nbr);
 % nb = find(M_.lead_lag_incidence(1,:));
 % GAM2(:, nb) = -g1(:,M_.lead_lag_incidence(1,nb));
-GAM2(:, kstate(k,1)) = -a(:,kstate(k,4));
-Dg2(:, kstate(k,1), :) = -da(:,kstate(k,4),:);
-naux = 0;
-k = find(kstate(:,2) < M_.maximum_endo_lag+1 & kstate(:,4));
-kk = find(kstate(:,2) < M_.maximum_endo_lag+1 );
-if ~isempty(k),
-    ax(:,k)= -a(:,kstate(k,4));
-    ax = ax(:,kk);
-    dax(:,k,:)= -da(:,kstate(k,4),:);
-    dax = dax(:,kk,:);
-    naux = size(ax,2);
-    GAM2 = [GAM2 ax];
-    Dg2 = cat(2,Dg2,dax);
-end
-
-GAM0 = blkdiag(GAM0,eye(naux));
-Dg0 = cat(1,Dg0,zeros(naux+nauxe,M_.endo_nbr,param_nbr));
-Dg0 = cat(2,Dg0,zeros(m+nauxe,naux,param_nbr));
-Dg0 = cat(2,zeros(m+nauxe,nauxe,param_nbr),Dg0);
-
-GAM2 = [GAM2 ; A(M_.endo_nbr+1:end,:)];
-Dg2 = cat(1,Dg2,zeros(naux+nauxe,m,param_nbr));
-Dg2 = cat(2,zeros(m+nauxe,nauxe,param_nbr),Dg2);
-GAM2 = [zeros(m+nauxe,nauxe) [GAM2; zeros(nauxe,m)]];
-GAM0 = [[zeros(m,nauxe);(eye(nauxe))] [GAM0; zeros(nauxe,m)]];
+GAM2(:, kstate(k2,1)) = -a(:,kstate(k2,4));
+Dg2(:, kstate(k2,1), :) = -da(:,kstate(k2,4),:);
+if nargout > 5,
+    D2g2 = zeros(M_.endo_nbr,M_.endo_nbr,param_nbr,param_nbr);
+    D2g2(:, kstate(k2,1), :, :) = -d2a(:,kstate(k2,4),:,:);
+end% naux = 0;
+% k = find(kstate(:,2) < M_.maximum_endo_lag+1 & kstate(:,4));
+% kk = find(kstate(:,2) < M_.maximum_endo_lag+1 );
+% if ~isempty(k),
+%     ax(:,k)= -a(:,kstate(k,4));
+%     ax = ax(:,kk);
+%     dax(:,k,:)= -da(:,kstate(k,4),:);
+%     dax = dax(:,kk,:);
+%     naux = size(ax,2);
+%     GAM2 = [GAM2 ax];
+%     Dg2 = cat(2,Dg2,dax);
+% end
+% 
+% GAM0 = blkdiag(GAM0,eye(naux));
+% Dg0 = cat(1,Dg0,zeros(naux+nauxe,M_.endo_nbr,param_nbr));
+% Dg0 = cat(2,Dg0,zeros(m+nauxe,naux,param_nbr));
+% Dg0 = cat(2,zeros(m+nauxe,nauxe,param_nbr),Dg0);
+% 
+% GAM2 = [GAM2 ; A(M_.endo_nbr+1:end,:)];
+% Dg2 = cat(1,Dg2,zeros(naux+nauxe,m,param_nbr));
+% Dg2 = cat(2,zeros(m+nauxe,nauxe,param_nbr),Dg2);
+% GAM2 = [zeros(m+nauxe,nauxe) [GAM2; zeros(nauxe,m)]];
+% GAM0 = [[zeros(m,nauxe);(eye(nauxe))] [GAM0; zeros(nauxe,m)]];
 
 GAM3 = -g1(:,length(yy0)+1:end);
 % GAM3 = -g1(oo_.dr.order_var,length(yy0)+1:end);
-GAM3 = [GAM3; zeros(naux+nauxe,size(GAM3,2))];
+% GAM3 = [GAM3; zeros(naux+nauxe,size(GAM3,2))];
 % Dg3 = -gp(oo_.dr.order_var,length(yy0)+1:end,:);
 Dg3 = -gp(:,length(yy0)+1:end,:);
-Dg3 = cat(1,Dg3,zeros(naux+nauxe,size(GAM3,2),param_nbr));
-
-auxe = zeros(0,1);
-k0 = kstate(find(kstate(:,2) >= M_.maximum_endo_lag+2),:);;
-i0 = find(k0(:,2) == M_.maximum_endo_lag+2);
-for i=M_.maximum_endo_lag+3:M_.maximum_endo_lag+2+M_.maximum_endo_lead,
-    i1 = find(k0(:,2) == i);
-    n1 = size(i1,1);
-    j = zeros(n1,1);
-    for j1 = 1:n1
-        j(j1) = find(k0(i0,1)==k0(i1(j1),1));
-    end
-    auxe = [auxe; i0(j)];
-    i0 = i1;
+if nargout>5,
+    D2g3 = -g22(:,length(yy0)+1:end,:,:);
+    clear g22 d2a
 end
-auxe = [(1:size(auxe,1))' auxe(end:-1:1)];
-n_ir1 = size(auxe,1);
-nr = m + n_ir1;
-GAM1 = [[GAM1 zeros(size(GAM1,1),naux)]; zeros(naux+nauxe,m+nauxe)];
-GAM1(m+1:end,:) = sparse(auxe(:,1),auxe(:,2),ones(n_ir1,1),n_ir1,nr);
-Dg1 = cat(2,Dg1,zeros(M_.endo_nbr,naux,param_nbr));
-Dg1 = cat(1,Dg1,zeros(naux+nauxe,m+nauxe,param_nbr));
-
-Ax = A;
-A1 = A;
-Bx = B;
-B1 = B;
-for j=1:M_.maximum_endo_lead-1,
-    A1 = A1*A;
-    B1 = A*B1;
-    k = find(kstate(:,2) == M_.maximum_endo_lag+2+j );
-    Ax = [A1(kstate(k,1),:); Ax];
-    Bx = [B1(kstate(k,1),:); Bx];
-end
-Ax = [zeros(m+nauxe,nauxe) Ax];
-A0 = A;
-A=Ax; clear Ax A1;
-B0=B;
-B = Bx; clear Bx B1;
+% Dg3 = cat(1,Dg3,zeros(naux+nauxe,size(GAM3,2),param_nbr));
+
+% auxe = zeros(0,1);
+% k0 = kstate(find(kstate(:,2) >= M_.maximum_endo_lag+2),:);;
+% i0 = find(k0(:,2) == M_.maximum_endo_lag+2);
+% for i=M_.maximum_endo_lag+3:M_.maximum_endo_lag+2+M_.maximum_endo_lead,
+%     i1 = find(k0(:,2) == i);
+%     n1 = size(i1,1);
+%     j = zeros(n1,1);
+%     for j1 = 1:n1
+%         j(j1) = find(k0(i0,1)==k0(i1(j1),1));
+%     end
+%     auxe = [auxe; i0(j)];
+%     i0 = i1;
+% end
+% auxe = [(1:size(auxe,1))' auxe(end:-1:1)];
+% n_ir1 = size(auxe,1);
+% nr = m + n_ir1;
+% GAM1 = [[GAM1 zeros(size(GAM1,1),naux)]; zeros(naux+nauxe,m+nauxe)];
+% GAM1(m+1:end,:) = sparse(auxe(:,1),auxe(:,2),ones(n_ir1,1),n_ir1,nr);
+% Dg1 = cat(2,Dg1,zeros(M_.endo_nbr,naux,param_nbr));
+% Dg1 = cat(1,Dg1,zeros(naux+nauxe,m+nauxe,param_nbr));
+
+% Ax = A;
+% A1 = A;
+% Bx = B;
+% B1 = B;
+% for j=1:M_.maximum_endo_lead-1,
+%     A1 = A1*A;
+%     B1 = A*B1;
+%     k = find(kstate(:,2) == M_.maximum_endo_lag+2+j );
+%     Ax = [A1(kstate(k,1),:); Ax];
+%     Bx = [B1(kstate(k,1),:); Bx];
+% end
+% Ax = [zeros(m+nauxe,nauxe) Ax];
+% A0 = A;
+% A=Ax; clear Ax A1;
+% B0=B;
+% B = Bx; clear Bx B1;
 
 m = size(A,1);
 n = size(B,2);
@@ -312,25 +356,35 @@ else % generalized sylvester equation
         d(:,:,j) = Dg2(:,:,j)-elem(:,:,j)*A;
     end
     xx=sylvester3mr(a,b,c,d);
+    H=zeros(m*m+m*(m+1)/2,param_nbr+length(indexo));
+    if nargout>1,
+        dOm = zeros(m,m,param_nbr+length(indexo));
+        dA=zeros(m,m,param_nbr+length(indexo));
+        dB=zeros(m,n,param_nbr);
+    end
     if ~isempty(indexo),
-        dSig = zeros(M_.exo_nbr,M_.exo_nbr);
+        dSig = zeros(M_.exo_nbr,M_.exo_nbr,length(indexo));
         for j=1:length(indexo)
-            dSig(indexo(j),indexo(j)) = 2*sqrt(M_.Sigma_e(indexo(j),indexo(j)));
-            y = B*dSig*B';
-            y = y(nauxe+1:end,nauxe+1:end);
-            H(:,j) = [zeros((m-nauxe)^2,1); dyn_vech(y)];
+            dSig(indexo(j),indexo(j),j) = 2*sqrt(M_.Sigma_e(indexo(j),indexo(j)));
+            y = B*dSig(:,:,j)*B';
+%             y = y(nauxe+1:end,nauxe+1:end);
+%             H(:,j) = [zeros((m-nauxe)^2,1); dyn_vech(y)];
+            H(:,j) = [zeros(m^2,1); dyn_vech(y)];
             if nargout>1,
                 dOm(:,:,j) = y;
             end
-            dSig(indexo(j),indexo(j)) = 0;
+%             dSig(indexo(j),indexo(j)) = 0;
         end
     end
     for j=1:param_nbr,
         x = xx(:,:,j);
         y = inva * (Dg3(:,:,j)-(elem(:,:,j)-GAM1*x)*B);
+        if nargout>1,
+            dB(:,:,j) = y;
+        end
         y = y*M_.Sigma_e*B'+B*M_.Sigma_e*y';
-        x = x(nauxe+1:end,nauxe+1:end);
-        y = y(nauxe+1:end,nauxe+1:end);
+%         x = x(nauxe+1:end,nauxe+1:end);
+%         y = y(nauxe+1:end,nauxe+1:end);
         if nargout>1,
             dA(:,:,j+length(indexo)) = x;
             dOm(:,:,j+length(indexo)) = y;
@@ -353,7 +407,120 @@ else % generalized sylvester equation
 
 end
 
+if nargout > 5,
+    tot_param_nbr = param_nbr + length(indexo);
+    tot_param_nbr_2 = tot_param_nbr*(tot_param_nbr+1)/2;
+    d = zeros(m,m,param_nbr_2);
+    d2A = zeros(m,m,tot_param_nbr,tot_param_nbr);
+    d2Om = zeros(m,m,tot_param_nbr,tot_param_nbr);
+    d2B = zeros(m,n,tot_param_nbr,tot_param_nbr);
+    cc=triu(ones(param_nbr,param_nbr));
+    [i2,j2]=find(cc);
+    cc = blkdiag( zeros(length(indexo),length(indexo)), cc);
+    [ipar2]=find(cc);
+    ctot=triu(ones(tot_param_nbr,tot_param_nbr));
+    ctot(1:length(indexo),1:length(indexo))=eye(length(indexo));
+    [itot2, jtot2]=find(ctot);
+    jcount=0;
+    for j=1:param_nbr,
+        for i=j:param_nbr,
+        elem1 = (D2g0(:,:,j,i)-D2g1(:,:,j,i)*A);
+        elem1 = D2g2(:,:,j,i)-elem1*A;
+        elemj0 = Dg0(:,:,j)-Dg1(:,:,j)*A;
+        elemi0 = Dg0(:,:,i)-Dg1(:,:,i)*A;
+        elem2 = -elemj0*xx(:,:,i)-elemi0*xx(:,:,j);
+        elem2 = elem2 + ( Dg1(:,:,j)*xx(:,:,i) + Dg1(:,:,i)*xx(:,:,j) )*A;
+        elem2 = elem2 + GAM1*( xx(:,:,i)*xx(:,:,j) + xx(:,:,j)*xx(:,:,i));
+        jcount=jcount+1;
+        d(:,:,jcount) = elem1+elem2;
+        end
+    end
+    xx2=sylvester3mr(a,b,c,d);
+    jcount = 0;
+    for j=1:param_nbr,
+        for i=j:param_nbr,
+        jcount=jcount+1;
+        x = xx2(:,:,jcount);
+        elem1 = (D2g0(:,:,j,i)-D2g1(:,:,j,i)*A);
+        elem1 = elem1 -( Dg1(:,:,j)*xx(:,:,i) + Dg1(:,:,i)*xx(:,:,j) );
+        elemj0 = Dg0(:,:,j)-Dg1(:,:,j)*A-GAM1*xx(:,:,j);
+        elemi0 = Dg0(:,:,i)-Dg1(:,:,i)*A-GAM1*xx(:,:,i);
+        elem0 = elemj0*dB(:,:,i)+elemi0*dB(:,:,j);
+        y = inva * (D2g3(:,:,j,i)-elem0-(elem1-GAM1*x)*B);
+        d2B(:,:,j+length(indexo),i+length(indexo)) = y;
+        d2B(:,:,i+length(indexo),j+length(indexo)) = y;
+        y = y*M_.Sigma_e*B'+B*M_.Sigma_e*y'+ ...
+            dB(:,:,j)*M_.Sigma_e*dB(:,:,i)'+dB(:,:,i)*M_.Sigma_e*dB(:,:,j)';
+%         x = x(nauxe+1:end,nauxe+1:end);
+%         y = y(nauxe+1:end,nauxe+1:end);
+        d2A(:,:,j+length(indexo),i+length(indexo)) = x;
+        d2A(:,:,i+length(indexo),j+length(indexo)) = x;
+        d2Om(:,:,j+length(indexo),i+length(indexo)) = y;
+        d2Om(:,:,i+length(indexo),j+length(indexo)) = y;
+        end
+    end    
+    if ~isempty(indexo),
+        d2Sig = zeros(M_.exo_nbr,M_.exo_nbr,length(indexo));
+        for j=1:length(indexo)
+            d2Sig(indexo(j),indexo(j),j) = 2;
+            y = B*d2Sig(:,:,j)*B';
+            d2Om(:,:,j,j) = y;
+%             y = y(nauxe+1:end,nauxe+1:end);
+%             H(:,j) = [zeros((m-nauxe)^2,1); dyn_vech(y)];
+%             H(:,j) = [zeros(m^2,1); dyn_vech(y)];
+%             dOm(:,:,j) = y;
+            for i=1:param_nbr, 
+                y = dB(:,:,i)*dSig(:,:,j)*B'+B*dSig(:,:,j)*dB(:,:,i)';
+                d2Om(:,:,j,i+length(indexo)) = y;
+                d2Om(:,:,i+length(indexo),j) = y;
+            end
+        end
+    end
+end
 
 return
 
+function g22 = get_2nd_deriv(gpp,m,n,i,j),
+
+g22=zeros(m,n);
+is=find(gpp(:,3)==i);
+is=is(find(gpp(is,4)==j));
+
+if ~isempty(is),
+    g22(gpp(is,1),gpp(is,2))=gpp(is,5);
+end
+return
+
+function g22 = get_all_2nd_derivs(gpp,m,n,npar),
+
+g22=zeros(m,n,npar,npar);
+% c=ones(npar,npar);
+% c=triu(c);
+% ic=find(c);
+
+for is=1:length(gpp),
+%     d=zeros(npar,npar);
+%     d(gpp(is,3),gpp(is,4))=1;
+%     indx = find(ic==find(d));
+    g22(gpp(is,1),gpp(is,2),gpp(is,3),gpp(is,4))=gpp(is,5);
+    g22(gpp(is,1),gpp(is,2),gpp(is,4),gpp(is,3))=gpp(is,5);
+end
+
+return
+
+function r22 = get_all_resid_2nd_derivs(rpp,m,npar),
+
+r22=zeros(m,npar,npar);
+% c=ones(npar,npar);
+% c=triu(c);
+% ic=find(c);
+
+for is=1:length(rpp),
+%     d=zeros(npar,npar);
+%     d(rpp(is,2),rpp(is,3))=1;
+%     indx = find(ic==find(d));
+    r22(rpp(is,1),rpp(is,2),rpp(is,3))=rpp(is,4);
+    r22(rpp(is,1),rpp(is,3),rpp(is,2))=rpp(is,4);
+end
 
+return
\ No newline at end of file
diff --git a/matlab/get_Hessian.m b/matlab/get_Hessian.m
new file mode 100644
index 0000000000000000000000000000000000000000..a42b3c81250b89ad85163dcfa46eff7f40094fe3
--- /dev/null
+++ b/matlab/get_Hessian.m
@@ -0,0 +1,258 @@
+function [Hess] = get_Hessian(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,D2T,D2Yss,D2Om,D2H,D2P,start,mf,kalman_tol,riccati_tol)
+% function [Hess] = get_Hessian(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,D2T,D2Yss,D2Om,D2H,D2P,start,mf,kalman_tol,riccati_tol)
+%
+% computes the hessian matrix of the log-likelihood function of
+% a state space model (notation as in kalman_filter.m in DYNARE
+% Thanks to  Nikolai Iskrev
+%
+% NOTE: the derivative matrices (DT,DR ...) are 3-dim. arrays with last
+% dimension equal to the number of structural parameters
+% NOTE: the derivative matrices (D2T,D2Om ...) are 4-dim. arrays with last
+% two dimensions equal to the number of structural parameters
+
+% Copyright (C) 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/>.
+
+
+    k = size(DT,3);                                 % number of structural parameters
+    smpl = size(Y,2);                               % Sample size.
+    pp   = size(Y,1);                               % Maximum number of observed variables.
+    mm   = size(T,2);                               % Number of state variables.
+    a    = zeros(mm,1);                             % State vector.
+    Om   = R*Q*transpose(R);                        % Variance of R times the vector of structural innovations.
+    t    = 0;                                       % Initialization of the time index.
+    oldK = 0;
+    notsteady   = 1;                                % Steady state flag.
+    F_singular  = 1;
+
+    Hess  = zeros(k,k);                             % Initialization of the Hessian
+    Da    = zeros(mm,k);                             % State vector.
+    Dv = zeros(length(mf),k);
+    D2a    = zeros(mm,k,k);                             % State vector.
+    D2v = zeros(length(mf),k,k);
+
+    C = zeros(length(mf),mm);
+    for ii=1:length(mf); C(ii,mf(ii))=1;end         % SELECTION MATRIX IN MEASUREMENT EQ. (FOR WHEN IT IS NOT CONSTANT)
+    dC = zeros(length(mf),mm,k);
+    d2C = zeros(length(mf),mm,k,k);
+    
+    s   = zeros(pp,1);                      % CONSTANT TERM IN MEASUREMENT EQ. (FOR WHEN IT IS NOT CONSTANT)
+    ds  = zeros(pp,1,k);
+    d2s = zeros(pp,1,k,k);
+    
+%     for ii = 1:k
+%         DOm = DR(:,:,ii)*Q*transpose(R) + R*DQ(:,:,ii)*transpose(R) + R*Q*transpose(DR(:,:,ii)); 
+%     end
+    
+    while notsteady & t<smpl
+        t  = t+1;
+        v  = Y(:,t)-a(mf);
+        F  = P(mf,mf) + H;
+        if rcond(F) < kalman_tol
+            if ~all(abs(F(:))<kalman_tol)
+                return
+            else
+                a = T*a;
+                P = T*P*transpose(T)+Om;
+            end
+        else
+            F_singular = 0;
+            iF     = inv(F);
+            K      = P(:,mf)*iF;
+
+            [DK,DF,DP1] = computeDKalman(T,DT,DOm,P,DP,DH,mf,iF,K);
+            [D2K,D2F,D2P1] = computeD2Kalman(T,DT,D2T,D2Om,P,DP,D2P,DH,mf,iF,K,DK);
+            tmp = (a+K*v);
+
+    for ii = 1:k
+        Dv(:,ii)   = -Da(mf,ii) - DYss(mf,ii);
+      %  dai = da(:,:,ii);
+        dKi  = DK(:,:,ii);
+        diFi = -iF*DF(:,:,ii)*iF;
+        dtmpi = Da(:,ii)+dKi*v+K*Dv(:,ii);
+            
+        
+        for jj = 1:ii
+            dFj    = DF(:,:,jj);
+            diFj   = -iF*DF(:,:,jj)*iF;
+            dKj  = DK(:,:,jj);
+            d2Kij  = D2K(:,:,jj,ii);
+            d2Fij  = D2F(:,:,jj,ii);
+            d2iFij = -diFi*dFj*iF -iF*d2Fij*iF -iF*dFj*diFi;
+            dtmpj = Da(:,jj)+dKj*v+K*Dv(:,jj);
+
+            d2vij  = -D2Yss(mf,jj,ii)  - D2a(mf,jj,ii); 
+            d2tmpij = D2a(:,jj,ii) + d2Kij*v + dKj*Dv(:,ii) + dKi*Dv(:,jj) + K*d2vij;
+            D2a(:,jj,ii) = D2T(:,:,jj,ii)*tmp + DT(:,:,jj)*dtmpi + DT(:,:,ii)*dtmpj + T*d2tmpij;            
+
+            Hesst(ii,jj) = getHesst_ij(v,Dv(:,ii),Dv(:,jj),d2vij,iF,diFi,diFj,d2iFij,dFj,d2Fij);
+        end
+        Da(:,ii)   = DT(:,:,ii)*tmp + T*dtmpi;
+    end
+%                     vecDPmf = reshape(DP(mf,mf,:),[],k);
+%                     iPmf = inv(P(mf,mf));
+                    if t>=start
+                        Hess = Hess + Hesst;
+                    end
+            a      = T*(a+K*v);                   
+            P      = T*(P-K*P(mf,:))*transpose(T)+Om;
+            DP     = DP1;
+            D2P     = D2P1;
+        end
+        notsteady = max(max(abs(K-oldK))) > riccati_tol;
+        oldK = K;
+    end
+
+    if F_singular
+        error('The variance of the forecast error remains singular until the end of the sample')
+    end
+
+    
+    if t < smpl
+        t0 = t+1;
+        while t < smpl
+            t = t+1;
+            v = Y(:,t)-a(mf);
+            tmp = (a+K*v);
+            for ii = 1:k,
+                Dv(:,ii)   = -Da(mf,ii)-DYss(mf,ii);
+                dKi  = DK(:,:,ii);
+                diFi = -iF*DF(:,:,ii)*iF;
+                dtmpi = Da(:,ii)+dKi*v+K*Dv(:,ii);
+                
+                for jj = 1:ii,
+                    dFj    = DF(:,:,jj);
+                    diFj   = -iF*DF(:,:,jj)*iF;
+                    dKj  = DK(:,:,jj);
+                    d2Kij  = D2K(:,:,jj,ii);
+                    d2Fij  = D2F(:,:,jj,ii);
+                    d2iFij = -diFi*dFj*iF -iF*d2Fij*iF -iF*dFj*diFi;
+                    dtmpj = Da(:,jj)+dKj*v+K*Dv(:,jj);
+                    
+                    d2vij  = -D2Yss(mf,jj,ii)  - D2a(mf,jj,ii);
+                    d2tmpij = D2a(:,jj,ii) + d2Kij*v + dKj*Dv(:,ii) + dKi*Dv(:,jj) + K*d2vij;
+                    D2a(:,jj,ii) = D2T(:,:,jj,ii)*tmp + DT(:,:,jj)*dtmpi + DT(:,:,ii)*dtmpj + T*d2tmpij;            
+                    
+                    Hesst(ii,jj) = getHesst_ij(v,Dv(:,ii),Dv(:,jj),d2vij,iF,diFi,diFj,d2iFij,dFj,d2Fij);
+                end
+                Da(:,ii)   = DT(:,:,ii)*tmp + T*dtmpi;
+            end
+            if t>=start
+                Hess = Hess + Hesst;
+            end
+            a = T*(a+K*v);
+        end
+%         Hess = Hess + .5*(smpl+t0-1)*(vecDPmf' * kron(iPmf,iPmf) * vecDPmf);
+        %         for ii = 1:k;
+        %             for jj = 1:ii
+        %              H(ii,jj) = trace(iPmf*(.5*DP(mf,mf,ii)*iPmf*DP(mf,mf,jj) + Dv(:,ii)*Dv(:,jj)'));
+        %             end
+        %         end
+    end
+    
+Hess = Hess + tril(Hess,-1)';
+
+Hess = -Hess/2;  
+% end of main function    
+
+function Hesst_ij = getHesst_ij(e,dei,dej,d2eij,iS,diSi,diSj,d2iSij,dSj,d2Sij);
+% computes (i,j) term in the Hessian
+
+Hesst_ij = trace(diSi*dSj + iS*d2Sij) + e'*d2iSij*e + 2*(dei'*diSj*e + dei'*iS*dej + e'*diSi*dej + e'*iS*d2eij);
+
+% end of getHesst_ij
+
+function [DK,DF,DP1] = computeDKalman(T,DT,DOm,P,DP,DH,mf,iF,K)
+
+            k      = size(DT,3);
+            tmp    = P-K*P(mf,:);
+
+for ii = 1:k
+    DF(:,:,ii)  = DP(mf,mf,ii) + DH(:,:,ii); 
+    DiF(:,:,ii) = -iF*DF(:,:,ii)*iF;
+    DK(:,:,ii)  = DP(:,mf,ii)*iF + P(:,mf)*DiF(:,:,ii);
+    Dtmp        = DP(:,:,ii) - DK(:,:,ii)*P(mf,:) - K*DP(mf,:,ii);
+    DP1(:,:,ii) = DT(:,:,ii)*tmp*T' + T*Dtmp*T' + T*tmp*DT(:,:,ii)' + DOm(:,:,ii);
+end
+
+% end of computeDKalman
+
+function [d2K,d2S,d2P1] = computeD2Kalman(A,dA,d2A,d2Om,P0,dP0,d2P0,DH,mf,iF,K0,dK0);
+% computes the second derivatives of the Kalman matrices
+% note: A=T in main func.
+        
+            k      = size(dA,3);
+            tmp    = P0-K0*P0(mf,:);
+[ns,no] = size(K0);
+
+% CPC = C*P0*C'; CPC = .5*(CPC+CPC');iF = inv(CPC);
+% APC = A*P0*C';
+% APA = A*P0*A';
+
+
+d2K  = zeros(ns,no,k,k);
+d2S  = zeros(no,no,k,k);
+d2P1 = zeros(ns,ns,k,k);
+
+for ii = 1:k
+    dAi = dA(:,:,ii);
+    dFi = dP0(mf,mf,ii);
+    d2Omi = d2Om(:,:,ii);
+    diFi = -iF*dFi*iF;
+    dKi = dK0(:,:,ii);
+    for jj = 1:k
+        dAj = dA(:,:,jj);
+        dFj = dP0(mf,mf,jj);
+        d2Omj = d2Om(:,:,jj);
+        dFj = dP0(mf,mf,jj);
+        diFj = -iF*dFj*iF;
+        dKj = dK0(:,:,jj);
+
+        d2Aij = d2A(:,:,jj,ii);
+        d2Pij = d2P0(:,:,jj,ii);
+        d2Omij = d2Om(:,:,jj,ii);
+       
+    % second order
+    
+    d2Fij = d2Pij(mf,mf) ;
+    
+%     d2APC = d2Aij*P0*C' + A*d2Pij*C' + A*P0*d2Cij' + dAi*dPj*C' + dAj*dPi*C' + A*dPj*dCi' + A*dPi*dCj' + dAi*P0*dCj' + dAj*P0*dCi';
+    d2APC = d2Pij(:,mf);
+    
+    d2iF = -diFi*dFj*iF -iF*d2Fij*iF -iF*dFj*diFi;
+    
+    d2Kij= d2Pij(:,mf)*iF + P0(:,mf)*d2iF + dP0(:,mf,jj)*diFi + dP0(:,mf,ii)*diFj;
+        
+    d2KCP = d2Kij*P0(mf,:) + K0*d2Pij(mf,:) + dKi*dP0(mf,:,jj) + dKj*dP0(mf,:,ii) ;
+    
+    dtmpi        = dP0(:,:,ii) - dK0(:,:,ii)*P0(mf,:) - K0*dP0(mf,:,ii);
+    dtmpj        = dP0(:,:,jj) - dK0(:,:,jj)*P0(mf,:) - K0*dP0(mf,:,jj);
+    d2tmp = d2Pij - d2KCP;
+
+    d2AtmpA = d2Aij*tmp*A' + A*d2tmp*A' + A*tmp*d2Aij' + dAi*dtmpj*A' + dAj*dtmpi*A' + A*dtmpj*dAi' + A*dtmpi*dAj' + dAi*tmp*dAj' + dAj*tmp*dAi';
+
+    d2K(:,:,ii,jj)  = d2Kij; %#ok<NASGU>
+    d2P1(:,:,ii,jj) = d2AtmpA  + d2Omij;  %#ok<*NASGU>
+    d2S(:,:,ii,jj)  = d2Fij;
+%     d2iS(:,:,ii,jj) = d2iF;
+    end
+end
+
+% end of computeD2Kalman
+
+
+            
\ No newline at end of file
diff --git a/matlab/lpdfgam.m b/matlab/lpdfgam.m
index 0a8505960f3b2a453c4276e940fb387f7176a612..1246f2a9e9586bda1933ba0c52660a31e7dde15f 100644
--- a/matlab/lpdfgam.m
+++ b/matlab/lpdfgam.m
@@ -1,4 +1,4 @@
-function  ldens = lpdfgam(x,a,b);
+function  [ldens,Dldens,D2ldens] = lpdfgam(x,a,b);
 % Evaluates the logged GAMMA PDF at x.
 %
 % INPUTS     
@@ -37,4 +37,22 @@ if length(a)==1
     ldens(idx) = -gammaln(a) - a*log(b) + (a-1)*log(x(idx)) - x(idx)/b ;
 else
     ldens(idx) = -gammaln(a(idx)) - a(idx).*log(b(idx)) + (a(idx)-1).*log(x(idx)) - x(idx)./b(idx) ;
+end
+
+
+
+if nargout >1 
+    if length(a)==1
+        Dldens(idx) = (a-1)./(x(idx)) - ones(length(idx),1)/b ;
+    else
+        Dldens(idx) = (a(idx)-1)./(x(idx)) - ones(length(idx),1)./b(idx) ;
+    end
+end
+
+if nargout == 3
+    if length(a)==1
+        D2ldens(idx) = -(a-1)./(x(idx)).^2;
+    else
+        D2ldens(idx) = -(a(idx)-1)./(x(idx)).^2;
+    end
 end
\ No newline at end of file
diff --git a/matlab/lpdfgbeta.m b/matlab/lpdfgbeta.m
index 33ff3df48ae142b639c6053f9917151f682139f0..27a3c7fde7a17a4f77e2692ccac9e5a1b64ad4cc 100644
--- a/matlab/lpdfgbeta.m
+++ b/matlab/lpdfgbeta.m
@@ -1,4 +1,4 @@
-function ldens = lpdfgbeta(x,a,b,aa,bb);
+function [ldens,Dldens,D2ldens] = lpdfgbeta(x,a,b,aa,bb);
 % Evaluates the logged BETA PDF at x. 
 %
 % INPUTS 
@@ -38,4 +38,22 @@ if length(a)==1
     ldens(idx) = -betaln(a,b) + (a-1)*log(x(idx)-aa) + (b-1)*log(bb-x(idx)) - (a+b-1)*log(bb-aa) ;
 else
     ldens(idx) = -betaln(a(idx),b(idx)) + (a(idx)-1).*log(x(idx)-aa(idx)) + (b(idx)-1).*log(bb(idx)-x(idx)) - (a(idx)+b(idx)-1).*log(bb(idx)-aa(idx));
+end
+
+
+if nargout >1 
+    if length(a)==1
+        Dldens(idx) = (a-1)./(x(idx)-aa) - (b-1)./(bb-x(idx)) ;
+    else
+        Dldens(idx) = (a(idx)-1)./(x(idx)-aa(idx)) - (b(idx)-1)./(bb(idx)-x(idx));
+    end
+end
+
+
+if nargout == 3 
+    if length(a)==1
+        D2ldens(idx) = -(a-1)./(x(idx)-aa).^2 - (b-1)./(bb-x(idx)).^2 ;
+    else
+        D2ldens(idx) = -(a(idx)-1)./(x(idx)-aa(idx)).^2 - (b(idx)-1)./(bb(idx)-x(idx)).^2;
+    end
 end
\ No newline at end of file
diff --git a/matlab/lpdfig1.m b/matlab/lpdfig1.m
index 9c17277a84231977a6ab49a0f01860aa652abc2a..c0bb8ff770ff0b8994d5a444a197ac71405cdeab 100644
--- a/matlab/lpdfig1.m
+++ b/matlab/lpdfig1.m
@@ -1,4 +1,4 @@
-function ldens = lpdfig1(x,s,nu)
+function [ldens,Dldens,D2ldens] = lpdfig1(x,s,nu)
 % Evaluates the logged INVERSE-GAMMA-1 PDF at x.
 %
 % X ~ IG1(s,nu) if X = sqrt(Y) where Y ~ IG2(s,nu) and Y = inv(Z) with Z ~ G(nu/2,2/s) (Gamma distribution) 
@@ -41,4 +41,21 @@ if length(s)==1
     ldens(idx) = log(2) - gammaln(.5*nu) - .5*nu*(log(2)-log(s)) - (nu+1)*log(x(idx)) - .5*s./(x(idx).*x(idx)) ;
 else
     ldens(idx) = log(2) - gammaln(.5*nu(idx)) - .5*nu(idx).*(log(2)-log(s(idx))) - (nu(idx)+1).*log(x(idx)) - .5*s(idx)./(x(idx).*x(idx)) ;
-end
\ No newline at end of file
+end
+
+if nargout >1 
+    if length(s)==1
+        Dldens(idx) = - (nu+1)./(x(idx)) + s./(x(idx).^3) ;
+    else
+        Dldens(idx) = - (nu(idx)+1)./(x(idx)) + s(idx)./(x(idx).^3) ;
+    end
+end
+
+if nargout == 3 
+    if length(s)==1
+        D2ldens(idx) =  (nu+1)./(x(idx).^2) - 3*s(idx)./(x(idx).^4) ;
+    else
+        D2ldens(idx) =  (nu(idx)+1)./(x(idx).^2) - 3*s(idx)./(x(idx).^4) ;
+    end
+end
+
diff --git a/matlab/lpdfig2.m b/matlab/lpdfig2.m
index 7f95581811503bc5e641e7e8dfe44ca65d2367df..63e458db000ac9759a0c5d77f852af751e16b6ea 100644
--- a/matlab/lpdfig2.m
+++ b/matlab/lpdfig2.m
@@ -1,4 +1,4 @@
-function ldens = lpdfig2(x,s,nu)
+function [ldens,Dldens,D2ldens] = lpdfig2(x,s,nu)
 % Evaluates the logged INVERSE-GAMMA-2 PDF at x.
 %
 % X ~ IG2(s,nu) if X = inv(Z) where Z ~ G(nu/2,2/s) (Gamma distribution) 
@@ -41,4 +41,20 @@ if length(s)==1
     ldens(idx) = -gammaln(.5*nu) - (.5*nu)*(log(2)-log(s)) - .5*(nu+2)*log(x(idx)) -.5*s./x(idx);
 else
     ldens(idx) = -gammaln(.5*nu(idx)) - (.5*nu(idx)).*(log(2)-log(s(idx))) - .5*(nu(idx)+2).*log(x(idx)) -.5*s(idx)./x(idx);
+end
+
+if nargout >1 
+    if length(s)==1
+        Dldens(idx) = - .5*(nu+2)./(x(idx)) + .5*s./x(idx).^2;
+    else
+        Dldens(idx) = - .5*(nu(idx)+2)./(x(idx)) + .5*s(idx)./x(idx).^2;
+    end
+end
+
+if nargout == 3 
+    if length(s)==1
+        D2ldens(idx) = .5*(nu+2)./(x(idx)).^2 - s./x(idx).^3;
+    else
+        D2ldens(idx) = .5*(nu(idx)+2)./(x(idx)).^2 - s(idx)./x(idx).^3;
+    end
 end
\ No newline at end of file
diff --git a/matlab/lpdfnorm.m b/matlab/lpdfnorm.m
index 0e232884a4e1a6d8481d0e6946e0fa6cb11ab3d7..9afdd2b6df2076e4210793049e29164f948c4482 100644
--- a/matlab/lpdfnorm.m
+++ b/matlab/lpdfnorm.m
@@ -1,4 +1,4 @@
-function  ldens = lpdfnorm(x,a,b)
+function  [ldens,Dldens,D2ldens] = lpdfnorm(x,a,b)
 % Evaluates the logged UNIVARIATE GAUSSIAN PDF at x.
 %
 % INPUTS 
@@ -32,4 +32,12 @@ function  ldens = lpdfnorm(x,a,b)
 
 if nargin<3, b=1; end
 if nargin<2, a=0; end
-ldens = -log(b) -.5*log(2*pi) - .5*((x-a)./b).*((x-a)./b) ;
\ No newline at end of file
+ldens = -log(b) -.5*log(2*pi) - .5*((x-a)./b).*((x-a)./b) ;
+
+if nargout >1 
+    Dldens =  - (1/b)*((x-a)/b) ;
+end
+
+if nargout == 3 
+    D2ldens =  - (1/b)^2 ;
+end
\ No newline at end of file
diff --git a/matlab/mode_check.m b/matlab/mode_check.m
index 8e14a3e204c8e7d3cb7710b0289ec149185b03f7..36d516c89b7c60c87bc97a22eae9321af58f7526 100644
--- a/matlab/mode_check.m
+++ b/matlab/mode_check.m
@@ -1,4 +1,4 @@
-function mode_check(func,x,fval,hessian,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults)
+function mode_check(fun,x,hessian,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults)
 
 % function mode_check(x,fval,hessian,gend,data,lb,ub)
 % Checks the maximum likelihood mode
@@ -91,7 +91,7 @@ for plt = 1:nbplt,
         end
         for i=1:length(z)
             xx(kk) = z(i);
-            [fval, exit_flag] = feval(fun,x,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
+            [fval, exit_flag] = feval(fun,xx,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
             if exit_flag
                 y(i,1) = fval;
             else
diff --git a/matlab/newrat.m b/matlab/newrat.m
index bd59de7b9a384e66168b2c094817e87724dbf7a0..619dd90a008d53eaff34bc5db62607d459e25964 100644
--- a/matlab/newrat.m
+++ b/matlab/newrat.m
@@ -46,7 +46,7 @@ icount=0;
 nx=length(x);
 xparam1=x;
 %ftol0=1.e-6;
-htol_base = max(1.e-5, ftol0);
+htol_base = max(1.e-7, ftol0);
 flagit=0;  % mode of computation of hessian in each iteration
 ftol=ftol0;
 gtol=1.e-3;
diff --git a/matlab/priordens.m b/matlab/priordens.m
index c1b64fe1ea7ed5decffb0ee56b401a7ca6406400..62cbda5d7f634e6b13a31681a9a9f02be72f6ff7 100644
--- a/matlab/priordens.m
+++ b/matlab/priordens.m
@@ -1,4 +1,4 @@
-function logged_prior_density = priordens(x, pshape, p6, p7, p3, p4,initialization)
+function [logged_prior_density, dlprior, d2lprior] = priordens(x, pshape, p6, p7, p3, p4,initialization)
 % Computes a prior density for the structural parameters of DSGE models
 %
 % INPUTS 
@@ -81,6 +81,11 @@ if tt1
     if isinf(logged_prior_density)
         return
     end
+    if nargout == 2,
+        [tmp, dlprior(id1)]=lpdfgbeta(x(id1),p6(id1),p7(id1),p3(id1),p4(id1));
+    elseif nargout == 3
+        [tmp, dlprior(id1), d2lprior(id1)]=lpdfgbeta(x(id1),p6(id1),p7(id1),p3(id1),p4(id1));
+    end
 end
 
 if tt2
@@ -88,10 +93,20 @@ if tt2
     if isinf(logged_prior_density)
         return
     end
+    if nargout == 2,
+        [tmp, dlprior(id2)]=lpdfgam(x(id2)-p3(id2),p6(id2),p7(id2));
+    elseif nargout == 3
+        [tmp, dlprior(id2), d2lprior(id2)]=lpdfgam(x(id2)-p3(id2),p6(id2),p7(id2));
+    end
 end
 
 if tt3
     logged_prior_density = logged_prior_density + sum(lpdfnorm(x(id3),p6(id3),p7(id3))) ;
+    if nargout == 2,
+        [tmp, dlprior(id3)]=lpdfnorm(x(id3),p6(id3),p7(id3));
+    elseif nargout == 3
+        [tmp, dlprior(id3), d2lprior(id3)]=lpdfnorm(x(id3),p6(id3),p7(id3));
+    end
 end
 
 if tt4
@@ -99,6 +114,11 @@ if tt4
     if isinf(logged_prior_density)
         return
     end
+    if nargout == 2,
+        [tmp, dlprior(id4)]=lpdfig1(x(id4)-p3(id4),p6(id4),p7(id4));
+    elseif nargout == 3
+        [tmp, dlprior(id4), d2lprior(id4)]=lpdfig1(x(id4)-p3(id4),p6(id4),p7(id4));
+    end
 end
 
 if tt5
@@ -107,6 +127,12 @@ if tt5
         return
     end
     logged_prior_density = logged_prior_density + sum(log(1./(p4(id5)-p3(id5)))) ;
+    if nargout >1,
+        dlprior(id5)=zeros(length(id5),1);
+    end
+    if nargout == 3
+        d2lprior(id5)=zeros(length(id5),1);
+    end
 end
 
 if tt6
@@ -114,4 +140,13 @@ if tt6
     if isinf(logged_prior_density)
         return
     end
+    if nargout == 2,
+        [tmp, dlprior(id6)]=lpdfig2(x(id6)-p3(id6),p6(id6),p7(id6));
+    elseif nargout == 3
+        [tmp, dlprior(id6), d2lprior(id6)]=lpdfig2(x(id6)-p3(id6),p6(id6),p7(id6));
+    end
+end
+
+if nargout==3,
+    d2lprior = diag(d2lprior);
 end
\ No newline at end of file
diff --git a/matlab/score.m b/matlab/score.m
index 9a41214afc259377575ccfc0bb1edf36cec26edf..dba04b84ce53796f428528e1e080ec5d76338a24 100644
--- a/matlab/score.m
+++ b/matlab/score.m
@@ -98,7 +98,7 @@ function [DLIK] = score(T,R,Q,H,P,Y,DT,DYss,DOm,DH,DP,start,mf,kalman_tol,riccat
             a = T*(a+K*v);
         end
         for ii = 1:k
-            DLIK(ii,1)  = DLIK(ii,1) + (smpl-t0+1)*trace( iF*DF(:,:,ii) );
+%             DLIK(ii,1)  = DLIK(ii,1) + (smpl-t0+1)*trace( iF*DF(:,:,ii) );
         end
         
     end