diff --git a/matlab/DsgeLikelihood.m b/matlab/DsgeLikelihood.m
index 14b1fbd7ccc8af2e6f115cd0c3439f829498d055..e879c609fe0e9653c4bb4d61eeacb1cb75f81fe7 100644
--- a/matlab/DsgeLikelihood.m
+++ b/matlab/DsgeLikelihood.m
@@ -170,7 +170,7 @@ function [fval,cost_flag,ys,trend_coeff,info] = DsgeLikelihood(xparam1,gend,data
if kalman_algo ~= 2
kalman_algo = 1;
end
- Pstar = lyapunov_symm(T,R*Q*R',options_.qz_criterium);
+ Pstar = lyapunov_symm(T,R*Q*R',options_.qz_criterium,options_.lyapunov_complex_threshold);
Pinf = [];
elseif options_.lik_init == 2 % Old Diffuse Kalman filter
if kalman_algo ~= 2
diff --git a/matlab/DsgeLikelihood_hh.m b/matlab/DsgeLikelihood_hh.m
index 6838c9ceadda3d972174d8ddc2da8e98355ef256..82d1a37849d3b371e53ebec6474d6699d76b22fd 100644
--- a/matlab/DsgeLikelihood_hh.m
+++ b/matlab/DsgeLikelihood_hh.m
@@ -171,7 +171,7 @@ function [fval,llik,cost_flag,ys,trend_coeff,info] = DsgeLikelihood_hh(xparam1,g
if kalman_algo ~= 2
kalman_algo = 1;
end
- Pstar = lyapunov_symm(T,R*Q*R',options_.qz_criterium);
+ Pstar = lyapunov_symm(T,R*Q*R',options_.qz_criterium,options_.lyapunov_complex_threshold);
Pinf = [];
elseif options_.lik_init == 2 % Old Diffuse Kalman filter
if kalman_algo ~= 2
diff --git a/matlab/DsgeSmoother.m b/matlab/DsgeSmoother.m
index 8eba116f10efb623714d16f3083a26b91d0ed1b9..096e53ae1b8ae81f4090dd36d71ddd9cecb02b00 100644
--- a/matlab/DsgeSmoother.m
+++ b/matlab/DsgeSmoother.m
@@ -113,7 +113,7 @@ function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T,R,P,PK,d,
if kalman_algo ~= 2
kalman_algo = 1;
end
- Pstar = lyapunov_symm(T,R*Q*transpose(R),options_.qz_criterium);
+ Pstar = lyapunov_symm(T,R*Q*transpose(R),options_.qz_criterium,options_.lyapunov_complex_threshold);
Pinf = [];
elseif options_.lik_init == 2 % Old Diffuse Kalman filter
if kalman_algo ~= 2
diff --git a/matlab/DsgeVarLikelihood.m b/matlab/DsgeVarLikelihood.m
index bff61250a5b1c6d97783606cba7ec88322c45130..a31ab9dc4b07a5963ba73128a1b3efba46ab9b8b 100644
--- a/matlab/DsgeVarLikelihood.m
+++ b/matlab/DsgeVarLikelihood.m
@@ -137,7 +137,7 @@ end
%------------------------------------------------------------------------------
% 3. theoretical moments (second order)
%------------------------------------------------------------------------------
-tmp0 = lyapunov_symm(T,R*Q*R',options_.qz_criterium);% I compute the variance-covariance matrix
+tmp0 = lyapunov_symm(T,R*Q*R',options_.qz_criterium,options_.lyapunov_complex_threshold);% I compute the variance-covariance matrix
mf = bayestopt_.mf1; % of the restricted state vector.
% Get the non centered second order moments
diff --git a/matlab/UnivariateSpectralDensity.m b/matlab/UnivariateSpectralDensity.m
index d8e1a630ed402ff6f0d26bbcf17a839bd65dae64..818063ad91712766d05d5d9ea81553e1b374e76e 100644
--- a/matlab/UnivariateSpectralDensity.m
+++ b/matlab/UnivariateSpectralDensity.m
@@ -96,7 +96,7 @@ for i=M_.maximum_lag:-1:2
end
Gamma = zeros(nvar,cutoff+1);
[A,B] = kalman_transition_matrix(dr,ikx',1:nx,dr.transition_auxiliary_variables,M_.exo_nbr);
-[vx, u] = lyapunov_symm(A,B*M_.Sigma_e*B',options_.qz_criterium);
+[vx, u] = lyapunov_symm(A,B*M_.Sigma_e*B',options_.qz_criterium,options_.lyapunov_complex_threshold);
iky = iv(ivar);
if ~isempty(u)
iky = iky(find(any(abs(ghx(iky,:)*u) < options_.Schur_vec_tol,2)));
diff --git a/matlab/calib_obj.m b/matlab/calib_obj.m
index 28c79686aa90695f5e4cc7a875ec08f2eb40c773..cd44b61eba138cb0171b7b3f13a58d262a2013c8 100644
--- a/matlab/calib_obj.m
+++ b/matlab/calib_obj.m
@@ -28,11 +28,11 @@ function f=calib_obj(M_.Sigma_e,A,ghu1,ghx,ghu,targets,var_weights,iy,nar)
b=ghu1*M_.Sigma_e*ghu1';
vx = [];
if isempty(vx)
- vx = lyapunov_symm(A,b,options_.qz_criterium);
+ vx = lyapunov_symm(A,b,options_.qz_criterium,options_.lyapunov_complex_threshold);
else
[vx,status] = bicgstab_(@f_var,b(:),vx(:),1e-8,50,A,nx);
if status
- vx = lyapunov_symm(A,b,options_.qz_criterium);
+ vx = lyapunov_symm(A,b,options_.qz_criterium,options_.lyapunov_complex_threshold);
else
vx=reshape(vx,nx,nx);
end
diff --git a/matlab/calib_obj2.m b/matlab/calib_obj2.m
index af0a034929de9bf41ce5d98e2e71b2aef0cd3cd9..0e3e8a03eec91b6f27f5943c5b7220f7e290e715 100644
--- a/matlab/calib_obj2.m
+++ b/matlab/calib_obj2.m
@@ -26,7 +26,7 @@ function objective=calib_obj2(M_.Sigma_e,A,ghu1,ghx,ghu,targets,var_weights,iy,n
M_.Sigma_e=diag(M_.Sigma_e);
nx = size(ghx,2);
b=ghu1*M_.Sigma_e*ghu1';
- vx = lyapunov_symm(A,b,options_.qz_criterium);
+ vx = lyapunov_symm(A,b,options_.qz_criterium,options_.lyapunov_complex_threshold);
oo_.gamma_y{1} = ghx*vx*ghx'+ ghu*M_.Sigma_e*ghu';
if ~isempty(targets{1})
objective{1} = sqrt(oo_.gamma_y{1}(iy{1}));
diff --git a/matlab/dsge_posterior_kernel.m b/matlab/dsge_posterior_kernel.m
index 4bf7e500d396340ebe896ff637eefe297e0c6487..f6d52d2b562240d4985341d98bb457c27cbcf5fb 100644
--- a/matlab/dsge_posterior_kernel.m
+++ b/matlab/dsge_posterior_kernel.m
@@ -170,7 +170,7 @@ function [fval,cost_flag,ys,trend_coeff,info] = DsgeLikelihood(xparam1,gend,data
if kalman_algo ~= 2
kalman_algo = 1;
end
- Pstar = lyapunov_symm(T,R*Q*R',options_.qz_criterium);
+ Pstar = lyapunov_symm(T,R*Q*R',options_.qz_criterium,options_.lyapunov_complex_threshold);
Pinf = [];
elseif options_.lik_init == 2 % Old Diffuse Kalman filter
if kalman_algo ~= 2
diff --git a/matlab/get_variance_of_endogenous_variables.m b/matlab/get_variance_of_endogenous_variables.m
index 7d509aaa2279b9bf6793031a8162e6cf0430d6a0..21542d2c4d2325bca205ca0ef3358bf056d16eb5 100644
--- a/matlab/get_variance_of_endogenous_variables.m
+++ b/matlab/get_variance_of_endogenous_variables.m
@@ -45,7 +45,7 @@ function [vx1,i_ns] = get_variance_of_endogenous_variables(dr,i_var)
[A,B] = kalman_transition_matrix(dr,nstatic+(1:npred),1:nc,dr.transition_auxiliary_variables,M_.exo_nbr);
- [vx,u] = lyapunov_symm(A,B*Sigma_e*B',options_.qz_criterium);
+ [vx,u] = lyapunov_symm(A,B*Sigma_e*B',options_.qz_criterium,options_.lyapunov_complex_threshold);
if size(u,2) > 0
i_stat = find(any(abs(ghx*u) < options_.Schur_vec_tol,2));
diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m
index 35a471566e7dfd2cd83fb80817ace59e93cea09e..e646754da0252e43697029578e97cc03f4b5f38c 100644
--- a/matlab/global_initialization.m
+++ b/matlab/global_initialization.m
@@ -44,6 +44,7 @@ function global_initialization()
options_.Schur_vec_tol = 1e-8; % used to find nonstationary variables
% in Schur decomposition of the
% transition matrix
+ options_.lyapunov_complex_threshold = 1e-15;
options_.solve_tolf = eps^(1/3);
options_.solve_tolx = eps^(2/3);
options_.solve_maxit = 500;
@@ -112,7 +113,7 @@ function global_initialization()
options_.drop = 100;
options_.simul_algo = 0;
options_.model_mode = 0;
- % if mjdgges.dll (or .mexw32 or ....) doesn't exit, matlab/qz is added to the path.
+ % if mjdgges.dll (or .mexw32 or ....) doesn't exist, matlab/qz is added to the path.
% There exists now qz/mjdgges.m that contains the calls to the old Sims code
% Hence, if mjdgges.m is visible exist(...)==2,
% this means that the DLL isn't avaiable and use_qzdiv is set to 1
@@ -151,6 +152,7 @@ function global_initialization()
options_.loglinear = 0;
options_.markowitz = 0.5;
options_.mh_conf_sig = 0.90;
+ options_.prior_interval = 0.90;
options_.mh_drop = 0.5;
options_.mh_jscale = 0.2;
options_.mh_init_scale = 2*options_.mh_jscale;
diff --git a/matlab/lyapunov_symm.m b/matlab/lyapunov_symm.m
index 35dad11f08f570e99527842713092719d42b3186..f5c3a6d9384fdf798ea300041a9cfe4f52944345 100644
--- a/matlab/lyapunov_symm.m
+++ b/matlab/lyapunov_symm.m
@@ -1,13 +1,17 @@
-function [x,u] = lyapunov_symm(a,b,qz_criterium,method)
+function [x,u] = lyapunov_symm(a,b,qz_criterium,lyapunov_complex_threshold,method)
% Solves the Lyapunov equation x-a*x*a' = b, for b and x symmetric matrices.
% If a has some unit roots, the function computes only the solution of the stable subsystem.
%
% INPUTS:
-% a [double] n*n matrix.
-% b [double] n*n matrix
-% method [integer] Scalar, if method=0 [default] then U, T, n and k are not persistent.
-% method=1 then U, T, n and k are declared as persistent variables and the schur decomposition is triggered.
-% method=2 then U, T, n and k are declared as persistent variables and the schur decomposition is not performed.
+% a [double] n*n matrix.
+% b [double] n*n matrix.
+% qz_criterium [double] scalar, unit root threshold for eigenvalues in a.
+% lyapunov_complex_threshold [double] scalar, complex block threshold for the upper triangular matrix T.
+% method [integer] Scalar, if method=0 [default] then U, T, n and k are not persistent.
+% method=1 then U, T, n and k are declared as persistent
+% variables and the schur decomposition is triggered.
+% method=2 then U, T, n and k are declared as persistent
+% variables and the schur decomposition is not performed.
% OUTPUTS
% x: [double] m*m solution matrix of the lyapunov equation, where m is the dimension of the stable subsystem.
% u: [double] Schur vectors associated with unit roots
@@ -34,7 +38,7 @@ function [x,u] = lyapunov_symm(a,b,qz_criterium,method)
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
- if nargin<4
+ if nargin<5
method = 0;
end
if method
@@ -64,25 +68,21 @@ function [x,u] = lyapunov_symm(a,b,qz_criterium,method)
if exist('ordschur','builtin')
% Selects stable roots
[U,T] = ordschur(U,T,e1);
- % B = U(:,k+1:end)'*b*U(:,k+1:end);
T = T(k+1:end,k+1:end);
elseif k > 0
% Problem for Matlab version that don't have ordschur
error(['lyapunov_sym: you need a Matlab version > 6.5 to handle models' ...
' with unit roots'])
- %else
- % No unit root
- %n = length(e1);
- %B=u'*b*u;
end
end
- B = U(:,k+1:end)'*b*U(:,k+1:end);
+
+ B = U(:,k+1:end)'*b*U(:,k+1:end);
x = zeros(n,n);
i = n;
-
+
while i >= 2
- if T(i,i-1) == 0
+ if abs(T(i,i-1))<lyapunov_complex_threshold
if i == n
c = zeros(n,1);
else
@@ -105,8 +105,8 @@ function [x,u] = lyapunov_symm(a,b,qz_criterium,method)
T(i-1,i-1)*T(1:i,i+1:end)*x(i+1:end,i-1) + ...
T(i-1,i)*T(1:i,i+1:end)*x(i+1:end,i);
end
- q = [eye(i)-T(1:i,1:i)*T(i,i) - T(1:i,1:i)*T(i,i-1) ; ...
- -T(1:i,1:i)*T(i-1,i) eye(i) - T(1:i,1:i)*T(i-1,i-1) ];
+ q = [ eye(i)-T(1:i,1:i)*T(i,i) , -T(1:i,1:i)*T(i,i-1) ; ...
+ -T(1:i,1:i)*T(i-1,i) , eye(i)-T(1:i,1:i)*T(i-1,i-1) ];
z = q\[ B(1:i,i)+c ; B(1:i,i-1) + c1 ];
x(1:i,i) = z(1:i);
x(1:i,i-1) = z(i+1:end);
diff --git a/matlab/th_autocovariances.m b/matlab/th_autocovariances.m
index 7f6caa8e152e3c01d8c667c07e2f38cf87e16d56..b60a692cb28891b10dc4ee94f8276ae4a9a2ec53 100644
--- a/matlab/th_autocovariances.m
+++ b/matlab/th_autocovariances.m
@@ -93,7 +93,7 @@ function [Gamma_y,ivar] = th_autocovariances(dr,ivar,M_,options_,nodecomposition
% state space representation for state variables only
[A,B] = kalman_transition_matrix(dr,ipred,1:nx,dr.transition_auxiliary_variables,M_.exo_nbr);
if options_.order == 2 | options_.hp_filter == 0
- [vx, u] = lyapunov_symm(A,B*M_.Sigma_e*B',options_.qz_criterium);
+ [vx, u] = lyapunov_symm(A,B*M_.Sigma_e*B',options_.qz_criterium,options_.lyapunov_complex_threshold);
iky = iv(ivar);
if ~isempty(u)
iky = iky(find(all(abs(ghx(iky,:)*u) < options_.Schur_vec_tol,2)));
@@ -132,10 +132,10 @@ function [Gamma_y,ivar] = th_autocovariances(dr,ivar,M_,options_,nodecomposition
b1 = b1*cs;
b2(:,exo_names_orig_ord) = ghu(iky,:);
b2 = b2*cs;
- vx = lyapunov_symm(A,b1*b1',options_.qz_criterium,1);
+ vx = lyapunov_symm(A,b1*b1',options_.qz_criterium,options_.lyapunov_complex_threshold,1);
vv = diag(aa*vx*aa'+b2*b2');
for i=1:M_.exo_nbr
- vx1 = lyapunov_symm(A,b1(:,i)*b1(:,i)',options_.qz_criterium,2);
+ vx1 = lyapunov_symm(A,b1(:,i)*b1(:,i)',options_.qz_criterium,options_.lyapunov_complex_threshold,2);
Gamma_y{nar+2}(:,i) = abs(diag(aa*vx1*aa'+b2(:,i)*b2(:,i)'))./vv;
end
end