diff --git a/matlab/partial_information/PCL_Part_info_irf.m b/matlab/partial_information/PCL_Part_info_irf.m index 2b8e40dee..81277b14f 100644 --- a/matlab/partial_information/PCL_Part_info_irf.m +++ b/matlab/partial_information/PCL_Part_info_irf.m @@ -1,11 +1,11 @@ -function [irfmat,irfst]=PCL_Part_info_irf( H, varobs, M_, dr, irfpers,ii) +function y=PCL_Part_info_irf( H, varobs, ivar, 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) 2006-2010 Dynare Team +% Copyright (C) 2001-20010 Dynare Team % % This file is part of Dynare. % @@ -49,16 +49,6 @@ function [irfmat,irfst]=PCL_Part_info_irf( H, varobs, M_, dr, irfpers,ii) 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 @@ -86,8 +76,8 @@ function [irfmat,irfst]=PCL_Part_info_irf( H, varobs, M_, dr, irfpers,ii) % 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; + 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; @@ -97,13 +87,11 @@ function [irfmat,irfst]=PCL_Part_info_irf( H, varobs, M_, dr, irfpers,ii) 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 @@ -138,25 +126,22 @@ function [irfmat,irfst]=PCL_Part_info_irf( H, varobs, M_, dr, irfpers,ii) 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=zeros(size(dr.PI_TT1 ,1),irfpers+1); + irfst=zeros(size(GG,1),irfpers+1); + irfst(:,1)=stderr(ii)*imp(:,ii); + for jj=2:irfpers+1; + irfst(:,jj)=GG*irfst(:,jj-1); 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'); + end + y=zeros(M_.endo_nbr,irfpers); + y(:,:)=irfmat(:,1:irfpers); + + 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 06484b6aa..f0e42ea18 100644 --- a/matlab/partial_information/PCL_Part_info_moments.m +++ b/matlab/partial_information/PCL_Part_info_moments.m @@ -1,11 +1,11 @@ -function [irfmat,irfst]=PCL_Part_info_moments( H, varobs, dr,ivar) +function AutoCOR_YRk=PCL_Part_info_irmoments( 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) 2006-2010 Dynare Team +% Copyright (C) 2001-20010 Dynare Team % % This file is part of Dynare. % @@ -30,7 +30,7 @@ function [irfmat,irfst]=PCL_Part_info_moments( H, varobs, dr,ivar) global M_ options_ oo_ warning_old_state = warning; warning off - + [junk,OBS] = ismember(varobs,M_.endo_names,'rows'); G1=dr.PI_ghx; @@ -38,7 +38,6 @@ function [irfmat,irfst]=PCL_Part_info_moments( H, varobs, dr,ivar) 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) @@ -80,7 +79,6 @@ function [irfmat,irfst]=PCL_Part_info_moments( H, varobs, dr,ivar) 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) @@ -101,7 +99,6 @@ function [irfmat,irfst]=PCL_Part_info_moments( H, varobs, dr,ivar) 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); @@ -141,7 +138,6 @@ function [irfmat,irfst]=PCL_Part_info_moments( H, varobs, dr,ivar) 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'); @@ -166,7 +162,6 @@ function [irfmat,irfst]=PCL_Part_info_moments( H, varobs, dr,ivar) 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'; diff --git a/matlab/stoch_simul.m b/matlab/stoch_simul.m index fb3809a16..63a47e10a 100644 --- a/matlab/stoch_simul.m +++ b/matlab/stoch_simul.m @@ -30,9 +30,9 @@ elseif options_.order == 3 end if options_.partial_information == 1 || options_.ACES_solver == 1 - ACES_solver = 1; + PI_PCL_solver = 1; else - ACES_solver = 0; + PI_PCL_solver = 0; end TeX = options_.TeX; @@ -50,7 +50,7 @@ end check_model; -if ACES_solver +if PI_PCL_solver [oo_.dr, info] = PCL_resol(oo_.steady_state,0); else [oo_.dr, info] = resol(oo_.steady_state,0); @@ -82,18 +82,26 @@ if ~options_.noprint disp(' ') disp(' SOLUTION UNDER PARTIAL INFORMATION') disp(' ') - disp(' OBSERVED VARIABLES') - for i=1:size(options_.varobs,1) - disp([' ' options_.varobs(i,:)]) + + if isfield(options_,'varobs')&& ~isempty(options_.varobs) + PCL_varobs=options_.varobs; + disp(' OBSERVED VARIABLES') + else + PCL_varobs=var_list; + disp(' VAROBS LIST NOT SPECIFIED') + disp(' ASSUMED OBSERVED VARIABLES') + end + for i=1:size(PCL_varobs,1) + disp([' ' PCL_varobs(i,:)]) end end disp(' ') - if options_.order <= 2 && ~ACES_solver + if options_.order <= 2 && ~PI_PCL_solver disp_dr(oo_.dr,options_.order,var_list); end end -if options_.periods > 0 && ~ACES_solver +if options_.periods > 0 && ~PI_PCL_solver if options_.periods < options_.drop disp(['STOCH_SIMUL error: The horizon of simulation is shorter' ... ' than the number of observations to be DROPed']) @@ -105,8 +113,8 @@ if options_.periods > 0 && ~ACES_solver end if options_.nomoments == 0 - if ACES_solver - PCL_Part_info_moments (0, options_.varobs, oo_.dr, i_var); + if PI_PCL_solver + PCL_Part_info_moments (0, PCL_varobs, oo_.dr, i_var); elseif options_.periods == 0 disp_th_moments(oo_.dr,var_list); else @@ -133,8 +141,8 @@ if options_.irf end for i=1:M_.exo_nbr if SS(i,i) > 1e-13 - if ACES_solver - y=PCL_Part_info_irf (0, options_.varobs, M_, oo_.dr, options_.irf, i); + if PI_PCL_solver + y=PCL_Part_info_irf (0, PCL_varobs, i_var, M_, oo_.dr, options_.irf, i); else y=irf(oo_.dr,cs(M_.exo_names_orig_ord,i), options_.irf, options_.drop, ... options_.replic, options_.order);