1function P=lyapunov_solver(T,R,Q,DynareOptions) % --*-- Unitary tests --*-- 2% function P=lyapunov_solver(T,R,Q,DynareOptions) 3% Solves the Lyapunov equation P-T*P*T' = R*Q*R' arising in a state-space 4% system, where P is the variance of the states 5% 6% Inputs 7% T [double] n*n matrix. 8% R [double] n*m matrix. 9% Q [double] m*m matrix. 10% DynareOptions [structure] Dynare options 11% 12% Outputs 13% P [double] n*n matrix. 14% 15% Algorithms 16% Default, if none of the other algorithms is selected: 17% Reordered Schur decomposition (Bartels-Stewart algorithm) 18% DynareOptions.lyapunov_fp == true 19% iteration-based fixed point algorithm 20% DynareOptions.lyapunov_db == true 21% doubling algorithm 22% DynareOptions.lyapunov_srs == true 23% Square-root solver for discrete-time Lyapunov equations (requires Matlab System Control toolbox 24% or Octave control package) 25 26% Copyright (C) 2016-2017 Dynare Team 27% 28% This file is part of Dynare. 29% 30% Dynare is free software: you can redistribute it and/or modify 31% it under the terms of the GNU General Public License as published by 32% the Free Software Foundation, either version 3 of the License, or 33% (at your option) any later version. 34% 35% Dynare is distributed in the hope that it will be useful, 36% but WITHOUT ANY WARRANTY; without even the implied warranty of 37% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 38% GNU General Public License for more details. 39% 40% You should have received a copy of the GNU General Public License 41% along with Dynare. If not, see <http://www.gnu.org/licenses/>. 42 43if DynareOptions.lyapunov_fp 44 P = lyapunov_symm(T,R*Q*R',DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, 3, DynareOptions.debug); 45elseif DynareOptions.lyapunov_db 46 [P, errorflag] = disclyap_fast(T,R*Q*R',DynareOptions.lyapunov_doubling_tol); 47 if errorflag %use Schur-based method 48 P = lyapunov_symm(T,R*Q*R',DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, [], DynareOptions.debug); 49 end 50elseif DynareOptions.lyapunov_srs 51 % works only with Matlab System Control toolbox or Octave control package, 52 if isoctave 53 if ~user_has_octave_forge_package('control') 54 error('lyapunov=square_root_solver not available; you must install the control package from Octave Forge') 55 end 56 else 57 if ~user_has_matlab_license('control_toolbox') 58 error('lyapunov=square_root_solver not available; you must install the control system toolbox') 59 end 60 end 61 chol_Q = R*chol(Q,'lower'); 62 R_P = dlyapchol(T,chol_Q); 63 P = R_P' * R_P; 64else 65 P = lyapunov_symm(T,R*Q*R',DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, [], DynareOptions.debug); 66end 67 68%@test:1 69%$ t = NaN(10,1); 70%$ options_.lyapunov_complex_threshold = 1e-15; 71%$ options_.qz_zero_threshold = 1e-6; 72%$ options_.qz_criterium=1-options_.qz_zero_threshold; 73%$ options_.lyapunov_fixed_point_tol = 1e-10; 74%$ options_.lyapunov_doubling_tol = 1e-16; 75%$ options_.debug=false; 76%$ 77%$ n_small=8; 78%$ m_small=10; 79%$ T_small=randn(n_small,n_small); 80%$ T_small=0.99*T_small/max(abs(eigs(T_small))); 81%$ tmp2=randn(m_small,m_small); 82%$ Q_small=tmp2*tmp2'; 83%$ R_small=randn(n_small,m_small); 84%$ 85%$ n_large=9; 86%$ m_large=11; 87%$ T_large=randn(n_large,n_large); 88%$ T_large=0.99*T_large/max(abs(eigs(T_large))); 89%$ tmp2=randn(m_large,m_large); 90%$ Q_large=tmp2*tmp2'; 91%$ R_large=randn(n_large,m_large); 92%$ 93%$ % DynareOptions.lyapunov_fp == 1 94%$ options_.lyapunov_fp = true; 95%$ try 96%$ Pstar1_small = lyapunov_solver(T_small,R_small,Q_small,options_); 97%$ Pstar1_large = lyapunov_solver(T_large,R_large,Q_large,options_); 98%$ t(1) = 1; 99%$ catch 100%$ t(1) = 0; 101%$ end 102%$ 103%$ % Dynareoptions.lyapunov_db == 1 104%$ options_.lyapunov_fp = false; 105%$ options_.lyapunov_db = true; 106%$ try 107%$ Pstar2_small = lyapunov_solver(T_small,R_small,Q_small,options_); 108%$ Pstar2_large = lyapunov_solver(T_large,R_large,Q_large,options_); 109%$ t(2) = 1; 110%$ catch 111%$ t(2) = 0; 112%$ end 113%$ 114%$ % Dynareoptions.lyapunov_srs == 1 115%$ if (isoctave && user_has_octave_forge_package('control')) || (~isoctave && user_has_matlab_license('control_toolbox')) 116%$ options_.lyapunov_db = false; 117%$ options_.lyapunov_srs = true; 118%$ try 119%$ Pstar3_small = lyapunov_solver(T_small,R_small,Q_small,options_); 120%$ Pstar3_large = lyapunov_solver(T_large,R_large,Q_large,options_); 121%$ t(3) = 1; 122%$ catch 123%$ t(3) = 0; 124%$ end 125%$ else 126%$ t(3) = 1; 127%$ end 128%$ 129%$ % Standard 130%$ options_.lyapunov_srs = false; 131%$ try 132%$ Pstar4_small = lyapunov_solver(T_small,R_small,Q_small,options_); 133%$ Pstar4_large = lyapunov_solver(T_large,R_large,Q_large,options_); 134%$ t(4) = 1; 135%$ catch 136%$ t(4) = 0; 137%$ end 138%$ 139%$ % Test the results. 140%$ 141%$ if max(max(abs(Pstar1_small-Pstar2_small)))>1e-8 142%$ t(5) = 0; 143%$ else 144%$ t(5) = 1; 145%$ end 146%$ 147%$ if (isoctave && user_has_octave_forge_package('control')) || (~isoctave && user_has_matlab_license('control_toolbox')) 148%$ if max(max(abs(Pstar1_small-Pstar3_small)))>1e-8 149%$ t(6) = 0; 150%$ else 151%$ t(6) = 1; 152%$ end 153%$ else 154%$ t(6) = 1; 155%$ end 156%$ 157%$ if max(max(abs(Pstar1_small-Pstar4_small)))>1e-8 158%$ t(7) = 0; 159%$ else 160%$ t(7) = 1; 161%$ end 162%$ 163%$ if max(max(abs(Pstar1_large-Pstar2_large)))>1e-8 164%$ t(8) = 0; 165%$ else 166%$ t(8) = 1; 167%$ end 168%$ 169%$ if (isoctave && user_has_octave_forge_package('control')) || (~isoctave && user_has_matlab_license('control_toolbox')) 170%$ if max(max(abs(Pstar1_large-Pstar3_large)))>1e-8 171%$ t(9) = 0; 172%$ else 173%$ t(9) = 1; 174%$ end 175%$ else 176%$ t(9) = 1; 177%$ end 178%$ 179%$ if max(max(abs(Pstar1_large-Pstar4_large)))>1e-8 180%$ t(10) = 0; 181%$ else 182%$ t(10) = 1; 183%$ end 184%$ 185%$ T = all(t); 186%@eof:1 187