disc_riccati_fast

 function Z=disc_riccati_fast(F,D,R,H,ch)
 
 Solves discrete Riccati Equation: 
 Z=FZF' - FZD'inv(DZD'+R)DZF' + H
 Using the Doubling Algorithm 

 George Perendia: based on the doubling algorithm for Riccati   
 and the disclyap_fast function provided by Prof. Joe Pearlman 
 V.1 19/5/2006
 V.2 22/10/06
 =================================================================

0001 function Z=disc_riccati_fast(F,D,R,H,ch)
0002 % function Z=disc_riccati_fast(F,D,R,H,ch)
0003 %
0004 % Solves discrete Riccati Equation:
0005 % Z=FZF' - FZD'inv(DZD'+R)DZF' + H
0006 % Using the Doubling Algorithm
0007 %
0008 % George Perendia: based on the doubling algorithm for Riccati
0009 % and the disclyap_fast function provided by Prof. Joe Pearlman
0010 % V.1 19/5/2006
0011 % V.2 22/10/06
0012 % =================================================================
0013 
0014 % Copyright (C) 2006-2011 Dynare Team
0015 %
0016 % This file is part of Dynare.
0017 %
0018 % Dynare is free software: you can redistribute it and/or modify
0019 % it under the terms of the GNU General Public License as published by
0020 % the Free Software Foundation, either version 3 of the License, or
0021 % (at your option) any later version.
0022 %
0023 % Dynare is distributed in the hope that it will be useful,
0024 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0025 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0026 % GNU General Public License for more details.
0027 %
0028 % You should have received a copy of the GNU General Public License
0029 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0030 
0031 if nargin == 4 || isempty( ch ) == 1 
0032     flag_ch = 0; 
0033 else 
0034     flag_ch = 1; 
0035 end 
0036 
0037 
0038 % intialisation
0039 tol = 1e-10; % iteration convergence threshold
0040 P0=H; 
0041 X0=F;
0042 if ~any(R) % i.e. ==0
0043     warning('Dangerously evading inversion of zero matrix!');
0044     Y0=zeros(size(R));
0045 else
0046     Y0=D'*inv(R)*D;
0047 end
0048 POYO=P0*Y0;
0049 I=eye(size(POYO));
0050 clear POYO;
0051 
0052 % iterate Riccati equation solution
0053 matd=1; 
0054 count=0;
0055 while matd > tol && count < 100
0056     INVPY=inv(I+P0*Y0);
0057     P1=X0*INVPY*P0*X0'+ P0; 
0058     Y1=X0'*Y0*INVPY*X0+ Y0; 
0059     X1=X0*INVPY*X0; 
0060     matd=sum( sum(abs( P1 - P0 ))); 
0061     %    P0=(P1+P1')/2
0062     P0=P1; 
0063     X0=X1;
0064     Y0=Y1;
0065     count=count+1;
0066     %    matd;
0067 end 
0068 
0069 Z=(P0+P0')/2;
0070 %Z=P0
0071 % check if the convergence took place
0072 if count==100
0073     matd
0074     error('Riccati not converged fast enough!');
0075     %    error.identifier='Riccati not converged!'
0076     %    error
0077 end
0078 %if count >5
0079 %    disp('Riccati count= ');
0080 %    count
0081 %end
0082 
0083 clear X0 X1 Y0 Y1 P1 I INVPY; 
0084 
0085 % Check that X is positive definite
0086 if flag_ch==1 
0087     [C,p]=chol(Z); 
0088     if p ~= 0 
0089         error('Z is not positive definite')
0090     end 
0091 end