From 2c5b1fed2d43ffdeb60207295e0a03610864b487 Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Sun, 18 Dec 2016 18:23:42 +0100 Subject: [PATCH] Use master function lyapunov_solver.m to call individual solvers --- matlab/DsgeSmoother.m | 20 +- matlab/UnivariateSpectralDensity.m | 2 +- matlab/dsge_likelihood.m | 32 +-- matlab/dsge_var_likelihood.m | 2 +- matlab/getJJ.m | 6 +- matlab/get_variance_of_endogenous_variables.m | 2 +- matlab/lyapunov_solver.m | 186 ++++++++++++++++++ matlab/lyapunov_symm.m | 140 +------------ matlab/non_linear_dsge_likelihood.m | 2 +- matlab/partial_information/disclyap_fast.m | 23 +-- matlab/th_autocovariances.m | 6 +- matlab/thet2tau.m | 2 +- 12 files changed, 223 insertions(+), 200 deletions(-) create mode 100644 matlab/lyapunov_solver.m diff --git a/matlab/DsgeSmoother.m b/matlab/DsgeSmoother.m index 4dd5711af..f9e459e00 100644 --- a/matlab/DsgeSmoother.m +++ b/matlab/DsgeSmoother.m @@ -156,15 +156,7 @@ if options_.lik_init == 1 % Kalman filter if kalman_algo ~= 2 kalman_algo = 1; end - if options_.lyapunov_fp == 1 - 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_.qz_criterium,options_.lyapunov_complex_threshold, 4, R, options_.debug); - else - Pstar = lyapunov_symm(T,R*Q*R',options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold, [], [], options_.debug); - end; + Pstar=lyapunov_solver(T,R,Q,options_); Pinf = []; elseif options_.lik_init == 2 % Old Diffuse Kalman filter if kalman_algo ~= 2 @@ -210,15 +202,7 @@ elseif options_.lik_init == 5 % Old diffuse Kalman filter only for th end 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_.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_.qz_criterium,options_.lyapunov_complex_threshold, 4, R_tmp, options_.debug); - else - 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_tmp=lyapunov_solver(T_tmp,R_tmp,Q,DynareOptions); Pstar(stable, stable) = Pstar_tmp; Pinf = []; end diff --git a/matlab/UnivariateSpectralDensity.m b/matlab/UnivariateSpectralDensity.m index 7a69585cb..fa0a9e22c 100644 --- a/matlab/UnivariateSpectralDensity.m +++ b/matlab/UnivariateSpectralDensity.m @@ -95,7 +95,7 @@ for i=M_.maximum_lag:-1:2 end [A,B] = kalman_transition_matrix(oo_.dr,ikx',1:nx,M_.exo_nbr); -[vx, u] = lyapunov_symm(A,B*M_.Sigma_e*B',options_.lyapunov_fixed_point_tol,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 f14682115..8a4a7d6b6 100644 --- a/matlab/dsge_likelihood.m +++ b/matlab/dsge_likelihood.m @@ -111,7 +111,7 @@ function [fval,info,exit_flag,DLIK,Hess,SteadyState,trend_coeff,Model,DynareOpti %! @sp 2 %! @strong{This function calls:} %! @sp 1 -%! @ref{dynare_resolve}, @ref{lyapunov_symm}, @ref{compute_Pinf_Pstar}, @ref{kalman_filter_d}, @ref{missing_observations_kalman_filter_d}, @ref{univariate_kalman_filter_d}, @ref{kalman_steady_state}, @ref{getH}, @ref{kalman_filter}, @ref{score}, @ref{AHessian}, @ref{missing_observations_kalman_filter}, @ref{univariate_kalman_filter}, @ref{priordens} +%! @ref{dynare_resolve}, @ref{lyapunov_symm}, @ref{lyapunov_solver}, @ref{compute_Pinf_Pstar}, @ref{kalman_filter_d}, @ref{missing_observations_kalman_filter_d}, @ref{univariate_kalman_filter_d}, @ref{kalman_steady_state}, @ref{getH}, @ref{kalman_filter}, @ref{score}, @ref{AHessian}, @ref{missing_observations_kalman_filter}, @ref{univariate_kalman_filter}, @ref{priordens} %! @end deftypefn %@eod: @@ -351,15 +351,7 @@ switch DynareOptions.lik_init % Use standard kalman filter except if the univariate filter is explicitely choosen. kalman_algo = 1; end - if DynareOptions.lyapunov_fp == 1 - 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.qz_criterium,DynareOptions.lyapunov_complex_threshold, 4, R, DynareOptions.debug); - else - Pstar = lyapunov_symm(T,R*Q*R',DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, [], [], DynareOptions.debug); - end; + Pstar=lyapunov_solver(T,R,Q,DynareOptions); Pinf = []; a = zeros(mm,1); Zflag = 0; @@ -471,7 +463,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.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, [], [], DynareOptions.debug); + Pstar=lyapunov_solver(T,R,Q,DynareOptions); end Pinf = []; a = zeros(mm,1); @@ -491,15 +483,7 @@ switch DynareOptions.lik_init end 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.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.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.qz_criterium,DynareOptions.lyapunov_complex_threshold, [], [], DynareOptions.debug); - end + Pstar_tmp=lyapunov_solver(T_tmp,R_tmp,Q,DynareOptions); Pstar(stable, stable) = Pstar_tmp; Pinf = []; a = zeros(mm,1); @@ -580,14 +564,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.lyapunov_fixed_point_tol,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.lyapunov_fixed_point_tol,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); @@ -607,7 +591,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.lyapunov_fixed_point_tol,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; @@ -621,7 +605,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.lyapunov_fixed_point_tol,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 ca3a86d2f..646120c5f 100644 --- a/matlab/dsge_var_likelihood.m +++ b/matlab/dsge_var_likelihood.m @@ -187,7 +187,7 @@ end %------------------------------------------------------------------------------ % Compute the theoretical second order moments -tmp0 = lyapunov_symm(T,R*Q*R',DynareOptions.lyapunov_fixed_point_tol,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 51dcf32e2..b7686a035 100644 --- a/matlab/getJJ.m +++ b/matlab/getJJ.m @@ -102,7 +102,7 @@ else % return % end m = length(A); - GAM = lyapunov_symm(A,B*M_.Sigma_e*B',options_.lyapunov_fixed_point_tol,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, @@ -116,7 +116,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_.lyapunov_fixed_point_tol,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; @@ -141,7 +141,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_.lyapunov_fixed_point_tol,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 ef2fd2b8a..005ceb8e0 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_.lyapunov_fixed_point_tol,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_solver.m b/matlab/lyapunov_solver.m new file mode 100644 index 000000000..a7b96c997 --- /dev/null +++ b/matlab/lyapunov_solver.m @@ -0,0 +1,186 @@ +function P=lyapunov_solver(T,R,Q,DynareOptions) % --*-- Unitary tests --*-- +% function P=lyapunov_solver(T,R,Q,DynareOptions) +% Solves the Lyapunov equation P-T*P*T' = R*Q*R' arising in a state-space +% system, where P is the variance of the states +% +% Inputs +% T [double] n*n matrix. +% R [double] n*m matrix. +% Q [double] m*m matrix. +% DynareOptions [structure] Dynare options +% +% Outputs +% P [double] n*n matrix. +% +% Algorithms +% Default, if none of the other algorithms is selected: +% Reordered Schur decomposition (Bartels-Stewart algorithm) +% DynareOptions.lyapunov_fp == 1 +% iteration-based fixed point algorithm +% DynareOptions.lyapunov_db == 1 +% doubling algorithm +% DynareOptions.lyapunov_srs == 1 +% Square-root solver for discrete-time Lyapunov equations (requires Matlab System Control toolbox +% or Octave control package) + +% Copyright (C) 2016 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 . + +if DynareOptions.lyapunov_fp == 1 + P = 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 + [P, errorflag] = disclyap_fast(T,R*Q*R',DynareOptions.lyapunov_doubling_tol); + if errorflag %use Schur-based method + P = lyapunov_symm(T,R*Q*R',DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, [], DynareOptions.debug); + end +elseif DynareOptions.lyapunov_srs == 1 + % works only with Matlab System Control toolbox or Octave control package, + if isoctave + if ~user_has_octave_forge_package('control') + error('lyapunov=square_root_solver not available; you must install the control package from Octave Forge') + end + else + if ~user_has_matlab_license('control_toolbox') + error('lyapunov=square_root_solver not available; you must install the control system toolbox') + end + end + chol_Q = R*chol(Q,'lower'); + R_P = dlyapchol(T,chol_Q); + P = R_P' * R_P; +else + P = lyapunov_symm(T,R*Q*R',DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, [], DynareOptions.debug); +end; + +%@test:1 +%$ t = NaN(10,1); +%$ options_.lyapunov_complex_threshold = 1e-15; +%$ options_.qz_zero_threshold = 1e-6; +%$ options_.qz_criterium=1-options_.qz_zero_threshold; +%$ options_.lyapunov_fixed_point_tol = 1e-10; +%$ options_.lyapunov_doubling_tol = 1e-16; +%$ options_.debug=0; +%$ +%$ n_small=8; +%$ m_small=10; +%$ T_small=randn(n_small,n_small); +%$ T_small=0.99*T_small/max(abs(eigs(T_small))); +%$ tmp2=randn(m_small,m_small); +%$ Q_small=tmp2*tmp2'; +%$ R_small=randn(n_small,m_small); +%$ +%$ n_large=9; +%$ m_large=11; +%$ T_large=randn(n_large,n_large); +%$ T_large=0.99*T_large/max(abs(eigs(T_large))); +%$ tmp2=randn(m_large,m_large); +%$ Q_large=tmp2*tmp2'; +%$ R_large=randn(n_large,m_large); +%$ +%$ % DynareOptions.lyapunov_fp == 1 +%$ options_.lyapunov_fp = 1; +%$ try +%$ Pstar1_small = lyapunov_solver(T_small,R_small,Q_small,options_); +%$ Pstar1_large = lyapunov_solver(T_large,R_large,Q_large,options_); +%$ t(1) = 1; +%$ catch +%$ t(1) = 0; +%$ end +%$ +%$ % Dynareoptions.lyapunov_db == 1 +%$ options_.lyapunov_fp = 0; +%$ options_.lyapunov_db = 1; +%$ try +%$ Pstar2_small = lyapunov_solver(T_small,R_small,Q_small,options_); +%$ Pstar2_large = lyapunov_solver(T_large,R_large,Q_large,options_); +%$ t(2) = 1; +%$ catch +%$ t(2) = 0; +%$ end +%$ +%$ % Dynareoptions.lyapunov_srs == 1 +%$ if (isoctave && user_has_octave_forge_package('control')) || (~isoctave && user_has_matlab_license('control_toolbox')) +%$ options_.lyapunov_db = 0; +%$ options_.lyapunov_srs = 1; +%$ try +%$ Pstar3_small = lyapunov_solver(T_small,R_small,Q_small,options_); +%$ Pstar3_large = lyapunov_solver(T_large,R_large,Q_large,options_); +%$ t(3) = 1; +%$ catch +%$ t(3) = 0; +%$ end +%$ else +%$ t(3) = 1; +%$ end +%$ +%$ % Standard +%$ options_.lyapunov_srs = 0; +%$ try +%$ Pstar4_small = lyapunov_solver(T_small,R_small,Q_small,options_); +%$ Pstar4_large = lyapunov_solver(T_large,R_large,Q_large,options_); +%$ t(4) = 1; +%$ catch +%$ t(4) = 0; +%$ end +%$ +%$ % Test the results. +%$ +%$ if max(max(abs(Pstar1_small-Pstar2_small)))>1e-8 +%$ t(5) = 0; +%$ else +%$ t(5) = 1; +%$ end +%$ +%$ if (isoctave && user_has_octave_forge_package('control')) || (~isoctave && user_has_matlab_license('control_toolbox')) +%$ if max(max(abs(Pstar1_small-Pstar3_small)))>1e-8 +%$ t(6) = 0; +%$ else +%$ t(6) = 1; +%$ end +%$ else +%$ t(6) = 1; +%$ end +%$ +%$ if max(max(abs(Pstar1_small-Pstar4_small)))>1e-8 +%$ t(7) = 0; +%$ else +%$ t(7) = 1; +%$ end +%$ +%$ if max(max(abs(Pstar1_large-Pstar2_large)))>1e-8 +%$ t(8) = 0; +%$ else +%$ t(8) = 1; +%$ end +%$ +%$ if (isoctave && user_has_octave_forge_package('control')) || (~isoctave && user_has_matlab_license('control_toolbox')) +%$ if max(max(abs(Pstar1_large-Pstar3_large)))>1e-8 +%$ t(9) = 0; +%$ else +%$ t(9) = 1; +%$ end +%$ else +%$ t(9) = 1; +%$ end +%$ +%$ if max(max(abs(Pstar1_large-Pstar4_large)))>1e-8 +%$ t(10) = 0; +%$ else +%$ t(10) = 1; +%$ end +%$ +%$ T = all(t); +%@eof:1 \ No newline at end of file diff --git a/matlab/lyapunov_symm.m b/matlab/lyapunov_symm.m index 81208c358..694415d71 100644 --- a/matlab/lyapunov_symm.m +++ b/matlab/lyapunov_symm.m @@ -1,4 +1,4 @@ -function [x,u] = lyapunov_symm(a,b,lyapunov_fixed_point_tol,qz_criterium,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,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. % @@ -20,12 +20,12 @@ function [x,u] = lyapunov_symm(a,b,lyapunov_fixed_point_tol,qz_criterium,lyapuno % % ALGORITHM % Uses reordered Schur decomposition (Bartels-Stewart algorithm) -% [method<3] or a fixed point algorithm (method==4) +% [method<3] or a fixed point algorithm (method==3) % % SPECIAL REQUIREMENTS % None -% Copyright (C) 2006-2014 Dynare Team +% Copyright (C) 2006-2016 Dynare Team % % This file is part of Dynare. % @@ -46,7 +46,7 @@ if nargin<6 || isempty(method) method = 0; end -if nargin<8 +if nargin<7 debug = 0; end @@ -90,22 +90,7 @@ if method == 3 return; end; end; -elseif method == 4 - % works only with Matlab System Control toolbox or octave the control package, - if isoctave - if ~user_has_octave_forge_package('control') - error('lyapunov=square_root_solver not available; you must install the control package from Octave Forge') - end - else - if ~user_has_matlab_license('control_toolbox') - error('lyapunov=square_root_solver not available; you must install the control system toolbox') - end - end - chol_b = R*chol(b,'lower'); - Rx = dlyapchol(a,chol_b); - x = Rx' * Rx; - return; -end; +end if method persistent U T k n @@ -178,117 +163,4 @@ if i == 1 x(1,1) = (B(1,1)+c)/(1-T(1,1)*T(1,1)); end x = U(:,k+1:end)*x*U(:,k+1:end)'; -u = U(:,1:k); - -%@test:1 -%$ t = NaN(10,1); -%$ lyapunov_complex_threshold = 1e-15; -%$ qz_zero_threshold = 1e-6; -%$ qz_criterium=1-qz_zero_threshold; -%$ lyapunov_fixed_point_tol = 1e-10; -%$ lyapunov_doubling_tol = 1e-16; -%$ -%$ n_small=8; -%$ m_small=10; -%$ T_small=randn(n_small,n_small); -%$ T_small=0.99*T_small/max(abs(eigs(T_small))); -%$ tmp2=randn(m_small,m_small); -%$ Q_small=tmp2*tmp2'; -%$ R_small=randn(n_small,m_small); -%$ -%$ n_large=9; -%$ m_large=11; -%$ T_large=randn(n_large,n_large); -%$ T_large=0.99*T_large/max(abs(eigs(T_large))); -%$ tmp2=randn(m_large,m_large); -%$ Q_large=tmp2*tmp2'; -%$ R_large=randn(n_large,m_large); -%$ -%$ % DynareOptions.lyapunov_fp == 1 -%$ try -%$ 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; -%$ end -%$ -%$ % Dynareoptions.lyapunov_db == 1 -%$ try -%$ Pstar2_small = disclyap_fast(T_small,R_small*Q_small*R_small',lyapunov_doubling_tol); -%$ Pstar2_large = disclyap_fast(T_large,R_large*Q_large*R_large',lyapunov_doubling_tol); -%$ t(2) = 1; -%$ catch -%$ t(2) = 0; -%$ end -%$ -%$ % 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,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; -%$ end -%$ else -%$ t(3) = 1; -%$ end -%$ -%$ % Standard -%$ try -%$ 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; -%$ end -%$ -%$ % Test the results. -%$ -%$ if max(max(abs(Pstar1_small-Pstar2_small)))>1e-8 -%$ t(5) = 0; -%$ else -%$ t(5) = 1; -%$ end -%$ -%$ if (isoctave && user_has_octave_forge_package('control')) || (~isoctave && user_has_matlab_license('control_toolbox')) -%$ if max(max(abs(Pstar1_small-Pstar3_small)))>1e-8 -%$ t(6) = 0; -%$ else -%$ t(6) = 1; -%$ end -%$ else -%$ t(6) = 1; -%$ end -%$ -%$ if max(max(abs(Pstar1_small-Pstar4_small)))>1e-8 -%$ t(7) = 0; -%$ else -%$ t(7) = 1; -%$ end -%$ -%$ if max(max(abs(Pstar1_large-Pstar2_large)))>1e-8 -%$ t(8) = 0; -%$ else -%$ t(8) = 1; -%$ end -%$ -%$ if (isoctave && user_has_octave_forge_package('control')) || (~isoctave && user_has_matlab_license('control_toolbox')) -%$ if max(max(abs(Pstar1_large-Pstar3_large)))>1e-8 -%$ t(9) = 0; -%$ else -%$ t(9) = 1; -%$ end -%$ else -%$ t(9) = 1; -%$ end -%$ -%$ if max(max(abs(Pstar1_large-Pstar4_large)))>1e-8 -%$ t(10) = 0; -%$ else -%$ t(10) = 1; -%$ end -%$ -%$ T = all(t); -%@eof:1 \ No newline at end of file +u = U(:,1:k); \ No newline at end of file diff --git a/matlab/non_linear_dsge_likelihood.m b/matlab/non_linear_dsge_likelihood.m index 34ed485c4..f8ac92e23 100644 --- a/matlab/non_linear_dsge_likelihood.m +++ b/matlab/non_linear_dsge_likelihood.m @@ -306,7 +306,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,:)',DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold,[],[],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 461efc26f..3e101a7a6 100644 --- a/matlab/partial_information/disclyap_fast.m +++ b/matlab/partial_information/disclyap_fast.m @@ -1,25 +1,25 @@ -function [X,exitflag]=disclyap_fast(G,V,tol,ch) +function [X,exitflag]=disclyap_fast(G,V,tol,check_flag) % 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 +% - G [double] first input matrix +% - V [double] second input matrix +% - tol [scalar] tolerance criterion +% - check_flag empty of non-empty if non-empty: check positive-definiteness % Outputs: -% - X [double] solution matrix -% - exitflag [scalar] 0 if solution is found, 1 otherwise +% - 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 % -% If ch is defined then the code will check if the resulting X +% If check_flag is defined then the code will check if the resulting X % is positive definite and generate an error message if it is not % % Joe Pearlman and Alejandro Justiniano % 3/5/2005 -% Copyright (C) 2010-2015 Dynare Team +% Copyright (C) 2010-2016 Dynare Team % % This file is part of Dynare. % @@ -36,16 +36,13 @@ function [X,exitflag]=disclyap_fast(G,V,tol,ch) % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . -if nargin <= 3 || isempty( ch ) == 1 +if nargin <= 3 || isempty( check_flag ) == 1 flag_ch = 0; else flag_ch = 1; end -s=size(G,1); exitflag=0; -%tol = 1e-16; - P0=V; A0=G; diff --git a/matlab/th_autocovariances.m b/matlab/th_autocovariances.m index 1102c2f90..966ece834 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_.lyapunov_fixed_point_tol,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 && ~options_.bandpass.indicator b1 = b1*cs; b2(:,exo_names_orig_ord) = ghu(iky,:); b2 = b2*cs; - vx = lyapunov_symm(A,b1*b1',options_.lyapunov_fixed_point_tol,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_.lyapunov_fixed_point_tol,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 d938c74d0..929193ca4 100644 --- a/matlab/thet2tau.m +++ b/matlab/thet2tau.m @@ -55,7 +55,7 @@ elseif flagmoments==-1 tau=[oo_.dr.ys(oo_.dr.order_var); g1(:)]; else - GAM = lyapunov_symm(A,B*M_.Sigma_e*B',options_.lyapunov_fixed_point_tol,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,