diff --git a/matlab/partial_information/PCL_Part_info_irf.m b/matlab/partial_information/PCL_Part_info_irf.m
index 94949c444..4eebceea4 100644
--- a/matlab/partial_information/PCL_Part_info_irf.m
+++ b/matlab/partial_information/PCL_Part_info_irf.m
@@ -1,162 +1,162 @@
-function [irfmat,irfst]=PCL_Part_info_irf( H, varobs, M_, dr, irfpers,ii)
-% sets up parameters and calls part-info kalman filter
-% developed by G Perendia, July 2006 for implementation from notes by Prof. Joe Pearlman to
-% suit partial information RE solution in accordance with, and based on, the
-% Pearlman, Currie and Levine 1986 solution.
-% 22/10/06 - Version 2 for new Riccati with 4 params instead 5
-
-% Copyright (C) 2001-20010 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 .
-
-
-% Recall that the state space is given by the
-% predetermined variables s(t-1), x(t-1)
-% and the jump variables x(t).
-% The jump variables have dimension NETA
-
-
- [junk,OBS] = ismember(varobs,M_.endo_names,'rows');
-
- G1=dr.PI_ghx;
- impact=dr.PI_ghu;
- nmat=dr.PI_nmat;
- CC=dr.PI_CC;
- NX=M_.exo_nbr; % no of exogenous varexo shock variables.
- FL_RANK=dr.PI_FL_RANK;
- NY=M_.endo_nbr;
- if isempty(OBS)
- NOBS=NY;
- LL=eye(NY,NY);
- else %and if no obsevations specify OBS=[0] but this is not going to work properly
- NOBS=length(OBS);
- LL=zeros(NOBS,NY);
- for i=1:NOBS
- LL(i,OBS(i))=1;
- end
- end
-
- if exist( 'irfpers')==1
- if ~isempty(irfpers)
- if irfpers<=0, irfpers=20, end;
- else
- irfpers=20;
- end
- else
- irfpers=20;
- end
-
- ss=size(G1,1);
- pd=ss-size(nmat,1);
- SDX=M_.Sigma_e^0.5; % =SD,not V-COV, of Exog shocks or M_.Sigma_e^0.5 num_exog x num_exog matrix
- if isempty(H)
- H=M_.H;
- end
- VV=H; % V-COV of observation errors.
- MM=impact*SDX; % R*(Q^0.5) in standard KF notation
- % observation vector indices
- % mapping to endogenous variables.
-
- L1=LL*dr.PI_TT1;
- L2=LL*dr.PI_TT2;
-
- MM1=MM(1:ss-FL_RANK,:);
- U11=MM1*MM1';
- % SDX
-
- U22=0;
- % determine K1 and K2 observation mapping matrices
- % This uses the fact that measurements are given by L1*s(t)+L2*x(t)
- % and s(t) is expressed in the dynamics as
- % H1*eps(t)+G11*s(t-1)+G12*x(t-1)+G13*x(t).
- % Thus the observations o(t) can be written in the form
- % o(t)=K1*[eps(t)' s(t-1)' x(t-1)']' + K2*x(t) where
- % K1=[L1*H1 L1*G11 L1*G12] K2=L1*G13+L2
-
- G12=G1(NX+1:ss-2*FL_RANK,:);
- KK1=L1*G12;
- K1=KK1(:,1:ss-FL_RANK);
- K2=KK1(:,ss-FL_RANK+1:ss)+L2;
-
- %pre calculate time-invariant factors
- A11=G1(1:pd,1:pd);
- A22=G1(pd+1:end, pd+1:end);
- A12=G1(1:pd, pd+1:end);
- A21=G1(pd+1:end,1:pd);
- Lambda= nmat*A12+A22;
- %A11_A12Nmat= A11-A12*nmat % test
- I_L=inv(Lambda);
- BB=A12*inv(A22);
- FF=K2*inv(A22);
- QQ=BB*U22*BB' + U11;
- UFT=U22*FF';
- % kf_param structure:
- AA=A11-BB*A21;
- CCCC=A11-A12*nmat; % F in new notation
- DD=K1-FF*A21; % H in new notation
- EE=K1-K2*nmat;
- RR=FF*UFT+VV;
- if ~any(RR)
- % if zero add some dummy measurement err. variance-covariances
- % with diagonals 0.000001. This would not be needed if we used
- % the slow solver, or the generalised eigenvalue approach,
- % but these are both slower.
- RR=eye(size(RR,1))*1.0e-6;
- end
- SS=BB*UFT;
- VKLUFT=VV+K2*I_L*UFT;
- ALUFT=A12*I_L*UFT;
- FULKV=FF*U22*I_L'*K2'+VV;
- FUBT=FF*U22*BB';
- nmat=nmat;
- % initialise pshat
- AQDS=AA*QQ*DD'+SS;
- DQDR=DD*QQ*DD'+RR;
- I_DQDR=inv(DQDR);
- AQDQ=AQDS*I_DQDR;
- ff=AA-AQDQ*DD;
- hh=AA*QQ*AA'-AQDQ*AQDS';%*(DD*QQ*AA'+SS');
- rr=DD*QQ*DD'+RR;
- ZSIG0=disc_riccati_fast(ff,DD,rr,hh);
- PP=ZSIG0 +QQ;
-
- exo_names=M_.exo_names(M_.exo_names_orig_ord,:);
-
- DPDR=DD*PP*DD'+RR;
- I_DPDR=inv(DPDR);
- PDIDPDRD=PP*DD'*I_DPDR*DD;
- %GG=[CCCC (AA-CCCC)*(eye(ss-FL_RANK)-PP*DD'*I_DQDR*DD); zeros(ss-FL_RANK) AA*(eye(ss-FL_RANK)-PP*DD'*I_DQDR*DD)];
- GG=[CCCC (AA-CCCC)*(eye(ss-FL_RANK)-PDIDPDRD); zeros(ss-FL_RANK) AA*(eye(ss-FL_RANK)-PDIDPDRD)];
- imp=[impact(1:ss-FL_RANK,:); impact(1:ss-FL_RANK,:)];
- % Calculate IRFs of observable series
- %The extra term in leads to
- %LL0=[EE (H-EE)(I-PH^T(HPH^T+V)^{-1}H)].
- %Then in the case of observing all variables without noise (V=0), this
- % leads to LL0=[EE 0].
- I_PD=(eye(ss-FL_RANK)-PDIDPDRD);
- LL0=[ EE (DD-EE)*I_PD];
- %OVV = [ zeros( size(dr.PI_TT1,1), NX ) dr.PI_TT1 dr.PI_TT2];
- VV = [ dr.PI_TT1 dr.PI_TT2];
- stderr=diag(M_.Sigma_e^0.5);
- irfmat=zeros(size(dr.PI_TT1 ,1),irfpers);
- irfst=zeros(size(GG,1),irfpers);
- irfst(:,1)=stderr(ii)*imp(:,ii);
- for jj=2:irfpers;
- irfst(:,jj)=GG*irfst(:,jj-1); % xjj=f irfstjj-2
- irfmat(:,jj-1)=VV*irfst(NX+1:ss-FL_RANK,jj);
- %irfmat(:,jj)=LL0*irfst(:,jj);
- end
- save ([M_.fname '_PCL_PtInfoIRFs_' num2str(ii) '_' deblank(exo_names(ii,:))], 'irfmat','irfst');
+function [irfmat,irfst]=PCL_Part_info_irf( H, varobs, M_, dr, irfpers,ii)
+% sets up parameters and calls part-info kalman filter
+% developed by G Perendia, July 2006 for implementation from notes by Prof. Joe Pearlman to
+% suit partial information RE solution in accordance with, and based on, the
+% Pearlman, Currie and Levine 1986 solution.
+% 22/10/06 - Version 2 for new Riccati with 4 params instead 5
+
+% Copyright (C) 2001-20010 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 .
+
+
+% Recall that the state space is given by the
+% predetermined variables s(t-1), x(t-1)
+% and the jump variables x(t).
+% The jump variables have dimension NETA
+
+
+ [junk,OBS] = ismember(varobs,M_.endo_names,'rows');
+
+ G1=dr.PI_ghx;
+ impact=dr.PI_ghu;
+ nmat=dr.PI_nmat;
+ CC=dr.PI_CC;
+ NX=M_.exo_nbr; % no of exogenous varexo shock variables.
+ FL_RANK=dr.PI_FL_RANK;
+ NY=M_.endo_nbr;
+ if isempty(OBS)
+ NOBS=NY;
+ LL=eye(NY,NY);
+ else %and if no obsevations specify OBS=[0] but this is not going to work properly
+ NOBS=length(OBS);
+ LL=zeros(NOBS,NY);
+ for i=1:NOBS
+ LL(i,OBS(i))=1;
+ end
+ end
+
+ if exist( 'irfpers')==1
+ if ~isempty(irfpers)
+ if irfpers<=0, irfpers=20, end;
+ else
+ irfpers=20;
+ end
+ else
+ irfpers=20;
+ end
+
+ ss=size(G1,1);
+ pd=ss-size(nmat,1);
+ SDX=M_.Sigma_e^0.5; % =SD,not V-COV, of Exog shocks or M_.Sigma_e^0.5 num_exog x num_exog matrix
+ if isempty(H)
+ H=M_.H;
+ end
+ VV=H; % V-COV of observation errors.
+ MM=impact*SDX; % R*(Q^0.5) in standard KF notation
+ % observation vector indices
+ % mapping to endogenous variables.
+
+ L1=LL*dr.PI_TT1;
+ L2=LL*dr.PI_TT2;
+
+ MM1=MM(1:ss-FL_RANK,:);
+ U11=MM1*MM1';
+ % SDX
+
+ U22=0;
+ % determine K1 and K2 observation mapping matrices
+ % This uses the fact that measurements are given by L1*s(t)+L2*x(t)
+ % and s(t) is expressed in the dynamics as
+ % H1*eps(t)+G11*s(t-1)+G12*x(t-1)+G13*x(t).
+ % Thus the observations o(t) can be written in the form
+ % o(t)=K1*[eps(t)' s(t-1)' x(t-1)']' + K2*x(t) where
+ % K1=[L1*H1 L1*G11 L1*G12] K2=L1*G13+L2
+
+ G12=G1(NX+1:ss-2*FL_RANK,:);
+ KK1=L1*G12;
+ K1=KK1(:,1:ss-FL_RANK);
+ K2=KK1(:,ss-FL_RANK+1:ss)+L2;
+
+ %pre calculate time-invariant factors
+ A11=G1(1:pd,1:pd);
+ A22=G1(pd+1:end, pd+1:end);
+ A12=G1(1:pd, pd+1:end);
+ A21=G1(pd+1:end,1:pd);
+ Lambda= nmat*A12+A22;
+ %A11_A12Nmat= A11-A12*nmat % test
+ I_L=inv(Lambda);
+ BB=A12*inv(A22);
+ FF=K2*inv(A22);
+ QQ=BB*U22*BB' + U11;
+ UFT=U22*FF';
+ % kf_param structure:
+ AA=A11-BB*A21;
+ CCCC=A11-A12*nmat; % F in new notation
+ DD=K1-FF*A21; % H in new notation
+ EE=K1-K2*nmat;
+ RR=FF*UFT+VV;
+ if ~any(RR)
+ % if zero add some dummy measurement err. variance-covariances
+ % with diagonals 0.000001. This would not be needed if we used
+ % the slow solver, or the generalised eigenvalue approach,
+ % but these are both slower.
+ RR=eye(size(RR,1))*1.0e-6;
+ end
+ SS=BB*UFT;
+ VKLUFT=VV+K2*I_L*UFT;
+ ALUFT=A12*I_L*UFT;
+ FULKV=FF*U22*I_L'*K2'+VV;
+ FUBT=FF*U22*BB';
+ nmat=nmat;
+ % initialise pshat
+ AQDS=AA*QQ*DD'+SS;
+ DQDR=DD*QQ*DD'+RR;
+ I_DQDR=inv(DQDR);
+ AQDQ=AQDS*I_DQDR;
+ ff=AA-AQDQ*DD;
+ hh=AA*QQ*AA'-AQDQ*AQDS';%*(DD*QQ*AA'+SS');
+ rr=DD*QQ*DD'+RR;
+ ZSIG0=disc_riccati_fast(ff,DD,rr,hh);
+ PP=ZSIG0 +QQ;
+
+ exo_names=M_.exo_names(M_.exo_names_orig_ord,:);
+
+ DPDR=DD*PP*DD'+RR;
+ I_DPDR=inv(DPDR);
+ PDIDPDRD=PP*DD'*I_DPDR*DD;
+ %GG=[CCCC (AA-CCCC)*(eye(ss-FL_RANK)-PP*DD'*I_DQDR*DD); zeros(ss-FL_RANK) AA*(eye(ss-FL_RANK)-PP*DD'*I_DQDR*DD)];
+ GG=[CCCC (AA-CCCC)*(eye(ss-FL_RANK)-PDIDPDRD); zeros(ss-FL_RANK) AA*(eye(ss-FL_RANK)-PDIDPDRD)];
+ imp=[impact(1:ss-FL_RANK,:); impact(1:ss-FL_RANK,:)];
+ % Calculate IRFs of observable series
+ %The extra term in leads to
+ %LL0=[EE (H-EE)(I-PH^T(HPH^T+V)^{-1}H)].
+ %Then in the case of observing all variables without noise (V=0), this
+ % leads to LL0=[EE 0].
+ I_PD=(eye(ss-FL_RANK)-PDIDPDRD);
+ LL0=[ EE (DD-EE)*I_PD];
+ %OVV = [ zeros( size(dr.PI_TT1,1), NX ) dr.PI_TT1 dr.PI_TT2];
+ VV = [ dr.PI_TT1 dr.PI_TT2];
+ stderr=diag(M_.Sigma_e^0.5);
+ irfmat=zeros(size(dr.PI_TT1 ,1),irfpers);
+ irfst=zeros(size(GG,1),irfpers);
+ irfst(:,1)=stderr(ii)*imp(:,ii);
+ for jj=2:irfpers;
+ irfst(:,jj)=GG*irfst(:,jj-1); % xjj=f irfstjj-2
+ irfmat(:,jj-1)=VV*irfst(NX+1:ss-FL_RANK,jj);
+ %irfmat(:,jj)=LL0*irfst(:,jj);
+ end
+ save ([M_.fname '_PCL_PtInfoIRFs_' num2str(ii) '_' deblank(exo_names(ii,:))], 'irfmat','irfst');
diff --git a/matlab/partial_information/PCL_Part_info_moments.m b/matlab/partial_information/PCL_Part_info_moments.m
index 66650d49c..84555f2b5 100644
--- a/matlab/partial_information/PCL_Part_info_moments.m
+++ b/matlab/partial_information/PCL_Part_info_moments.m
@@ -1,200 +1,200 @@
-function [irfmat,irfst]=PCL_Part_info_moments( H, varobs, dr,ivar)
-% sets up parameters and calls part-info kalman filter
-% developed by G Perendia, July 2006 for implementation from notes by Prof. Joe Pearlman to
-% suit partial information RE solution in accordance with, and based on, the
-% Pearlman, Currie and Levine 1986 solution.
-% 22/10/06 - Version 2 for new Riccati with 4 params instead 5
-
-% Copyright (C) 2001-20010 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 .
-
- % Recall that the state space is given by the
- % predetermined variables s(t-1), x(t-1)
- % and the jump variables x(t).
- % The jump variables have dimension NETA
-
- global M_ options_ oo_
- warning_old_state = warning;
- warning off
-
- [junk,OBS] = ismember(varobs,M_.endo_names,'rows');
-
- G1=dr.PI_ghx;
- impact=dr.PI_ghu;
- nmat=dr.PI_nmat;
- CC=dr.PI_CC;
- NX=M_.exo_nbr; % no of exogenous varexo shock variables.
-% NETA=dr.nfwrd+dr.nboth; % total no of exp. errors set to no of forward looking equations
- FL_RANK=dr.PI_FL_RANK;
- NY=M_.endo_nbr;
- if isempty(OBS)
- NOBS=NY;
- LL=eye(NY,NY);
- else %and if no obsevations specify OBS=[0] but this is not going to work properly
- NOBS=length(OBS);
- LL=zeros(NOBS,NY);
- for i=1:NOBS
- LL(i,OBS(i))=1;
- end
- end
-
- if exist( 'irfpers')==1
- if ~isempty(irfpers)
- if irfpers<=0, irfpers=20, end;
- else
- irfpers=20;
- end
- else
- irfpers=20;
- end
-
- ss=size(G1,1);
-
- pd=ss-size(nmat,1);
- SDX=M_.Sigma_e^0.5; % =SD,not V-COV, of Exog shocks or M_.Sigma_e^0.5 num_exog x num_exog matrix
- if isempty(H)
- H=M_.H;
- end
- VV=H; % V-COV of observation errors.
- MM=impact*SDX; % R*(Q^0.5) in standard KF notation
- % observation vector indices
- % mapping to endogenous variables.
-
- L1=LL*dr.PI_TT1;
- L2=LL*dr.PI_TT2;
-
- MM1=MM(1:ss-FL_RANK,:);
- U11=MM1*MM1';
- % SDX
-
- U22=0;
- % determine K1 and K2 observation mapping matrices
- % This uses the fact that measurements are given by L1*s(t)+L2*x(t)
- % and s(t) is expressed in the dynamics as
- % H1*eps(t)+G11*s(t-1)+G12*x(t-1)+G13*x(t).
- % Thus the observations o(t) can be written in the form
- % o(t)=K1*[eps(t)' s(t-1)' x(t-1)']' + K2*x(t) where
- % K1=[L1*H1 L1*G11 L1*G12] K2=L1*G13+L2
-
- G12=G1(NX+1:ss-2*FL_RANK,:);
- KK1=L1*G12;
- K1=KK1(:,1:ss-FL_RANK);
- K2=KK1(:,ss-FL_RANK+1:ss)+L2;
-
- %pre calculate time-invariant factors
- A11=G1(1:pd,1:pd);
- A22=G1(pd+1:end, pd+1:end);
- A12=G1(1:pd, pd+1:end);
- A21=G1(pd+1:end,1:pd);
- Lambda= nmat*A12+A22;
- %A11_A12Nmat= A11-A12*nmat % test
- I_L=inv(Lambda);
- BB=A12*inv(A22);
- FF=K2*inv(A22);
- QQ=BB*U22*BB' + U11;
- UFT=U22*FF';
- % kf_param structure:
- AA=A11-BB*A21;
- CCCC=A11-A12*nmat; % F in new notation
- DD=K1-FF*A21; % H in new notation
- EE=K1-K2*nmat;
- RR=FF*UFT+VV;
- if ~any(RR)
- % if zero add some dummy measurement err. variance-covariances
- % with diagonals 0.000001. This would not be needed if we used
- % the slow solver, or the generalised eigenvalue approach,
- % but these are both slower.
- RR=eye(size(RR,1))*1.0e-6;
- end
- SS=BB*UFT;
- VKLUFT=VV+K2*I_L*UFT;
- ALUFT=A12*I_L*UFT;
- FULKV=FF*U22*I_L'*K2'+VV;
- FUBT=FF*U22*BB';
- nmat=nmat;
- % initialise pshat
- AQDS=AA*QQ*DD'+SS;
- DQDR=DD*QQ*DD'+RR;
- I_DQDR=inv(DQDR);
- AQDQ=AQDS*I_DQDR;
- ff=AA-AQDQ*DD;
- hh=AA*QQ*AA'-AQDQ*AQDS';%*(DD*QQ*AA'+SS');
- rr=DD*QQ*DD'+RR;
- ZSIG0=disc_riccati_fast(ff,DD,rr,hh);
- PP=ZSIG0 +QQ;
-
- exo_names=M_.exo_names(M_.exo_names_orig_ord,:);
-
- DPDR=DD*PP*DD'+RR;
- I_DPDR=inv(DPDR);
- %GG=[ CCCC, zeros(pd,NETA); -nmat*CCCC, zeros(NETA,NETA)];
- PDIDPDRD=PP*DD'*I_DPDR*DD;
- MSIG=disclyap_fast(CCCC, CCCC*PDIDPDRD*PP*CCCC');
-
- COV_P=[ PP, PP; PP, PP+MSIG]; % P0
-
- dr.PI_GG=[CCCC (AA-CCCC)*(eye(ss-FL_RANK)-PDIDPDRD); zeros(ss-FL_RANK) AA*(eye(ss-FL_RANK)-PDIDPDRD)];
-
- GAM= [ AA*(eye(ss-FL_RANK)-PDIDPDRD) zeros(ss-FL_RANK); (AA-CCCC)*(eye(ss-FL_RANK)-PDIDPDRD), CCCC];
-
- VV = [ dr.PI_TT1 dr.PI_TT2];
- nn=size(VV,1);
- COV_OMEGA= COV_P( end-nn+1:end, end-nn+1:end);
- COV_YR0= VV*COV_OMEGA*VV';
- diagCovYR0=diag(COV_YR0);
- labels = deblank(M_.endo_names(ivar,:));
-
- if options_.nomoments == 0
- z = [ sqrt(diagCovYR0(ivar)) diagCovYR0(ivar) ];
- title='MOMENTS OF SIMULATED VARIABLES';
- headers=strvcat('VARIABLE','STD. DEV.','VARIANCE');
- dyntable(title,headers,labels,z,size(labels,2)+2,16,10);
- end
- if options_.nocorr == 0
- diagSqrtCovYR0=sqrt(diagCovYR0);
- %COR_Y= diag(diagSqrtCovYR0)*COV_YR0*diag(diagSqrtCovYR0);
- DELTA=inv(diag(diagSqrtCovYR0));
- COR_Y= DELTA*COV_YR0*DELTA;
- title = 'CORRELATION OF SIMULATED VARIABLES';
- headers = strvcat('VARIABLE',M_.endo_names(ivar,:));
- dyntable(title,headers,labels,COR_Y(ivar,ivar),size(labels,2)+2,8,4);
- else
- COR_Y=[];
- end
-
- ar = options_.ar;
- options_ = set_default_option(options_,'ar',5);
- ar = options_.ar;
- if ar > 0
- COV_YRk= zeros(nn,ar);
- AutoCOR_YRk= zeros(nn,ar);
- for k=1:ar;
- COV_P=GAM*COV_P;
- COV_OMEGA= COV_P( end-nn+1:end, end-nn+1:end);
- COV_YRk = VV*COV_OMEGA*VV';
- AutoCOR_YRkMAT=DELTA*COV_YRk*DELTA;
- oo_.autocorr{k}=AutoCOR_YRkMAT(ivar,ivar);
- AutoCOR_YRk(:,k)= diag(COV_YRk)./diagCovYR0;
- end
- title = 'AUTOCORRELATION OF SIMULATED VARIABLES';
- headers = strvcat('VARIABLE',int2str([1:ar]'));
- dyntable(title,headers,labels,AutoCOR_YRk(ivar,:),size(labels,2)+2,8,4);
- else
- AutoCOR_YRk=[];
- end
- save ([M_.fname '_PCL_moments'], 'COV_YR0','AutoCOR_YRk', 'COR_Y');
- warning(warning_old_state);
+function [irfmat,irfst]=PCL_Part_info_moments( H, varobs, dr,ivar)
+% sets up parameters and calls part-info kalman filter
+% developed by G Perendia, July 2006 for implementation from notes by Prof. Joe Pearlman to
+% suit partial information RE solution in accordance with, and based on, the
+% Pearlman, Currie and Levine 1986 solution.
+% 22/10/06 - Version 2 for new Riccati with 4 params instead 5
+
+% Copyright (C) 2001-20010 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 .
+
+ % Recall that the state space is given by the
+ % predetermined variables s(t-1), x(t-1)
+ % and the jump variables x(t).
+ % The jump variables have dimension NETA
+
+ global M_ options_ oo_
+ warning_old_state = warning;
+ warning off
+
+ [junk,OBS] = ismember(varobs,M_.endo_names,'rows');
+
+ G1=dr.PI_ghx;
+ impact=dr.PI_ghu;
+ nmat=dr.PI_nmat;
+ CC=dr.PI_CC;
+ NX=M_.exo_nbr; % no of exogenous varexo shock variables.
+% NETA=dr.nfwrd+dr.nboth; % total no of exp. errors set to no of forward looking equations
+ FL_RANK=dr.PI_FL_RANK;
+ NY=M_.endo_nbr;
+ if isempty(OBS)
+ NOBS=NY;
+ LL=eye(NY,NY);
+ else %and if no obsevations specify OBS=[0] but this is not going to work properly
+ NOBS=length(OBS);
+ LL=zeros(NOBS,NY);
+ for i=1:NOBS
+ LL(i,OBS(i))=1;
+ end
+ end
+
+ if exist( 'irfpers')==1
+ if ~isempty(irfpers)
+ if irfpers<=0, irfpers=20, end;
+ else
+ irfpers=20;
+ end
+ else
+ irfpers=20;
+ end
+
+ ss=size(G1,1);
+
+ pd=ss-size(nmat,1);
+ SDX=M_.Sigma_e^0.5; % =SD,not V-COV, of Exog shocks or M_.Sigma_e^0.5 num_exog x num_exog matrix
+ if isempty(H)
+ H=M_.H;
+ end
+ VV=H; % V-COV of observation errors.
+ MM=impact*SDX; % R*(Q^0.5) in standard KF notation
+ % observation vector indices
+ % mapping to endogenous variables.
+
+ L1=LL*dr.PI_TT1;
+ L2=LL*dr.PI_TT2;
+
+ MM1=MM(1:ss-FL_RANK,:);
+ U11=MM1*MM1';
+ % SDX
+
+ U22=0;
+ % determine K1 and K2 observation mapping matrices
+ % This uses the fact that measurements are given by L1*s(t)+L2*x(t)
+ % and s(t) is expressed in the dynamics as
+ % H1*eps(t)+G11*s(t-1)+G12*x(t-1)+G13*x(t).
+ % Thus the observations o(t) can be written in the form
+ % o(t)=K1*[eps(t)' s(t-1)' x(t-1)']' + K2*x(t) where
+ % K1=[L1*H1 L1*G11 L1*G12] K2=L1*G13+L2
+
+ G12=G1(NX+1:ss-2*FL_RANK,:);
+ KK1=L1*G12;
+ K1=KK1(:,1:ss-FL_RANK);
+ K2=KK1(:,ss-FL_RANK+1:ss)+L2;
+
+ %pre calculate time-invariant factors
+ A11=G1(1:pd,1:pd);
+ A22=G1(pd+1:end, pd+1:end);
+ A12=G1(1:pd, pd+1:end);
+ A21=G1(pd+1:end,1:pd);
+ Lambda= nmat*A12+A22;
+ %A11_A12Nmat= A11-A12*nmat % test
+ I_L=inv(Lambda);
+ BB=A12*inv(A22);
+ FF=K2*inv(A22);
+ QQ=BB*U22*BB' + U11;
+ UFT=U22*FF';
+ % kf_param structure:
+ AA=A11-BB*A21;
+ CCCC=A11-A12*nmat; % F in new notation
+ DD=K1-FF*A21; % H in new notation
+ EE=K1-K2*nmat;
+ RR=FF*UFT+VV;
+ if ~any(RR)
+ % if zero add some dummy measurement err. variance-covariances
+ % with diagonals 0.000001. This would not be needed if we used
+ % the slow solver, or the generalised eigenvalue approach,
+ % but these are both slower.
+ RR=eye(size(RR,1))*1.0e-6;
+ end
+ SS=BB*UFT;
+ VKLUFT=VV+K2*I_L*UFT;
+ ALUFT=A12*I_L*UFT;
+ FULKV=FF*U22*I_L'*K2'+VV;
+ FUBT=FF*U22*BB';
+ nmat=nmat;
+ % initialise pshat
+ AQDS=AA*QQ*DD'+SS;
+ DQDR=DD*QQ*DD'+RR;
+ I_DQDR=inv(DQDR);
+ AQDQ=AQDS*I_DQDR;
+ ff=AA-AQDQ*DD;
+ hh=AA*QQ*AA'-AQDQ*AQDS';%*(DD*QQ*AA'+SS');
+ rr=DD*QQ*DD'+RR;
+ ZSIG0=disc_riccati_fast(ff,DD,rr,hh);
+ PP=ZSIG0 +QQ;
+
+ exo_names=M_.exo_names(M_.exo_names_orig_ord,:);
+
+ DPDR=DD*PP*DD'+RR;
+ I_DPDR=inv(DPDR);
+ %GG=[ CCCC, zeros(pd,NETA); -nmat*CCCC, zeros(NETA,NETA)];
+ PDIDPDRD=PP*DD'*I_DPDR*DD;
+ MSIG=disclyap_fast(CCCC, CCCC*PDIDPDRD*PP*CCCC');
+
+ COV_P=[ PP, PP; PP, PP+MSIG]; % P0
+
+ dr.PI_GG=[CCCC (AA-CCCC)*(eye(ss-FL_RANK)-PDIDPDRD); zeros(ss-FL_RANK) AA*(eye(ss-FL_RANK)-PDIDPDRD)];
+
+ GAM= [ AA*(eye(ss-FL_RANK)-PDIDPDRD) zeros(ss-FL_RANK); (AA-CCCC)*(eye(ss-FL_RANK)-PDIDPDRD), CCCC];
+
+ VV = [ dr.PI_TT1 dr.PI_TT2];
+ nn=size(VV,1);
+ COV_OMEGA= COV_P( end-nn+1:end, end-nn+1:end);
+ COV_YR0= VV*COV_OMEGA*VV';
+ diagCovYR0=diag(COV_YR0);
+ labels = deblank(M_.endo_names(ivar,:));
+
+ if options_.nomoments == 0
+ z = [ sqrt(diagCovYR0(ivar)) diagCovYR0(ivar) ];
+ title='MOMENTS OF SIMULATED VARIABLES';
+ headers=strvcat('VARIABLE','STD. DEV.','VARIANCE');
+ dyntable(title,headers,labels,z,size(labels,2)+2,16,10);
+ end
+ if options_.nocorr == 0
+ diagSqrtCovYR0=sqrt(diagCovYR0);
+ %COR_Y= diag(diagSqrtCovYR0)*COV_YR0*diag(diagSqrtCovYR0);
+ DELTA=inv(diag(diagSqrtCovYR0));
+ COR_Y= DELTA*COV_YR0*DELTA;
+ title = 'CORRELATION OF SIMULATED VARIABLES';
+ headers = strvcat('VARIABLE',M_.endo_names(ivar,:));
+ dyntable(title,headers,labels,COR_Y(ivar,ivar),size(labels,2)+2,8,4);
+ else
+ COR_Y=[];
+ end
+
+ ar = options_.ar;
+ options_ = set_default_option(options_,'ar',5);
+ ar = options_.ar;
+ if ar > 0
+ COV_YRk= zeros(nn,ar);
+ AutoCOR_YRk= zeros(nn,ar);
+ for k=1:ar;
+ COV_P=GAM*COV_P;
+ COV_OMEGA= COV_P( end-nn+1:end, end-nn+1:end);
+ COV_YRk = VV*COV_OMEGA*VV';
+ AutoCOR_YRkMAT=DELTA*COV_YRk*DELTA;
+ oo_.autocorr{k}=AutoCOR_YRkMAT(ivar,ivar);
+ AutoCOR_YRk(:,k)= diag(COV_YRk)./diagCovYR0;
+ end
+ title = 'AUTOCORRELATION OF SIMULATED VARIABLES';
+ headers = strvcat('VARIABLE',int2str([1:ar]'));
+ dyntable(title,headers,labels,AutoCOR_YRk(ivar,:),size(labels,2)+2,8,4);
+ else
+ AutoCOR_YRk=[];
+ end
+ save ([M_.fname '_PCL_moments'], 'COV_YR0','AutoCOR_YRk', 'COR_Y');
+ warning(warning_old_state);
diff --git a/matlab/partial_information/PI_gensys.m b/matlab/partial_information/PI_gensys.m
index 3088b81b3..9daec3ce2 100644
--- a/matlab/partial_information/PI_gensys.m
+++ b/matlab/partial_information/PI_gensys.m
@@ -1,244 +1,244 @@
-function [G1pi,C,impact,nmat,TT1,TT2,gev,eu, DD, E3, E5]=PI_gensys(a0,a1,a2,a3,c,PSI,NX,NETA,FL_RANK,M_,options_)
-% System given as
-% a0*E_t[y(t+1])+a1*y(t)=a2*y(t-1)+c+psi*eps(t)
-% with z an exogenous variable process.
-% Returned system is
-% [s(t)' x(t)' E_t x(t+1)']'=G1pi [s(t-1)' x(t-1)' x(t)]'+C+impact*eps(t),
-% and (a) the matrix nmat satisfying nmat*E_t z(t)+ E_t x(t+1)=0
-% (b) matrices TT1, TT2 that relate y(t) to these states: y(t)=[TT1 TT2][s(t)' x(t)']'.
-% Note that the dimension of the state vector = dim(a0)+NO_FL_EQS
-%
-% If div is omitted from argument list, a div>1 is calculated.
-% eu(1)=1 for existence, eu(2)=1 for uniqueness. eu(1)=-1 for
-% existence only with not-s.c. z; eu=[-2,-2] for coincident zeros.
-% Based on
-% Christopher A. Sims
-% Corrected 10/28/96 by CAS
-
-% Copyright (C) 1996-2010 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 .
-
-
-global lq_instruments;
-eu=[0;0];C=c;
-realsmall=1e-6;
-fixdiv=(nargin==6);
-n=size(a0,1);
-DD=[];E3=[]; E5=0;%[];
-%
-% Find SVD of a0, and create partitions of U, S and V
-%
-[U0,S0,V0] = svd(a0);
-
-FL_RANK=rank(S0);
-U01=U0(1:n,1:FL_RANK);
-U02=U0(1:n,FL_RANK+1:n);
-V01=V0(1:n,1:FL_RANK);
-V02=V0(1:n,FL_RANK+1:n);
-S01=S0(1:FL_RANK,1:FL_RANK);
-%
-% Define TT1, TT2
-%
-TT1=V02;
-TT2=V01;
-%
-%Invert S01
-%
-Sinv=eye(FL_RANK);
-for i=1:FL_RANK
- Sinv(i,i)=1/S01(i,i);
-end
-%
-%Set up the system matrix for the variables s(t)=V02*Y(t), x(t)=V01*Y(t) and E_t x(t+1)
-% and define z(t)'=[s(t)' x(t)]
-%
-FF=Sinv*U01'*a1*V02;
-UAVinv=inv(U02'*a1*V02);
-G11=UAVinv*U02'*a2*V02;
-G12=UAVinv*U02'*a2*V01;
-G13=-UAVinv*U02'*a1*V01;
-G31=-FF*G11+Sinv*U01'*a2*V02;
-G32=-FF*G12+Sinv*U01'*a2*V01;
-G33=-FF*G13-Sinv*U01'*a1*V01;
-H1=UAVinv*U02'*PSI;
-H3=-FF*H1+Sinv*U01'*PSI; % This should equal 0
-G21=zeros(FL_RANK,n-FL_RANK);
-G22=zeros(FL_RANK,FL_RANK);
-G23=eye(FL_RANK);
-%H2=zeros(FL_RANK,NX);
-num_inst=0;
-P1=H1;
-P3=H3;
-if(options_.ACES_solver==1 & isfield(lq_instruments,'names'))
- num_inst=size(lq_instruments.names,1);
- if ~isfield(lq_instruments,'inst_varexo_indices') & num_inst>0
- for i=1:num_inst
- i_tmp = strmatch(deblank(lq_instruments.names(i,:)),M_.exo_names,'exact');
- if isempty(i_tmp)
- error (['One of the specified instrument variables does not exist']) ;
- else
- i_var(i) = i_tmp;
- end
- end
- lq_instruments.inst_varexo_indices=i_var;
- elseif size(lq_instruments.inst_varexo_indices)>0
- i_var=lq_instruments.inst_varexo_indices;
- if ~num_inst
- num_inst=size(lq_instruments.inst_varexo_indices);
- end
- else
- i_var=[];
- num_inst=0;
- end
- if size(i_var,2)>0 & size(i_var,2)==num_inst
- N1=H1(:,i_var);
- N3=H3(:,i_var);
- x_var=zeros(NX,1);
- for i=1:NX
- if isempty(find(i_var==i))
- x_var(i)=i;
- end
- end
- x_var=nonzeros(x_var);
- P1=H1(:,x_var);
- P3=H3(:,x_var);
- NX=NX-num_inst;
- else
- error('WARNING: There are no instrumnets for ACES!');
- end
- lq_instruments.N1=N1;
- lq_instruments.N3=N3;
-elseif(options_.ACES_solver==1)
- error('WARNING: There are no instrumnets for ACES!');
-end
-% New Definitions
-Ze11=zeros(NX,NX);
-Ze12=zeros(NX,n-FL_RANK);
-Ze134=zeros(NX,FL_RANK);
-Ze31=zeros(FL_RANK,NX);
-
-% End of New Definitions
-
-%
-% System matrix is called 'G1pi'; Shock matrix is called 'impact'
-%
-
-G1pi=[Ze11 Ze12 Ze134 Ze134; P1 G11 G12 G13; Ze31 G21 G22 G23; P3 G31 G32 G33];
-
-impact=[eye(NX,NX); zeros(n+FL_RANK,NX)];
-
-if(options_.ACES_solver==1)
- E3=V02*[P1 G11 G12 G13];
- E3=E3+ [zeros(size(V01,1),size(E3,2)-size(V01,2)) V01];
- E5=-V02*N1;
- DD=[zeros(NX,size(N1,2));N1; zeros(FL_RANK,size(N1,2));N3];
- eu =[1; 1], nmat=[], gev=[];
- return; % do not check B&K compliancy
-end
-
-G0pi=eye(n+FL_RANK+NX);
-try
- [a b q z v]=qz(G0pi,G1pi);
-catch
- try
- lerror=lasterror;
- disp(['PI_Gensys: ' lerror.message]);
- if 0==strcmp('MATLAB:qz:matrixWithNaNInf',lerror.identifier)
- disp '** Unexpected Error PI_Gensys:qz: ** :';
- button=questdlg('Continue Y/N?','Unexpected Error in qz','No','Yes','Yes');
- switch button
- case 'No'
- error ('Terminated')
- %case 'Yes'
-
- end
- end
- G1pi=[];impact=[];nmat=[]; gev=[];
- eu=[-2;-2];
- return
- catch
- disp '** Unexpected Error in qz ** :';
- disp lerror.message;
- button=questdlg('Continue Y/N?','Unexpected Error in qz','No','Yes','Yes');
- switch button
- case 'No'
- error ('Terminated')
- case 'Yes'
- G1pi=[];impact=[];nmat=[]; gev=[];
- eu=[-2;-2];
- return
- end
- end
-end
-
-if ~fixdiv, div=1.01; end
-nunstab=0;
-zxz=0;
-nn=size(a,1);
-for i=1:nn
-% ------------------div calc------------
- if ~fixdiv
- if abs(a(i,i)) > 0
- divhat=abs(b(i,i))/abs(a(i,i));
- % bug detected by Vasco Curdia and Daria Finocchiaro, 2/25/2004 A root of
- % exactly 1.01 and no root between 1 and 1.02, led to div being stuck at 1.01
- % and the 1.01 root being misclassified as stable. Changing < to <= below fixes this.
- if 1+realsmalldiv*abs(a(i,i)));
- if abs(a(i,i))1 is calculated.
+% eu(1)=1 for existence, eu(2)=1 for uniqueness. eu(1)=-1 for
+% existence only with not-s.c. z; eu=[-2,-2] for coincident zeros.
+% Based on
+% Christopher A. Sims
+% Corrected 10/28/96 by CAS
+
+% Copyright (C) 1996-2010 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 .
+
+
+global lq_instruments;
+eu=[0;0];C=c;
+realsmall=1e-6;
+fixdiv=(nargin==6);
+n=size(a0,1);
+DD=[];E3=[]; E5=0;%[];
+%
+% Find SVD of a0, and create partitions of U, S and V
+%
+[U0,S0,V0] = svd(a0);
+
+FL_RANK=rank(S0);
+U01=U0(1:n,1:FL_RANK);
+U02=U0(1:n,FL_RANK+1:n);
+V01=V0(1:n,1:FL_RANK);
+V02=V0(1:n,FL_RANK+1:n);
+S01=S0(1:FL_RANK,1:FL_RANK);
+%
+% Define TT1, TT2
+%
+TT1=V02;
+TT2=V01;
+%
+%Invert S01
+%
+Sinv=eye(FL_RANK);
+for i=1:FL_RANK
+ Sinv(i,i)=1/S01(i,i);
+end
+%
+%Set up the system matrix for the variables s(t)=V02*Y(t), x(t)=V01*Y(t) and E_t x(t+1)
+% and define z(t)'=[s(t)' x(t)]
+%
+FF=Sinv*U01'*a1*V02;
+UAVinv=inv(U02'*a1*V02);
+G11=UAVinv*U02'*a2*V02;
+G12=UAVinv*U02'*a2*V01;
+G13=-UAVinv*U02'*a1*V01;
+G31=-FF*G11+Sinv*U01'*a2*V02;
+G32=-FF*G12+Sinv*U01'*a2*V01;
+G33=-FF*G13-Sinv*U01'*a1*V01;
+H1=UAVinv*U02'*PSI;
+H3=-FF*H1+Sinv*U01'*PSI; % This should equal 0
+G21=zeros(FL_RANK,n-FL_RANK);
+G22=zeros(FL_RANK,FL_RANK);
+G23=eye(FL_RANK);
+%H2=zeros(FL_RANK,NX);
+num_inst=0;
+P1=H1;
+P3=H3;
+if(options_.ACES_solver==1 & isfield(lq_instruments,'names'))
+ num_inst=size(lq_instruments.names,1);
+ if ~isfield(lq_instruments,'inst_varexo_indices') & num_inst>0
+ for i=1:num_inst
+ i_tmp = strmatch(deblank(lq_instruments.names(i,:)),M_.exo_names,'exact');
+ if isempty(i_tmp)
+ error (['One of the specified instrument variables does not exist']) ;
+ else
+ i_var(i) = i_tmp;
+ end
+ end
+ lq_instruments.inst_varexo_indices=i_var;
+ elseif size(lq_instruments.inst_varexo_indices)>0
+ i_var=lq_instruments.inst_varexo_indices;
+ if ~num_inst
+ num_inst=size(lq_instruments.inst_varexo_indices);
+ end
+ else
+ i_var=[];
+ num_inst=0;
+ end
+ if size(i_var,2)>0 & size(i_var,2)==num_inst
+ N1=H1(:,i_var);
+ N3=H3(:,i_var);
+ x_var=zeros(NX,1);
+ for i=1:NX
+ if isempty(find(i_var==i))
+ x_var(i)=i;
+ end
+ end
+ x_var=nonzeros(x_var);
+ P1=H1(:,x_var);
+ P3=H3(:,x_var);
+ NX=NX-num_inst;
+ else
+ error('WARNING: There are no instrumnets for ACES!');
+ end
+ lq_instruments.N1=N1;
+ lq_instruments.N3=N3;
+elseif(options_.ACES_solver==1)
+ error('WARNING: There are no instrumnets for ACES!');
+end
+% New Definitions
+Ze11=zeros(NX,NX);
+Ze12=zeros(NX,n-FL_RANK);
+Ze134=zeros(NX,FL_RANK);
+Ze31=zeros(FL_RANK,NX);
+
+% End of New Definitions
+
+%
+% System matrix is called 'G1pi'; Shock matrix is called 'impact'
+%
+
+G1pi=[Ze11 Ze12 Ze134 Ze134; P1 G11 G12 G13; Ze31 G21 G22 G23; P3 G31 G32 G33];
+
+impact=[eye(NX,NX); zeros(n+FL_RANK,NX)];
+
+if(options_.ACES_solver==1)
+ E3=V02*[P1 G11 G12 G13];
+ E3=E3+ [zeros(size(V01,1),size(E3,2)-size(V01,2)) V01];
+ E5=-V02*N1;
+ DD=[zeros(NX,size(N1,2));N1; zeros(FL_RANK,size(N1,2));N3];
+ eu =[1; 1], nmat=[], gev=[];
+ return; % do not check B&K compliancy
+end
+
+G0pi=eye(n+FL_RANK+NX);
+try
+ [a b q z v]=qz(G0pi,G1pi);
+catch
+ try
+ lerror=lasterror;
+ disp(['PI_Gensys: ' lerror.message]);
+ if 0==strcmp('MATLAB:qz:matrixWithNaNInf',lerror.identifier)
+ disp '** Unexpected Error PI_Gensys:qz: ** :';
+ button=questdlg('Continue Y/N?','Unexpected Error in qz','No','Yes','Yes');
+ switch button
+ case 'No'
+ error ('Terminated')
+ %case 'Yes'
+
+ end
+ end
+ G1pi=[];impact=[];nmat=[]; gev=[];
+ eu=[-2;-2];
+ return
+ catch
+ disp '** Unexpected Error in qz ** :';
+ disp lerror.message;
+ button=questdlg('Continue Y/N?','Unexpected Error in qz','No','Yes','Yes');
+ switch button
+ case 'No'
+ error ('Terminated')
+ case 'Yes'
+ G1pi=[];impact=[];nmat=[]; gev=[];
+ eu=[-2;-2];
+ return
+ end
+ end
+end
+
+if ~fixdiv, div=1.01; end
+nunstab=0;
+zxz=0;
+nn=size(a,1);
+for i=1:nn
+% ------------------div calc------------
+ if ~fixdiv
+ if abs(a(i,i)) > 0
+ divhat=abs(b(i,i))/abs(a(i,i));
+ % bug detected by Vasco Curdia and Daria Finocchiaro, 2/25/2004 A root of
+ % exactly 1.01 and no root between 1 and 1.02, led to div being stuck at 1.01
+ % and the 1.01 root being misclassified as stable. Changing < to <= below fixes this.
+ if 1+realsmalldiv*abs(a(i,i)));
+ if abs(a(i,i)).
-
-if nargin == 4 | isempty( ch ) == 1
- flag_ch = 0;
-else
- flag_ch = 1;
-end
-
-
-% intialisation
-tol = 1e-10; % iteration convergence threshold
-P0=H;
-X0=F;
-if ~any(R) % i.e. ==0
- warning('Dangerously evading inversion of zero matrix!');
- Y0=zeros(size(R));
-else
- Y0=D'*inv(R)*D;
-end
-POYO=P0*Y0;
-I=eye(size(POYO));
-clear POYO;
-
-% iterate Riccati equation solution
-matd=1;
-count=0;
-while matd > tol && count < 100
- INVPY=inv(I+P0*Y0);
- P1=X0*INVPY*P0*X0'+ P0;
- Y1=X0'*Y0*INVPY*X0+ Y0;
- X1=X0*INVPY*X0;
- matd=sum( sum(abs( P1 - P0 )));
-% P0=(P1+P1')/2
- P0=P1;
- X0=X1;
- Y0=Y1;
- count=count+1;
-% matd;
-end
-
-Z=(P0+P0')/2;
-%Z=P0
-% check if the convergence took place
-if count==100
- matd
- error('Riccati not converged fast enough!');
-% error.identifier='Riccati not converged!'
-% error
-end
-%if count >5
-% disp('Riccati count= ');
-% count
-%end
-
-clear X0 X1 Y0 Y1 P1 I INVPY;
-
-% Check that X is positive definite
-if flag_ch==1
- [C,p]=chol(Z);
- if p ~= 0
- error('Z is not positive definite')
- end
-end
+ function Z=disc_riccati_fast(F,D,R,H,ch)
+% 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
+% =================================================================
+
+% Copyright (C) 1996-2010 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 nargin == 4 | isempty( ch ) == 1
+ flag_ch = 0;
+else
+ flag_ch = 1;
+end
+
+
+% intialisation
+tol = 1e-10; % iteration convergence threshold
+P0=H;
+X0=F;
+if ~any(R) % i.e. ==0
+ warning('Dangerously evading inversion of zero matrix!');
+ Y0=zeros(size(R));
+else
+ Y0=D'*inv(R)*D;
+end
+POYO=P0*Y0;
+I=eye(size(POYO));
+clear POYO;
+
+% iterate Riccati equation solution
+matd=1;
+count=0;
+while matd > tol && count < 100
+ INVPY=inv(I+P0*Y0);
+ P1=X0*INVPY*P0*X0'+ P0;
+ Y1=X0'*Y0*INVPY*X0+ Y0;
+ X1=X0*INVPY*X0;
+ matd=sum( sum(abs( P1 - P0 )));
+% P0=(P1+P1')/2
+ P0=P1;
+ X0=X1;
+ Y0=Y1;
+ count=count+1;
+% matd;
+end
+
+Z=(P0+P0')/2;
+%Z=P0
+% check if the convergence took place
+if count==100
+ matd
+ error('Riccati not converged fast enough!');
+% error.identifier='Riccati not converged!'
+% error
+end
+%if count >5
+% disp('Riccati count= ');
+% count
+%end
+
+clear X0 X1 Y0 Y1 P1 I INVPY;
+
+% Check that X is positive definite
+if flag_ch==1
+ [C,p]=chol(Z);
+ if p ~= 0
+ error('Z is not positive definite')
+ end
+end
diff --git a/tests/partial_information/PItest3aHc0PCLsimModPiYrVarobsAll.mod b/tests/partial_information/PItest3aHc0PCLsimModPiYrVarobsAll.mod
index 68d32cee9..f4418db14 100644
--- a/tests/partial_information/PItest3aHc0PCLsimModPiYrVarobsAll.mod
+++ b/tests/partial_information/PItest3aHc0PCLsimModPiYrVarobsAll.mod
@@ -1,47 +1,47 @@
-% inflation target (pitarg) modelled as AR1
-% cy = 0.614479/0.769365 - obtained from the computed steady state from the nonlinear counterpart
-
-
-var pi mc mun muc c y n r g a pitarg ;
-varexo eps_g eps_a eps_e eps_m eps_targ;
-
-parameters beta xi hc wd sigma gamma rho_g rho_a rho_r rho_targ thetap cy varrho;
-beta = 0.99;
-xi = 0.5034;
-hc = 0.0;
-wd = 0.40;
-gamma = 0.5868;
-sigma = 4.0897;
-cy = 0.614479/0.769365;
-rho_g = 0.8325;
-rho_a = 0.9827;
-rho_r = 0.3529;
-thetap = 2.2161;
-varrho = 0.3853;
-rho_targ = 0.6133;
-
-model(linear);
-pi = (beta/(1+beta*gamma))*pi(+1)+(gamma/(1+beta*gamma))*pi(-1)+(((1-beta*xi)*(1-xi))/((1+beta*gamma)*xi))*(mc+eps_m);
-mc = mun-muc-a;
-mun = (c-hc*c(-1))/(1-hc)+wd*n/(1-wd)+muc;
-muc = ((1-varrho)*(1-sigma)-1)*(c-hc*c(-1))/(1-hc)-wd*varrho*(1-sigma)*n/(1-wd);
-muc(+1) = muc-(r-pi(+1));
-y = cy*c+(1-cy)*g;
-n = y-a;
-r = rho_r*r(-1)+thetap*(1-rho_r)*(pi(+1)-rho_targ*pitarg)+eps_e;
-g = rho_g*g(-1)+eps_g;
-a = rho_a*a(-1)+eps_a;
-pitarg = rho_targ*pitarg(-1)+eps_targ;
-end;
-
-shocks;
-var eps_g; stderr 3.8505;
-var eps_a; stderr 0.7573;
-var eps_e; stderr 0.2409;
-var eps_m; stderr 0.8329;
-var eps_targ; stderr 0.3978;
-end;
-
-varobs pi mc mun muc c y n r g a pitarg;
-stoch_simul(partial_information,irf=30)pi y r;//pi n c y r;
-
+% inflation target (pitarg) modelled as AR1
+% cy = 0.614479/0.769365 - obtained from the computed steady state from the nonlinear counterpart
+
+
+var pi mc mun muc c y n r g a pitarg ;
+varexo eps_g eps_a eps_e eps_m eps_targ;
+
+parameters beta xi hc wd sigma gamma rho_g rho_a rho_r rho_targ thetap cy varrho;
+beta = 0.99;
+xi = 0.5034;
+hc = 0.0;
+wd = 0.40;
+gamma = 0.5868;
+sigma = 4.0897;
+cy = 0.614479/0.769365;
+rho_g = 0.8325;
+rho_a = 0.9827;
+rho_r = 0.3529;
+thetap = 2.2161;
+varrho = 0.3853;
+rho_targ = 0.6133;
+
+model(linear);
+pi = (beta/(1+beta*gamma))*pi(+1)+(gamma/(1+beta*gamma))*pi(-1)+(((1-beta*xi)*(1-xi))/((1+beta*gamma)*xi))*(mc+eps_m);
+mc = mun-muc-a;
+mun = (c-hc*c(-1))/(1-hc)+wd*n/(1-wd)+muc;
+muc = ((1-varrho)*(1-sigma)-1)*(c-hc*c(-1))/(1-hc)-wd*varrho*(1-sigma)*n/(1-wd);
+muc(+1) = muc-(r-pi(+1));
+y = cy*c+(1-cy)*g;
+n = y-a;
+r = rho_r*r(-1)+thetap*(1-rho_r)*(pi(+1)-rho_targ*pitarg)+eps_e;
+g = rho_g*g(-1)+eps_g;
+a = rho_a*a(-1)+eps_a;
+pitarg = rho_targ*pitarg(-1)+eps_targ;
+end;
+
+shocks;
+var eps_g; stderr 3.8505;
+var eps_a; stderr 0.7573;
+var eps_e; stderr 0.2409;
+var eps_m; stderr 0.8329;
+var eps_targ; stderr 0.3978;
+end;
+
+varobs pi mc mun muc c y n r g a pitarg;
+stoch_simul(partial_information,irf=30)pi y r;//pi n c y r;
+
diff --git a/tests/partial_information/PItest3aHc0PCLsimModPiYrVarobsCNR.mod b/tests/partial_information/PItest3aHc0PCLsimModPiYrVarobsCNR.mod
index 07a324078..0a0f84645 100644
--- a/tests/partial_information/PItest3aHc0PCLsimModPiYrVarobsCNR.mod
+++ b/tests/partial_information/PItest3aHc0PCLsimModPiYrVarobsCNR.mod
@@ -1,47 +1,47 @@
-% inflation target (pitarg) modelled as AR1
-% cy = 0.614479/0.769365 - obtained from the computed steady state from the nonlinear counterpart
-
-
-var pi mc mun muc c y n r g a pitarg;
-varexo eps_g eps_a eps_e eps_m eps_targ;
-
-parameters beta xi hc wd sigma gamma rho_g rho_a rho_r rho_targ thetap cy varrho;
-beta = 0.99;
-xi = 0.5034;
-hc = 0.0;
-wd = 0.40;
-gamma = 0.5868;
-sigma = 4.0897;
-cy = 0.614479/0.769365;
-rho_g = 0.8325;
-rho_a = 0.9827;
-rho_r = 0.3529;
-thetap = 2.2161;
-varrho = 0.3853;
-rho_targ = 0.6133;
-
-model(linear);
-pi = (beta/(1+beta*gamma))*pi(+1)+(gamma/(1+beta*gamma))*pi(-1)+(((1-beta*xi)*(1-xi))/((1+beta*gamma)*xi))*(mc+eps_m);
-mc = mun-muc-a;
-mun = (c-hc*c(-1))/(1-hc)+wd*n/(1-wd)+muc;
-muc = ((1-varrho)*(1-sigma)-1)*(c-hc*c(-1))/(1-hc)-wd*varrho*(1-sigma)*n/(1-wd);
-muc(+1) = muc-(r-pi(+1));
-y = cy*c+(1-cy)*g;
-n = y-a;
-r = rho_r*r(-1)+thetap*(1-rho_r)*(pi(+1)-rho_targ*pitarg)+eps_e;
-g = rho_g*g(-1)+eps_g;
-a = rho_a*a(-1)+eps_a;
-pitarg = rho_targ*pitarg(-1)+eps_targ;
-end;
-
-shocks;
-var eps_g; stderr 3.8505;
-var eps_a; stderr 0.7573;
-var eps_e; stderr 0.2409;
-var eps_m; stderr 0.8329;
-var eps_targ; stderr 0.3978;
-end;
-
-varobs c n r ;
-stoch_simul(partial_information,irf=30)pi y r;//pi n c y r;
-
+% inflation target (pitarg) modelled as AR1
+% cy = 0.614479/0.769365 - obtained from the computed steady state from the nonlinear counterpart
+
+
+var pi mc mun muc c y n r g a pitarg;
+varexo eps_g eps_a eps_e eps_m eps_targ;
+
+parameters beta xi hc wd sigma gamma rho_g rho_a rho_r rho_targ thetap cy varrho;
+beta = 0.99;
+xi = 0.5034;
+hc = 0.0;
+wd = 0.40;
+gamma = 0.5868;
+sigma = 4.0897;
+cy = 0.614479/0.769365;
+rho_g = 0.8325;
+rho_a = 0.9827;
+rho_r = 0.3529;
+thetap = 2.2161;
+varrho = 0.3853;
+rho_targ = 0.6133;
+
+model(linear);
+pi = (beta/(1+beta*gamma))*pi(+1)+(gamma/(1+beta*gamma))*pi(-1)+(((1-beta*xi)*(1-xi))/((1+beta*gamma)*xi))*(mc+eps_m);
+mc = mun-muc-a;
+mun = (c-hc*c(-1))/(1-hc)+wd*n/(1-wd)+muc;
+muc = ((1-varrho)*(1-sigma)-1)*(c-hc*c(-1))/(1-hc)-wd*varrho*(1-sigma)*n/(1-wd);
+muc(+1) = muc-(r-pi(+1));
+y = cy*c+(1-cy)*g;
+n = y-a;
+r = rho_r*r(-1)+thetap*(1-rho_r)*(pi(+1)-rho_targ*pitarg)+eps_e;
+g = rho_g*g(-1)+eps_g;
+a = rho_a*a(-1)+eps_a;
+pitarg = rho_targ*pitarg(-1)+eps_targ;
+end;
+
+shocks;
+var eps_g; stderr 3.8505;
+var eps_a; stderr 0.7573;
+var eps_e; stderr 0.2409;
+var eps_m; stderr 0.8329;
+var eps_targ; stderr 0.3978;
+end;
+
+varobs c n r ;
+stoch_simul(partial_information,irf=30)pi y r;//pi n c y r;
+