diff --git a/matlab/DsgeSmoother.m b/matlab/DsgeSmoother.m index 488cd7434..74faa49e5 100644 --- a/matlab/DsgeSmoother.m +++ b/matlab/DsgeSmoother.m @@ -128,13 +128,13 @@ if options_.lik_init == 1 % Kalman filter kalman_algo = 1; end if options_.lyapunov_fp == 1 - Pstar = lyapunov_symm(T,R*Q*R',options_.lyapunov_fixed_point_tol,options_.lyapunov_complex_threshold, 3, [], options_.debug); + Pstar = lyapunov_symm(T,R*Q*R',options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold, 3, [], options_.debug); elseif options_.lyapunov_db == 1 Pstar = disclyap_fast(T,R*Q*R',options_.lyapunov_doubling_tol); elseif options_.lyapunov_srs == 1 - Pstar = lyapunov_symm(T,Q,options_.lyapunov_fixed_point_tol,options_.lyapunov_complex_threshold, 4, R, options_.debug); + Pstar = lyapunov_symm(T,Q,options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold, 4, R, options_.debug); else - Pstar = lyapunov_symm(T,R*Q*R',options_.qz_criterium,options_.lyapunov_complex_threshold, [], [], options_.debug); + Pstar = lyapunov_symm(T,R*Q*R',options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold, [], [], options_.debug); end; Pinf = []; elseif options_.lik_init == 2 % Old Diffuse Kalman filter @@ -171,13 +171,13 @@ elseif options_.lik_init == 5 % Old diffuse Kalman filter only for th R_tmp = R(stable, :); T_tmp = T(stable,stable); if options_.lyapunov_fp == 1 - Pstar_tmp = lyapunov_symm(T_tmp,R_tmp*Q*R_tmp',options_.lyapunov_fixed_point_tol,options_.lyapunov_complex_threshold, 3, [], options_.debug); + Pstar_tmp = lyapunov_symm(T_tmp,R_tmp*Q*R_tmp',options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold, 3, [], options_.debug); elseif options_.lyapunov_db == 1 Pstar_tmp = disclyap_fast(T_tmp,R_tmp*Q*R_tmp',options_.lyapunov_doubling_tol); elseif options_.lyapunov_srs == 1 - Pstar_tmp = lyapunov_symm(T_tmp,Q,options_.lyapunov_fixed_point_tol,options_.lyapunov_complex_threshold, 4, R_tmp, options_.debug); + Pstar_tmp = lyapunov_symm(T_tmp,Q,options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold, 4, R_tmp, options_.debug); else - Pstar_tmp = lyapunov_symm(T_tmp,R_tmp*Q*R_tmp',options_.qz_criterium,options_.lyapunov_complex_threshold, [], [], options_.debug); + Pstar_tmp = lyapunov_symm(T_tmp,R_tmp*Q*R_tmp',options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold, [], [], options_.debug); end Pstar(stable, stable) = Pstar_tmp; Pinf = []; diff --git a/matlab/UnivariateSpectralDensity.m b/matlab/UnivariateSpectralDensity.m index d2eb851e8..9aae99c4c 100644 --- a/matlab/UnivariateSpectralDensity.m +++ b/matlab/UnivariateSpectralDensity.m @@ -93,7 +93,7 @@ for i=M_.maximum_lag:-1:2 end Gamma = zeros(nvar,cutoff+1); [A,B] = kalman_transition_matrix(dr,ikx',1:nx,M_.exo_nbr); -[vx, u] = lyapunov_symm(A,B*M_.Sigma_e*B',options_.qz_criterium,options_.lyapunov_complex_threshold,[],[],options_.debug); +[vx, u] = lyapunov_symm(A,B*M_.Sigma_e*B',options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold,[],[],options_.debug); iky = iv(ivar); if ~isempty(u) iky = iky(find(any(abs(ghx(iky,:)*u) < options_.Schur_vec_tol,2))); diff --git a/matlab/dsge_likelihood.m b/matlab/dsge_likelihood.m index 53d04922c..4c663458b 100644 --- a/matlab/dsge_likelihood.m +++ b/matlab/dsge_likelihood.m @@ -370,13 +370,13 @@ switch DynareOptions.lik_init kalman_algo = 1; end if DynareOptions.lyapunov_fp == 1 - Pstar = lyapunov_symm(T,R*Q'*R',DynareOptions.lyapunov_fixed_point_tol,DynareOptions.lyapunov_complex_threshold, 3, [], DynareOptions.debug); + Pstar = lyapunov_symm(T,R*Q'*R',DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, 3, [], DynareOptions.debug); elseif DynareOptions.lyapunov_db == 1 Pstar = disclyap_fast(T,R*Q*R',DynareOptions.lyapunov_doubling_tol); elseif DynareOptions.lyapunov_srs == 1 - Pstar = lyapunov_symm(T,Q,DynareOptions.lyapunov_fixed_point_tol,DynareOptions.lyapunov_complex_threshold, 4, R, DynareOptions.debug); + Pstar = lyapunov_symm(T,Q,DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, 4, R, DynareOptions.debug); else - Pstar = lyapunov_symm(T,R*Q*R',DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, [], [], DynareOptions.debug); + Pstar = lyapunov_symm(T,R*Q*R',DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, [], [], DynareOptions.debug); end; Pinf = []; a = zeros(mm,1); @@ -472,7 +472,7 @@ switch DynareOptions.lik_init if err disp(['dsge_likelihood:: I am not able to solve the Riccati equation, so I switch to lik_init=1!']); DynareOptions.lik_init = 1; - Pstar = lyapunov_symm(T,R*Q*R',DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, [], [], DynareOptions.debug); + Pstar = lyapunov_symm(T,R*Q*R',DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, [], [], DynareOptions.debug); end Pinf = []; a = zeros(mm,1); @@ -493,13 +493,13 @@ switch DynareOptions.lik_init R_tmp = R(stable, :); T_tmp = T(stable,stable); if DynareOptions.lyapunov_fp == 1 - Pstar_tmp = lyapunov_symm(T_tmp,R_tmp*Q*R_tmp',DynareOptions.lyapunov_fixed_point_tol,DynareOptions.lyapunov_complex_threshold, 3, [], DynareOptions.debug); + Pstar_tmp = lyapunov_symm(T_tmp,R_tmp*Q*R_tmp',DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, 3, [], DynareOptions.debug); elseif DynareOptions.lyapunov_db == 1 Pstar_tmp = disclyap_fast(T_tmp,R_tmp*Q*R_tmp',DynareOptions.lyapunov_doubling_tol); elseif DynareOptions.lyapunov_srs == 1 - Pstar_tmp = lyapunov_symm(T_tmp,Q,DynareOptions.lyapunov_fixed_point_tol,DynareOptions.lyapunov_complex_threshold, 4, R_tmp, DynareOptions.debug); + Pstar_tmp = lyapunov_symm(T_tmp,Q,DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, 4, R_tmp, DynareOptions.debug); else - Pstar_tmp = lyapunov_symm(T_tmp,R_tmp*Q*R_tmp',DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, [], [], DynareOptions.debug); + Pstar_tmp = lyapunov_symm(T_tmp,R_tmp*Q*R_tmp',DynareOptions.qz_criterium,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, [], [], DynareOptions.debug); end Pstar(stable, stable) = Pstar_tmp; Pinf = []; @@ -579,14 +579,14 @@ if analytic_derivation, for i=1:EstimatedParameters.nvx k =EstimatedParameters.var_exo(i,1); DQ(k,k,i) = 2*sqrt(Q(k,k)); - dum = lyapunov_symm(T,DOm(:,:,i),DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold,[],[],DynareOptions.debug); + dum = lyapunov_symm(T,DOm(:,:,i),DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold,[],[],DynareOptions.debug); % kk = find(abs(dum) < 1e-12); % dum(kk) = 0; DP(:,:,i)=dum; if full_Hess for j=1:i, jcount=jcount+1; - dum = lyapunov_symm(T,dyn_unvech(D2Om(:,jcount)),DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold,[],[],DynareOptions.debug); + dum = lyapunov_symm(T,dyn_unvech(D2Om(:,jcount)),DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold,[],[],DynareOptions.debug); % kk = (abs(dum) < 1e-12); % dum(kk) = 0; D2P(:,jcount)=dyn_vech(dum); @@ -606,7 +606,7 @@ if analytic_derivation, offset = offset + EstimatedParameters.nvn; if DynareOptions.lik_init==1, for j=1:EstimatedParameters.np - dum = lyapunov_symm(T,DT(:,:,j+offset)*Pstar*T'+T*Pstar*DT(:,:,j+offset)'+DOm(:,:,j+offset),DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold,[],[],DynareOptions.debug); + dum = lyapunov_symm(T,DT(:,:,j+offset)*Pstar*T'+T*Pstar*DT(:,:,j+offset)'+DOm(:,:,j+offset),DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold,[],[],DynareOptions.debug); % kk = find(abs(dum) < 1e-12); % dum(kk) = 0; DP(:,:,j+offset)=dum; @@ -620,7 +620,7 @@ if analytic_derivation, D2Tij = reshape(D2T(:,jcount),size(T)); D2Omij = dyn_unvech(D2Om(:,jcount)); 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,[],[],DynareOptions.debug); + dum = lyapunov_symm(T,tmp,DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold,[],[],DynareOptions.debug); % dum(abs(dum)<1.e-12) = 0; D2P(:,jcount) = dyn_vech(dum); % D2P(:,:,j+offset,i) = dum; diff --git a/matlab/dsge_var_likelihood.m b/matlab/dsge_var_likelihood.m index a8a9e8710..ddd30bbe4 100644 --- a/matlab/dsge_var_likelihood.m +++ b/matlab/dsge_var_likelihood.m @@ -162,7 +162,7 @@ end %------------------------------------------------------------------------------ % Compute the theoretical second order moments -tmp0 = lyapunov_symm(T,R*Q*R',DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, [], [], DynareOptions.debug); +tmp0 = lyapunov_symm(T,R*Q*R',DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, [], [], DynareOptions.debug); mf = BayesInfo.mf1; % Get the non centered second order moments diff --git a/matlab/getJJ.m b/matlab/getJJ.m index de13063eb..4b8130342 100644 --- a/matlab/getJJ.m +++ b/matlab/getJJ.m @@ -64,7 +64,7 @@ else % return % end m = length(A); - GAM = lyapunov_symm(A,B*M_.Sigma_e*B',options_.qz_criterium,options_.lyapunov_complex_threshold,1,[],options_.debug); + GAM = lyapunov_symm(A,B*M_.Sigma_e*B',options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold,1,[],options_.debug); k = find(abs(GAM) < 1e-12); GAM(k) = 0; % if useautocorr, @@ -78,7 +78,7 @@ else % end % XX = lyapunov_symm_mr(A,BB,options_.qz_criterium,options_.lyapunov_complex_threshold,0); for j=1:length(indexo), - dum = lyapunov_symm(A,dOm(:,:,j),options_.qz_criterium,options_.lyapunov_complex_threshold,2,[],options_.debug); + dum = lyapunov_symm(A,dOm(:,:,j),options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold,2,[],options_.debug); % dum = XX(:,:,j); k = find(abs(dum) < 1e-12); dum(k) = 0; @@ -103,7 +103,7 @@ else end nexo = length(indexo); for j=1:length(indx), - dum = lyapunov_symm(A,dA(:,:,j+nexo)*GAM*A'+A*GAM*dA(:,:,j+nexo)'+dOm(:,:,j+nexo),options_.qz_criterium,options_.lyapunov_complex_threshold,2,[],options_.debug); + dum = lyapunov_symm(A,dA(:,:,j+nexo)*GAM*A'+A*GAM*dA(:,:,j+nexo)'+dOm(:,:,j+nexo),options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold,2,[],options_.debug); % dum = XX(:,:,j); k = find(abs(dum) < 1e-12); dum(k) = 0; diff --git a/matlab/get_variance_of_endogenous_variables.m b/matlab/get_variance_of_endogenous_variables.m index df03a1a48..ef2fd2b8a 100644 --- a/matlab/get_variance_of_endogenous_variables.m +++ b/matlab/get_variance_of_endogenous_variables.m @@ -45,7 +45,7 @@ n = length(i_var); [A,B] = kalman_transition_matrix(dr,nstatic+(1:nspred),1:nc,M_.exo_nbr); -[vx,u] = lyapunov_symm(A,B*Sigma_e*B',options_.qz_criterium,options_.lyapunov_complex_threshold, [], [], options_.debug); +[vx,u] = lyapunov_symm(A,B*Sigma_e*B',options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold, [], [], options_.debug); if size(u,2) > 0 i_stat = find(any(abs(ghx*u) < options_.Schur_vec_tol,2)); %only set those variances of objective function for which variance is finite diff --git a/matlab/lyapunov_symm.m b/matlab/lyapunov_symm.m index b694b1354..81208c358 100644 --- a/matlab/lyapunov_symm.m +++ b/matlab/lyapunov_symm.m @@ -1,14 +1,12 @@ -function [x,u] = lyapunov_symm(a,b,third_argument,lyapunov_complex_threshold,method, R, debug) % --*-- Unitary tests --*-- +function [x,u] = lyapunov_symm(a,b,lyapunov_fixed_point_tol,qz_criterium,lyapunov_complex_threshold,method, R, debug) % --*-- Unitary tests --*-- % 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. -% third_argument [double] scalar, if method <= 2 : -% qz_criterium = third_argument unit root threshold for eigenvalues in a, -% elseif method = 3 : -% tol =third_argument the convergence criteria for fixed_point algorithm. +% qz_criterium [double] unit root threshold for eigenvalues +% lyapunov_fixed_point_tol [double] convergence criteria for fixed_point algorithm. % 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 @@ -44,11 +42,11 @@ function [x,u] = lyapunov_symm(a,b,third_argument,lyapunov_complex_threshold,met % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . -if nargin<5 || isempty(method) +if nargin<6 || isempty(method) method = 0; end -if nargin<7 +if nargin<8 debug = 0; end @@ -57,7 +55,6 @@ if method == 3 if ~isempty(method1) method = method1; end; - tol = third_argument; if debug fprintf('lyapunov_symm:: [method=%d] \n',method); end @@ -73,7 +70,7 @@ if method == 3 end; at = a'; %fixed point iterations - while evol > tol && it_fp < max_it_fp; + while evol > lyapunov_fixed_point_tol && it_fp < max_it_fp; X_old = X; X = a * X * at + b; evol = max(sum(abs(X - X_old))); %norm_1 @@ -110,7 +107,6 @@ elseif method == 4 return; end; -qz_criterium = third_argument; if method persistent U T k n else @@ -210,8 +206,8 @@ u = U(:,1:k); %$ %$ % DynareOptions.lyapunov_fp == 1 %$ try -%$ Pstar1_small = lyapunov_symm(T_small,R_small*Q_small*R_small',lyapunov_fixed_point_tol,lyapunov_fixed_point_tol,3, [], 0); -%$ Pstar1_large = lyapunov_symm(T_large,R_large*Q_large*R_large',lyapunov_fixed_point_tol,lyapunov_fixed_point_tol,3, [], 0); +%$ Pstar1_small = lyapunov_symm(T_small,R_small*Q_small*R_small',lyapunov_fixed_point_tol,qz_criterium,lyapunov_complex_threshold,3, [], 0); +%$ Pstar1_large = lyapunov_symm(T_large,R_large*Q_large*R_large',lyapunov_fixed_point_tol,qz_criterium,lyapunov_complex_threshold,3, [], 0); %$ t(1) = 1; %$ catch %$ t(1) = 0; @@ -229,8 +225,8 @@ u = U(:,1:k); %$ % Dynareoptions.lyapunov_srs == 1 %$ if (isoctave && user_has_octave_forge_package('control')) || (~isoctave && user_has_matlab_license('control_toolbox')) %$ try -%$ Pstar3_small = lyapunov_symm(T_small,Q_small,lyapunov_fixed_point_tol,lyapunov_complex_threshold,4,R_small,0); -%$ Pstar3_large = lyapunov_symm(T_large,Q_large,lyapunov_fixed_point_tol,lyapunov_complex_threshold,4,R_large,0); +%$ Pstar3_small = lyapunov_symm(T_small,Q_small,lyapunov_fixed_point_tol,qz_criterium,lyapunov_complex_threshold,4,R_small,0); +%$ Pstar3_large = lyapunov_symm(T_large,Q_large,lyapunov_fixed_point_tol,qz_criterium,lyapunov_complex_threshold,4,R_large,0); %$ t(3) = 1; %$ catch %$ t(3) = 0; @@ -241,8 +237,8 @@ u = U(:,1:k); %$ %$ % Standard %$ try -%$ Pstar4_small = lyapunov_symm(T_small,R_small*Q_small*R_small',qz_criterium, lyapunov_complex_threshold, [], [], 0); -%$ Pstar4_large = lyapunov_symm(T_large,R_large*Q_large*R_large',qz_criterium, lyapunov_complex_threshold, [], [], 0); +%$ Pstar4_small = lyapunov_symm(T_small,R_small*Q_small*R_small',lyapunov_fixed_point_tol,qz_criterium, lyapunov_complex_threshold, [], [], 0); +%$ Pstar4_large = lyapunov_symm(T_large,R_large*Q_large*R_large',lyapunov_fixed_point_tol,qz_criterium, lyapunov_complex_threshold, [], [], 0); %$ t(4) = 1; %$ catch %$ t(4) = 0; diff --git a/matlab/non_linear_dsge_likelihood.m b/matlab/non_linear_dsge_likelihood.m index a9348034a..d80f61cc2 100644 --- a/matlab/non_linear_dsge_likelihood.m +++ b/matlab/non_linear_dsge_likelihood.m @@ -271,7 +271,7 @@ ReducedForm.mf1 = mf1; switch DynareOptions.particle.initialization case 1% Initial state vector covariance is the ergodic variance associated to the first order Taylor-approximation of the model. StateVectorMean = ReducedForm.constant(mf0); - StateVectorVariance = lyapunov_symm(ReducedForm.ghx(mf0,:),ReducedForm.ghu(mf0,:)*ReducedForm.Q*ReducedForm.ghu(mf0,:)',1e-12,1e-12,[],[],DynareOptions.debug); + StateVectorVariance = lyapunov_symm(ReducedForm.ghx(mf0,:),ReducedForm.ghu(mf0,:)*ReducedForm.Q*ReducedForm.ghu(mf0,:)',DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold,[],[],DynareOptions.debug); case 2% Initial state vector covariance is a monte-carlo based estimate of the ergodic variance (consistent with a k-order Taylor-approximation of the model). StateVectorMean = ReducedForm.constant(mf0); old_DynareOptionsperiods = DynareOptions.periods; diff --git a/matlab/partial_information/disclyap_fast.m b/matlab/partial_information/disclyap_fast.m index 2f748b6d2..461efc26f 100644 --- a/matlab/partial_information/disclyap_fast.m +++ b/matlab/partial_information/disclyap_fast.m @@ -1,6 +1,14 @@ -function X=disclyap_fast(G,V,tol,ch) +function [X,exitflag]=disclyap_fast(G,V,tol,ch) % function X=disclyap_fast(G,V,ch) -% +% Inputs: +% - G [double] first input matrix +% - V [double] second input matrix +% - tol [scalar] tolerance criterion +% - ch empty of non-empty if non-empty: check positive-definiteness +% Outputs: +% - X [double] solution matrix +% - exitflag [scalar] 0 if solution is found, 1 otherwise +% % Solve the discrete Lyapunov Equation % X=G*X*G'+V % Using the Doubling Algorithm @@ -11,7 +19,7 @@ function X=disclyap_fast(G,V,tol,ch) % Joe Pearlman and Alejandro Justiniano % 3/5/2005 -% Copyright (C) 2010-2012 Dynare Team +% Copyright (C) 2010-2015 Dynare Team % % This file is part of Dynare. % @@ -34,6 +42,7 @@ else flag_ch = 1; end s=size(G,1); +exitflag=0; %tol = 1e-16; @@ -41,13 +50,20 @@ P0=V; A0=G; matd=1; -while matd > tol +iter=1; +while matd > tol && iter< 2000 P1=P0+A0*P0*A0'; A1=A0*A0; matd=max( max( abs( P1 - P0 ) ) ); P0=P1; A0=A1; + iter=iter+1; end +if iter==5000 + X=NaN(P0); + exitflag=1; + return +end clear A0 A1 P1; X=(P0+P0')/2; @@ -56,6 +72,7 @@ X=(P0+P0')/2; if flag_ch==1 [C,p]=chol(X); if p ~= 0 + exitflag=1; error('X is not positive definite') end end \ No newline at end of file diff --git a/matlab/th_autocovariances.m b/matlab/th_autocovariances.m index eb3f8dc8a..6852f9654 100644 --- a/matlab/th_autocovariances.m +++ b/matlab/th_autocovariances.m @@ -129,7 +129,7 @@ end; % and compute 2nd order mean correction on stationary variables (in case of % HP filtering, this mean correction is computed *before* filtering) if options_.order == 2 || options_.hp_filter == 0 - [vx, u] = lyapunov_symm(A,B*M_.Sigma_e*B',options_.qz_criterium,options_.lyapunov_complex_threshold,[],[],options_.debug); + [vx, u] = lyapunov_symm(A,B*M_.Sigma_e*B',options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold,[],[],options_.debug); if options_.block == 0 iky = inv_order_var(ivar); else @@ -188,11 +188,11 @@ if options_.hp_filter == 0 b1 = b1*cs; b2(:,exo_names_orig_ord) = ghu(iky,:); b2 = b2*cs; - vx = lyapunov_symm(A,b1*b1',options_.qz_criterium,options_.lyapunov_complex_threshold,1,[],options_.debug); + vx = lyapunov_symm(A,b1*b1',options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold,1,[],options_.debug); vv = diag(aa*vx*aa'+b2*b2'); vv2 = 0; for i=1:M_.exo_nbr - vx1 = lyapunov_symm(A,b1(:,i)*b1(:,i)',options_.qz_criterium,options_.lyapunov_complex_threshold,2,[],options_.debug); + vx1 = lyapunov_symm(A,b1(:,i)*b1(:,i)',options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold,2,[],options_.debug); vx2 = abs(diag(aa*vx1*aa'+b2(:,i)*b2(:,i)')); Gamma_y{nar+2}(stationary_vars,i) = vx2; vv2 = vv2 +vx2; diff --git a/matlab/thet2tau.m b/matlab/thet2tau.m index 9a3214344..25a489432 100644 --- a/matlab/thet2tau.m +++ b/matlab/thet2tau.m @@ -51,7 +51,7 @@ elseif flagmoments==-1 tau=[oo_.dr.ys(oo_.dr.order_var); g1(:)]; else - GAM = lyapunov_symm(A,B*M_.Sigma_e*B',options_.qz_criterium,options_.lyapunov_complex_threshold,[],[],options_.debug); + GAM = lyapunov_symm(A,B*M_.Sigma_e*B',options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold,[],[],options_.debug); k = find(abs(GAM) < 1e-12); GAM(k) = 0; if useautocorr, diff --git a/tests/optimal_policy/mult_elimination_test.mod b/tests/optimal_policy/mult_elimination_test.mod index 69f064601..1a9b2117b 100644 --- a/tests/optimal_policy/mult_elimination_test.mod +++ b/tests/optimal_policy/mult_elimination_test.mod @@ -60,7 +60,7 @@ V0 = oo_.var(k4,k4); Atest = [dr2.M1(k3,:) dr2.M2(k3,:) dr2.M4(k3,:); eye(6) zeros(6,10);zeros(4,16)]; Btest = [dr2.M3(k3,:); zeros(6,4); eye(4)]; -V1=lyapunov_symm(Atest,Btest*M_.Sigma_e*Btest',1+1e-6,options_.lyapunov_complex_threshold); +V1=lyapunov_symm(Atest,Btest*M_.Sigma_e*Btest',options_.lyapunov_fixed_point_tol,1+1e-6,options_.lyapunov_complex_threshold); if max(max(abs(V1(1:6,1:6)-V0))) disp('Test OK')