From 323e6fdf4f7258e05f03279bc3e89b23f16e44dc Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?=
<stephane.adjemian@univ-lemans.fr>
Date: Fri, 9 Mar 2012 16:36:26 +0100
Subject: [PATCH] Added a second gstep option used as a parameter for the
routine computing the hessian matrix.
In some cases, for instance for the non linear filters, it helps to reduce this new gstep parameter
to get a positive definite hessian matrix. options_.gstep is now a 2*1 vector. The first element is
the old gstep parameter, the second element is the new gstep parameter. The step defined for the
computation of the hessian matrix is now:
h1=max(abs(x),sqrt(gstep(1))*ones(n,1))*eps^(1/6)*gstep(2);
---
matlab/dynare_solve.m | 2 +-
matlab/global_initialization.m | 4 +++-
matlab/hessian.m | 2 +-
matlab/solve1.m | 2 +-
4 files changed, 6 insertions(+), 4 deletions(-)
diff --git a/matlab/dynare_solve.m b/matlab/dynare_solve.m
index 11310e984..921a1d53e 100644
--- a/matlab/dynare_solve.m
+++ b/matlab/dynare_solve.m
@@ -105,7 +105,7 @@ elseif options_.solve_algo == 2 || options_.solve_algo == 4
if ~jacobian_flag
fjac = zeros(nn,nn) ;
- dh = max(abs(x),options_.gstep*ones(nn,1))*eps^(1/3);
+ dh = max(abs(x),options_.gstep(1)*ones(nn,1))*eps^(1/3);
for j = 1:nn
xdh = x ;
xdh(j) = xdh(j)+dh(j) ;
diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m
index 406cee251..0d50155c4 100644
--- a/matlab/global_initialization.m
+++ b/matlab/global_initialization.m
@@ -45,7 +45,9 @@ options_.dynatol.x = 1e-5;
options_.maxit_ = 10;
options_.slowc = 1;
options_.timing = 0;
-options_.gstep = 1e-2;
+options_.gstep = ones(2,1);
+options_.gstep(1) = 1e-2;
+options_.gstep(2) = 1.0;
options_.scalv = 1;
options_.debug = 0;
options_.initval_file = 0;
diff --git a/matlab/hessian.m b/matlab/hessian.m
index 3b8fdb8b9..9e1825cd5 100644
--- a/matlab/hessian.m
+++ b/matlab/hessian.m
@@ -39,7 +39,7 @@ if ~isa(func, 'function_handle')
func = str2func(func);
end
n=size(x,1);
-h1=max(abs(x),sqrt(gstep)*ones(n,1))*eps^(1/6);
+h1=max(abs(x),sqrt(gstep(1))*ones(n,1))*eps^(1/6)*gstep(2);
h_1=h1;
xh1=x+h1;
h1=xh1-x;
diff --git a/matlab/solve1.m b/matlab/solve1.m
index 6b3dcec81..104ccf332 100644
--- a/matlab/solve1.m
+++ b/matlab/solve1.m
@@ -77,7 +77,7 @@ for its = 1:maxit
fvec = fvec(j1);
fjac = fjac(j1,j2);
else
- dh = max(abs(x(j2)),options_.gstep*ones(nn,1))*eps^(1/3);
+ dh = max(abs(x(j2)),options_.gstep(1)*ones(nn,1))*eps^(1/3);
for j = 1:nn
xdh = x ;
--
GitLab