PCL_Part_info_moments

 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

0001 function  AutoCOR_YRk=PCL_Part_info_moments( H, varobs, dr,ivar)
0002 % sets up parameters and calls part-info kalman filter
0003 % developed by G Perendia, July 2006 for implementation from notes by Prof. Joe Pearlman to
0004 % suit partial information RE solution in accordance with, and based on, the
0005 % Pearlman, Currie and Levine 1986 solution.
0006 % 22/10/06 - Version 2 for new Riccati with 4 params instead 5
0007 
0008 % Copyright (C) 2006-2012 Dynare Team
0009 %
0010 % This file is part of Dynare.
0011 %
0012 % Dynare is free software: you can redistribute it and/or modify
0013 % it under the terms of the GNU General Public License as published by
0014 % the Free Software Foundation, either version 3 of the License, or
0015 % (at your option) any later version.
0016 %
0017 % Dynare is distributed in the hope that it will be useful,
0018 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0019 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0020 % GNU General Public License for more details.
0021 %
0022 % You should have received a copy of the GNU General Public License
0023 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0024 
0025 % Recall that the state space is given by the
0026 % predetermined variables s(t-1), x(t-1)
0027 % and the jump variables x(t).
0028 % The jump variables have dimension NETA
0029 
0030 global M_ options_ oo_
0031 warning_old_state = warning;
0032 warning off
0033 
0034 OBS = [];
0035 for i=1:rows(varobs)
0036     OBS = [OBS find(strcmp(deblank(varobs(i,:)), cellstr(M_.endo_names))) ];
0037 end
0038 NOBS = length(OBS);
0039 
0040 G1=dr.PI_ghx;
0041 impact=dr.PI_ghu;
0042 nmat=dr.PI_nmat;
0043 CC=dr.PI_CC;
0044 NX=M_.exo_nbr; % no of exogenous varexo shock variables.
0045 FL_RANK=dr.PI_FL_RANK;
0046 NY=M_.endo_nbr;
0047 LL = sparse(1:NOBS,OBS,ones(NOBS,1),NY,NY);
0048 
0049 if exist( 'irfpers')==1
0050     if ~isempty(irfpers)
0051         if irfpers<=0, irfpers=20, end;
0052     else
0053         irfpers=20;
0054     end
0055 else
0056     irfpers=20;
0057 end      
0058 
0059 ss=size(G1,1);
0060 
0061 pd=ss-size(nmat,1);
0062 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
0063 if isempty(H)
0064     H=M_.H;
0065 end
0066 VV=H; % V-COV of observation errors.
0067 MM=impact*SDX; % R*(Q^0.5) in standard KF notation
0068                % observation vector indices
0069                % mapping to endogenous variables.
0070 
0071 L1=LL*dr.PI_TT1;
0072 L2=LL*dr.PI_TT2;
0073 
0074 MM1=MM(1:ss-FL_RANK,:);
0075 U11=MM1*MM1';
0076 % SDX
0077 U22=0;
0078 % determine K1 and K2 observation mapping matrices
0079 % This uses the fact that measurements are given by L1*s(t)+L2*x(t)
0080 % and s(t) is expressed in the dynamics as
0081 % H1*eps(t)+G11*s(t-1)+G12*x(t-1)+G13*x(t).
0082 % Thus the observations o(t) can be written in the form
0083 % o(t)=K1*[eps(t)' s(t-1)' x(t-1)']' + K2*x(t) where
0084 % K1=[L1*H1 L1*G11 L1*G12] K2=L1*G13+L2
0085 
0086 G12=G1(NX+1:ss-2*FL_RANK,:);
0087 KK1=L1*G12;
0088 K1=KK1(:,1:ss-FL_RANK);
0089 K2=KK1(:,ss-FL_RANK+1:ss)+L2;
0090 
0091 %pre calculate time-invariant factors
0092 A11=G1(1:pd,1:pd);
0093 A22=G1(pd+1:end, pd+1:end);
0094 A12=G1(1:pd, pd+1:end);
0095 A21=G1(pd+1:end,1:pd);
0096 Lambda= nmat*A12+A22;
0097 I_L=inv(Lambda);
0098 BB=A12*inv(A22);
0099 FF=K2*inv(A22);       
0100 QQ=BB*U22*BB' + U11;        
0101 UFT=U22*FF';
0102 % kf_param structure:
0103 AA=A11-BB*A21;
0104 CCCC=A11-A12*nmat; % F in new notation
0105 DD=K1-FF*A21; % H in new notation
0106 EE=K1-K2*nmat;
0107 RR=FF*UFT+VV;
0108 if ~any(RR) 
0109     % if zero add some dummy measurement err. variance-covariances
0110     % with diagonals 0.000001. This would not be needed if we used
0111     % the slow solver, or the generalised eigenvalue approach,
0112     % but these are both slower.
0113     RR=eye(size(RR,1))*1.0e-6;
0114 end
0115 SS=BB*UFT;
0116 VKLUFT=VV+K2*I_L*UFT;
0117 ALUFT=A12*I_L*UFT;
0118 FULKV=FF*U22*I_L'*K2'+VV;
0119 FUBT=FF*U22*BB';
0120 nmat=nmat;
0121 % initialise pshat
0122 AQDS=AA*QQ*DD'+SS;
0123 DQDR=DD*QQ*DD'+RR;
0124 I_DQDR=inv(DQDR);
0125 AQDQ=AQDS*I_DQDR;
0126 ff=AA-AQDQ*DD;
0127 hh=AA*QQ*AA'-AQDQ*AQDS';%*(DD*QQ*AA'+SS');
0128 rr=DD*QQ*DD'+RR;
0129 ZSIG0=disc_riccati_fast(ff,DD,rr,hh);
0130 PP=ZSIG0 +QQ;
0131 
0132 exo_names=M_.exo_names(M_.exo_names_orig_ord,:);
0133 
0134 DPDR=DD*PP*DD'+RR;
0135 I_DPDR=inv(DPDR);
0136 PDIDPDRD=PP*DD'*I_DPDR*DD;
0137 MSIG=disclyap_fast(CCCC, CCCC*PDIDPDRD*PP*CCCC', options_.lyapunov_doubling_tol);
0138 
0139 COV_P=[ PP, PP; PP, PP+MSIG]; % P0
0140 
0141 dr.PI_GG=[CCCC (AA-CCCC)*(eye(ss-FL_RANK)-PDIDPDRD); zeros(ss-FL_RANK) AA*(eye(ss-FL_RANK)-PDIDPDRD)];
0142 
0143 GAM= [ AA*(eye(ss-FL_RANK)-PDIDPDRD) zeros(ss-FL_RANK); (AA-CCCC)*(eye(ss-FL_RANK)-PDIDPDRD),  CCCC];
0144 
0145 VV = [  dr.PI_TT1 dr.PI_TT2];
0146 nn=size(VV,1);
0147 COV_OMEGA= COV_P( end-nn+1:end, end-nn+1:end);
0148 COV_YR0= VV*COV_OMEGA*VV';
0149 diagCovYR0=diag(COV_YR0);
0150 labels = deblank(M_.endo_names(ivar,:));
0151 
0152 if options_.nomoments == 0
0153     z = [ sqrt(diagCovYR0(ivar)) diagCovYR0(ivar) ]; 
0154     title='THEORETICAL MOMENTS';
0155     headers=char('VARIABLE','STD. DEV.','VARIANCE');
0156     dyntable(title,headers,labels,z,size(labels,2)+2,16,10);
0157 end
0158 if options_.nocorr == 0
0159     diagSqrtCovYR0=sqrt(diagCovYR0);
0160     DELTA=inv(diag(diagSqrtCovYR0));
0161     COR_Y= DELTA*COV_YR0*DELTA;
0162     title = 'MATRIX OF CORRELATION';
0163     headers = char('VARIABLE',M_.endo_names(ivar,:));
0164     dyntable(title,headers,labels,COR_Y(ivar,ivar),size(labels,2)+2,8,4);
0165 else
0166     COR_Y=[];
0167 end
0168 
0169 ar = options_.ar;
0170 if ar > 0
0171     COV_YRk= zeros(nn,ar); 
0172     AutoCOR_YRk= zeros(nn,ar); 
0173     for k=1:ar;
0174         COV_P=GAM*COV_P;
0175         COV_OMEGA= COV_P( end-nn+1:end, end-nn+1:end);
0176         COV_YRk = VV*COV_OMEGA*VV';
0177         AutoCOR_YRkMAT=DELTA*COV_YRk*DELTA;
0178         oo_.autocorr{k}=AutoCOR_YRkMAT(ivar,ivar);
0179         AutoCOR_YRk(:,k)= diag(COV_YRk)./diagCovYR0;
0180     end
0181     title = 'COEFFICIENTS OF AUTOCORRELATION';
0182     headers = char('VARIABLE',int2str([1:ar]'));
0183     dyntable(title,headers,labels,AutoCOR_YRk(ivar,:),size(labels,2)+2,8,4);
0184 else
0185     AutoCOR_YRk=[];
0186 end
0187 save ([M_.fname '_PCL_moments'], 'COV_YR0','AutoCOR_YRk', 'COR_Y');
0188 warning(warning_old_state);