diff --git a/matlab/initial_estimation_checks.m b/matlab/initial_estimation_checks.m
index 0d6463f19..53da3c263 100644
--- a/matlab/initial_estimation_checks.m
+++ b/matlab/initial_estimation_checks.m
@@ -46,11 +46,11 @@ else
     info = d;
 end
 
-if DynareOptions.mode_compute==5
-    if ~strcmp(func2str(objective_function),'dsge_likelihood')
-        error('Options mode_compute=5 is not compatible with non linear filters or Dsge-VAR models!')
-    end
-end
+% if DynareOptions.mode_compute==5
+%     if ~strcmp(func2str(objective_function),'dsge_likelihood')
+%         error('Options mode_compute=5 is not compatible with non linear filters or Dsge-VAR models!')
+%     end
+% end
 
 if info(1) > 0
     disp('Error in computing likelihood for initial parameter values')
diff --git a/matlab/mr_hessian.m b/matlab/mr_hessian.m
index 99149969f..46a8db8de 100644
--- a/matlab/mr_hessian.m
+++ b/matlab/mr_hessian.m
@@ -11,7 +11,7 @@ function [hessian_mat, gg, htol1, ihh, hh_mat0, hh1] = mr_hessian(init,x,func,hf
 %    of the log-likelihood to compute outer product gradient
 %  x = parameter values
 %  hflag = 0, Hessian computed with outer product gradient, one point
-%  increments for partial derivatives in  gradients
+%  increments for partial derivatives in gradients
 %  hflag = 1, 'mixed' Hessian: diagonal elements computed with numerical second order derivatives
 %             with correlation structure as from outer product gradient;
 %             two point evaluation of derivatives for partial derivatives
@@ -55,6 +55,12 @@ end
 h2=BayesInfo.ub-BayesInfo.lb;
 hmax=BayesInfo.ub-x;
 hmax=min(hmax,x-BayesInfo.lb);
+if isempty(ff0),
+    outer_product_gradient=0;
+else
+    outer_product_gradient=1;
+end
+
 
 h1 = min(h1,0.5.*hmax);
 
@@ -64,9 +70,11 @@ end
 xh1=x;
 f1=zeros(size(f0,1),n);
 f_1=f1;
-ff1=zeros(size(ff0));
-ff_1=ff1;
-ggh=zeros(size(ff0,1),n);
+if outer_product_gradient
+    ff1=zeros(size(ff0));
+    ff_1=ff1;
+    ggh=zeros(size(ff0,1),n);
+end
 
 i=0;
 while i<n
@@ -124,20 +132,24 @@ while i<n
         end
     end
     f1(:,i)=fx;
-    if any(isnan(ffx))
-        ff1=ones(size(ff0)).*fx/length(ff0);
-    else
-        ff1=ffx;
+    if outer_product_gradient,
+        if any(isnan(ffx))
+            ff1=ones(size(ff0)).*fx/length(ff0);
+        else
+            ff1=ffx;
+        end
     end
     xh1(i)=x(i)-h1(i);
     [fx, ffx]=feval(func,xh1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
     f_1(:,i)=fx;
-    if any(isnan(ffx))
-        ff_1=ones(size(ff0)).*fx/length(ff0);
-    else
-        ff_1=ffx;
+    if outer_product_gradient,
+        if any(isnan(ffx))
+            ff_1=ones(size(ff0)).*fx/length(ff0);
+        else
+            ff_1=ffx;
+        end
+        ggh(:,i)=(ff1-ff_1)./(2.*h1(i));
     end
-    ggh(:,i)=(ff1-ff_1)./(2.*h1(i));
     xh1(i)=x(i);
     if hcheck && htol<1
         htol=min(1,max(min(abs(dx))*2,htol*10));
@@ -152,85 +164,93 @@ xh_1=xh1;
 
 gg=(f1'-f_1')./(2.*h1);
 
-if hflag==2
-    gg=(f1'-f_1')./(2.*h1);
-    hessian_mat = zeros(size(f0,1),n*n);
-    for i=1:n
-        if i > 1
-            k=[i:n:n*(i-1)];
-            hessian_mat(:,(i-1)*n+1:(i-1)*n+i-1)=hessian_mat(:,k);
+if outer_product_gradient,
+    if hflag==2
+        gg=(f1'-f_1')./(2.*h1);
+        hessian_mat = zeros(size(f0,1),n*n);
+        for i=1:n
+            if i > 1
+                k=[i:n:n*(i-1)];
+                hessian_mat(:,(i-1)*n+1:(i-1)*n+i-1)=hessian_mat(:,k);
+            end
+            hessian_mat(:,(i-1)*n+i)=(f1(:,i)+f_1(:,i)-2*f0)./(h1(i)*h_1(i));
+            temp=f1+f_1-f0*ones(1,n);
+            for j=i+1:n
+                xh1(i)=x(i)+h1(i);
+                xh1(j)=x(j)+h_1(j);
+                xh_1(i)=x(i)-h1(i);
+                xh_1(j)=x(j)-h_1(j);
+                temp1 = feval(func,xh1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
+                temp2 = feval(func,xh_1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
+                hessian_mat(:,(i-1)*n+j)=-(-temp1 -temp2+temp(:,i)+temp(:,j))./(2*h1(i)*h_1(j));
+                xh1(i)=x(i);
+                xh1(j)=x(j);
+                xh_1(i)=x(i);
+                xh_1(j)=x(j);
+                j=j+1;
+            end
+            i=i+1;
         end
-        hessian_mat(:,(i-1)*n+i)=(f1(:,i)+f_1(:,i)-2*f0)./(h1(i)*h_1(i));
-        temp=f1+f_1-f0*ones(1,n);
-        for j=i+1:n
-            xh1(i)=x(i)+h1(i);
-            xh1(j)=x(j)+h_1(j);
-            xh_1(i)=x(i)-h1(i);
-            xh_1(j)=x(j)-h_1(j);
-            temp1 = feval(func,xh1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
-            temp2 = feval(func,xh_1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
-            hessian_mat(:,(i-1)*n+j)=-(-temp1 -temp2+temp(:,i)+temp(:,j))./(2*h1(i)*h_1(j));
-            xh1(i)=x(i);
-            xh1(j)=x(j);
-            xh_1(i)=x(i);
-            xh_1(j)=x(j);
-            j=j+1;
-        end
-        i=i+1;
-    end
-elseif hflag==1
-    hessian_mat = zeros(size(f0,1),n*n);
-    for i=1:n
-        dum = (f1(:,i)+f_1(:,i)-2*f0)./(h1(i)*h_1(i));
-        if dum>eps
-            hessian_mat(:,(i-1)*n+i)=dum;
-        else
-            hessian_mat(:,(i-1)*n+i)=max(eps, gg(i)^2);
+    elseif hflag==1
+        hessian_mat = zeros(size(f0,1),n*n);
+        for i=1:n
+            dum = (f1(:,i)+f_1(:,i)-2*f0)./(h1(i)*h_1(i));
+            if dum>eps
+                hessian_mat(:,(i-1)*n+i)=dum;
+            else
+                hessian_mat(:,(i-1)*n+i)=max(eps, gg(i)^2);
+            end
         end
     end
-end
-
-gga=ggh.*kron(ones(size(ff1)),2.*h1');  % re-scaled gradient
-hh_mat=gga'*gga;  % rescaled outer product hessian
-hh_mat0=ggh'*ggh;  % outer product hessian
-A=diag(2.*h1);  % rescaling matrix
-% igg=inv(hh_mat);  % inverted rescaled outer product hessian
-ihh=A'*(hh_mat\A);  % inverted outer product hessian
-if hflag>0 && min(eig(reshape(hessian_mat,n,n)))>0
-    hh0 = A*reshape(hessian_mat,n,n)*A';  %rescaled second order derivatives
-    hh = reshape(hessian_mat,n,n);  %rescaled second order derivatives
-    sd0=sqrt(diag(hh0));   %rescaled 'standard errors' using second order derivatives
-    sd=sqrt(diag(hh_mat));  %rescaled 'standard errors' using outer product
-    hh_mat=hh_mat./(sd*sd').*(sd0*sd0');  %rescaled inverse outer product with 'true' std's
-    igg=inv(hh_mat);   % rescaled outer product hessian with 'true' std's
+    
+    gga=ggh.*kron(ones(size(ff1)),2.*h1');  % re-scaled gradient
+    hh_mat=gga'*gga;  % rescaled outer product hessian
+    hh_mat0=ggh'*ggh;  % outer product hessian
+    A=diag(2.*h1);  % rescaling matrix
+    % igg=inv(hh_mat);  % inverted rescaled outer product hessian
     ihh=A'*(hh_mat\A);  % inverted outer product hessian
-    hh_mat0=inv(A)'*hh_mat*inv(A);  % outer product hessian with 'true' std's
-    sd=sqrt(diag(ihh));   %standard errors
-    sdh=sqrt(1./diag(hh));   %diagonal standard errors
-    for j=1:length(sd)
-        sd0(j,1)=min(BayesInfo.p2(j), sd(j));  %prior std
-        sd0(j,1)=10^(0.5*(log10(sd0(j,1))+log10(sdh(j,1))));
+    if hflag>0 && min(eig(reshape(hessian_mat,n,n)))>0
+        hh0 = A*reshape(hessian_mat,n,n)*A';  %rescaled second order derivatives
+        hh = reshape(hessian_mat,n,n);  %rescaled second order derivatives
+        sd0=sqrt(diag(hh0));   %rescaled 'standard errors' using second order derivatives
+        sd=sqrt(diag(hh_mat));  %rescaled 'standard errors' using outer product
+        hh_mat=hh_mat./(sd*sd').*(sd0*sd0');  %rescaled inverse outer product with 'true' std's
+        igg=inv(hh_mat);   % rescaled outer product hessian with 'true' std's
+        ihh=A'*(hh_mat\A);  % inverted outer product hessian
+        hh_mat0=inv(A)'*hh_mat*inv(A);  % outer product hessian with 'true' std's
+        sd=sqrt(diag(ihh));   %standard errors
+        sdh=sqrt(1./diag(hh));   %diagonal standard errors
+        for j=1:length(sd)
+            sd0(j,1)=min(BayesInfo.p2(j), sd(j));  %prior std
+            sd0(j,1)=10^(0.5*(log10(sd0(j,1))+log10(sdh(j,1))));
+        end
+        ihh=ihh./(sd*sd').*(sd0*sd0');  %inverse outer product with modified std's
+        igg=inv(A)'*ihh*inv(A);  % inverted rescaled outer product hessian with modified std's
+        hh_mat=inv(igg);   % outer product rescaled hessian with modified std's
+        hh_mat0=inv(A)'*hh_mat*inv(A);  % outer product hessian with modified std's
+        %     sd0=sqrt(1./diag(hh0));   %rescaled 'standard errors' using second order derivatives
+        %     sd=sqrt(diag(igg));  %rescaled 'standard errors' using outer product
+        %     igg=igg./(sd*sd').*(sd0*sd0');  %rescaled inverse outer product with 'true' std's
+        %     hh_mat=inv(igg);   % rescaled outer product hessian with 'true' std's
+        %     ihh=A'*igg*A;  % inverted outer product hessian
+        %     hh_mat0=inv(A)'*hh_mat*inv(A);  % outer product hessian with 'true' std's
     end
-    ihh=ihh./(sd*sd').*(sd0*sd0');  %inverse outer product with modified std's
-    igg=inv(A)'*ihh*inv(A);  % inverted rescaled outer product hessian with modified std's
-    hh_mat=inv(igg);   % outer product rescaled hessian with modified std's
-    hh_mat0=inv(A)'*hh_mat*inv(A);  % outer product hessian with modified std's
-                                    %     sd0=sqrt(1./diag(hh0));   %rescaled 'standard errors' using second order derivatives
-                                    %     sd=sqrt(diag(igg));  %rescaled 'standard errors' using outer product
-                                    %     igg=igg./(sd*sd').*(sd0*sd0');  %rescaled inverse outer product with 'true' std's
-                                    %     hh_mat=inv(igg);   % rescaled outer product hessian with 'true' std's
-                                    %     ihh=A'*igg*A;  % inverted outer product hessian
-                                    %     hh_mat0=inv(A)'*hh_mat*inv(A);  % outer product hessian with 'true' std's
-end
-if hflag<2
-    hessian_mat=hh_mat0(:);
+    if hflag<2
+        hessian_mat=hh_mat0(:);
+    end
+    
+    if any(isnan(hessian_mat))
+        hh_mat0=eye(length(hh_mat0));
+        ihh=hh_mat0;
+        hessian_mat=hh_mat0(:);
+    end
+    hh1=h1;
+    save hess.mat hessian_mat
+else
+    hessian_mat=[];
+    ihh=[];
+    hh_mat0 = [];
+    hh1 = [];
 end
 
-if any(isnan(hessian_mat))
-    hh_mat0=eye(length(hh_mat0));
-    ihh=hh_mat0;
-    hessian_mat=hh_mat0(:);
-end
-hh1=h1;
 htol1=htol;
-save hess.mat hessian_mat
\ No newline at end of file
diff --git a/matlab/newrat.m b/matlab/newrat.m
index d907493ea..60c53f3b7 100644
--- a/matlab/newrat.m
+++ b/matlab/newrat.m
@@ -61,16 +61,22 @@ fval=fval0;
 
 % initialize mr_gstep and mr_hessian
 mr_hessian(1,x,[],[],[],DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
+outer_product_gradient=1;
 
 if isempty(hh)
     [dum, gg, htol0, igg, hhg, h1]=mr_hessian(0,x,func0,flagit,htol,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
-    hh0 = reshape(dum,nx,nx);
-    hh=hhg;
-    if min(eig(hh0))<0
-        hh0=hhg; %generalized_cholesky(hh0);
-    elseif flagit==2
-        hh=hh0;
-        igg=inv(hh);
+    if isempty(dum),
+        outer_product_gradient=0;
+        igg = 1e-4*eye(nx);
+    else
+        hh0 = reshape(dum,nx,nx);
+        hh=hhg;
+        if min(eig(hh0))<0
+            hh0=hhg; %generalized_cholesky(hh0);
+        elseif flagit==2
+            hh=hh0;
+            igg=inv(hh);
+        end
     end
     if htol0>htol
         htol=htol0;
@@ -192,6 +198,9 @@ while norm(gg)>gtol && check==0 && jit<nit
                 save m1.mat x fval0 nig
             end
             [dum, gg, htol0, igg, hhg, h1]=mr_hessian(0,xparam1,func0,flagit,htol,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
+            if isempty(dum),
+                outer_product_gradient=0;
+            end
             if htol0>htol
                 htol=htol0;
                 disp(' ')
@@ -199,26 +208,33 @@ while norm(gg)>gtol && check==0 && jit<nit
                 disp('Tolerance has to be relaxed')
                 disp(' ')
             end
-            hh0 = reshape(dum,nx,nx);
-            hh=hhg;
-            if flagit==2
-                if min(eig(hh0))<=0
-                    hh0=hhg; %generalized_cholesky(hh0);
-                else
-                    hh=hh0;
-                    igg=inv(hh);
+            if ~outer_product_gradient,
+                H = bfgsi1(H,gg-g(:,icount),xparam1-x(:,icount));
+                hh=inv(H);
+                hhg=hh;
+            else
+                hh0 = reshape(dum,nx,nx);
+                hh=hhg;
+                if flagit==2
+                    if min(eig(hh0))<=0
+                        hh0=hhg; %generalized_cholesky(hh0);
+                    else
+                        hh=hh0;
+                        igg=inv(hh);
+                    end
                 end
+                H = igg;
             end
         end
         disp(['Gradient norm  ',num2str(norm(gg))])
         ee=eig(hh);
         disp(['Minimum Hessian eigenvalue ',num2str(min(ee))])
         disp(['Maximum Hessian eigenvalue ',num2str(max(ee))])
-        if max(eig(hh))<0, disp('Negative definite Hessian! Local maximum!'), pause, end,
+        if max(eig(hh))<0, disp('Negative definite Hessian! Local maximum!'), pause(1), end,
         t=toc;
         disp(['Elapsed time for iteration ',num2str(t),' s.'])
         g(:,icount+1)=gg;
-        H = igg;
+
         save m1.mat x hh g hhg igg fval0 nig H
     end
 end