added new Partial Inforamtion code sent by G. Perendia
parent
9b20ae958c
commit
c5bc424dd4
|
@ -0,0 +1,160 @@
|
|||
function [irfmat,irfst]=PCL_Part_info_irf( H, OBS, 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 <http://www.gnu.org/licenses/>.
|
||||
|
||||
|
||||
% 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
|
||||
|
||||
|
||||
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=size(find(OBS),2);
|
||||
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+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); % 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');
|
|
@ -0,0 +1,197 @@
|
|||
function [irfmat,irfst]=PCL_Part_info_irf( H, OBS, dr,ivar, var_list)
|
||||
% 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 <http://www.gnu.org/licenses/>.
|
||||
|
||||
% 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
|
||||
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=size(find(OBS),2);
|
||||
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);
|
|
@ -0,0 +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 <http://www.gnu.org/licenses/>.
|
||||
|
||||
|
||||
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+realsmall<divhat & divhat<=div
|
||||
div=.5*(1+divhat);
|
||||
end
|
||||
end
|
||||
end
|
||||
% ----------------------------------------
|
||||
nunstab=nunstab+(abs(b(i,i))>div*abs(a(i,i)));
|
||||
if abs(a(i,i))<realsmall & abs(b(i,i))<realsmall
|
||||
zxz=1;
|
||||
end
|
||||
end
|
||||
div ;
|
||||
if ~zxz
|
||||
[a b q z]=qzdiv(div,a,b,q,z);
|
||||
end
|
||||
|
||||
gev=[diag(a) diag(b)];
|
||||
if zxz
|
||||
disp('Coincident zeros. Indeterminacy and/or nonexistence.')
|
||||
eu=[-2;-2];
|
||||
% correction added 7/29/2003. Otherwise the failure to set output
|
||||
% arguments leads to an error message and no output (including eu).
|
||||
nmat=[]; %;gev=[]
|
||||
return
|
||||
end
|
||||
if (FL_RANK ~= nunstab & options_.ACES_solver~=1)
|
||||
disp(['Number of unstable variables ' num2str(nunstab)]);
|
||||
disp( ['does not match number of expectational equations ' num2str(FL_RANK)]);
|
||||
nmat=[];% gev=[];
|
||||
eu=[-2;-2];
|
||||
return
|
||||
end
|
||||
|
||||
% New Definitions
|
||||
z1=z(:,1:n+NX)';
|
||||
z2=z(:,n+NX+1:n+NX+FL_RANK)';
|
||||
|
||||
% New N Matrix by J Pearlman
|
||||
z12=z2(:,1:n+NX);
|
||||
z22=z2(:,n+NX+1:n+NX+FL_RANK);
|
||||
% End of New Definitions
|
||||
|
||||
% modified by GP:
|
||||
nmat=real(inv(z22)*z12);
|
||||
eu=[1;1];
|
|
@ -0,0 +1,91 @@
|
|||
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 <http://www.gnu.org/licenses/>.
|
||||
|
||||
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
|
|
@ -0,0 +1,327 @@
|
|||
function [dr,info,M_,options_,oo_] = dr1_PI(dr,task,M_,options_,oo_)
|
||||
% function [dr,info,M_,options_,oo_] = dr1_PI(dr,task,M_,options_,oo_)
|
||||
% Computes the reduced form solution of a rational expectation model first
|
||||
% order
|
||||
% approximation of the Partial Information stochastic model solver around the deterministic steady state).
|
||||
% Prepares System as
|
||||
% A0*E_t[y(t+1])+A1*y(t)=A2*y(t-1)+c+psi*eps(t)
|
||||
% with z an exogenous variable process.
|
||||
% and calls PI_Gensys.m solver
|
||||
% based on Pearlman et al 1986 paper and derived from
|
||||
% C.Sims' gensys linear solver.
|
||||
% to return solution in format
|
||||
% [s(t)' x(t)' E_t x(t+1)']'=G1pi [s(t-1)' x(t-1)' x(t)]'+C+impact*eps(t),
|
||||
%
|
||||
% INPUTS
|
||||
% dr [matlab structure] Decision rules for stochastic simulations.
|
||||
% task [integer] if task = 0 then dr1 computes decision rules.
|
||||
% if task = 1 then dr1 computes eigenvalues.
|
||||
% M_ [matlab structure] Definition of the model.
|
||||
% options_ [matlab structure] Global options.
|
||||
% oo_ [matlab structure] Results
|
||||
%
|
||||
% OUTPUTS
|
||||
% dr [matlab structure] Decision rules for stochastic simulations.
|
||||
% info [integer] info=1: the model doesn't define current variables uniquely
|
||||
% info=2: problem in mjdgges.dll info(2) contains error code.
|
||||
% info=3: BK order condition not satisfied info(2) contains "distance"
|
||||
% absence of stable trajectory.
|
||||
% info=4: BK order condition not satisfied info(2) contains "distance"
|
||||
% indeterminacy.
|
||||
% info=5: BK rank condition not satisfied.
|
||||
% M_ [matlab structure]
|
||||
% options_ [matlab structure]
|
||||
% oo_ [matlab structure]
|
||||
%
|
||||
% ALGORITHM
|
||||
% ...
|
||||
%
|
||||
% SPECIAL REQUIREMENTS
|
||||
% none.
|
||||
%
|
||||
|
||||
% 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 <http://www.gnu.org/licenses/>.
|
||||
|
||||
info = 0;
|
||||
|
||||
options_ = set_default_option(options_,'loglinear',0);
|
||||
options_ = set_default_option(options_,'noprint',0);
|
||||
options_ = set_default_option(options_,'olr',0);
|
||||
options_ = set_default_option(options_,'olr_beta',1);
|
||||
options_ = set_default_option(options_,'qz_criterium',1.000001);
|
||||
|
||||
xlen = M_.maximum_endo_lead + M_.maximum_endo_lag + 1;
|
||||
|
||||
if (options_.aim_solver == 1)
|
||||
options_.aim_solver == 0;
|
||||
warning('You can not use AIM with Part Info solver. AIM ignored');
|
||||
end
|
||||
if (options_.order > 1)
|
||||
warning('You can not use order higher than 1 with Part Info solver. Order 1 assumed');
|
||||
options_.order =1;
|
||||
end
|
||||
|
||||
% expanding system for Optimal Linear Regulator
|
||||
if options_.ramsey_policy
|
||||
if isfield(M_,'orig_model')
|
||||
orig_model = M_.orig_model;
|
||||
M_.endo_nbr = orig_model.endo_nbr;
|
||||
M_.endo_names = orig_model.endo_names;
|
||||
M_.lead_lag_incidence = orig_model.lead_lag_incidence;
|
||||
M_.maximum_lead = orig_model.maximum_lead;
|
||||
M_.maximum_endo_lead = orig_model.maximum_endo_lead;
|
||||
M_.maximum_lag = orig_model.maximum_lag;
|
||||
M_.maximum_endo_lag = orig_model.maximum_endo_lag;
|
||||
end
|
||||
old_solve_algo = options_.solve_algo;
|
||||
% options_.solve_algo = 1;
|
||||
oo_.steady_state = dynare_solve('ramsey_static',oo_.steady_state,0,M_,options_,oo_,it_);
|
||||
options_.solve_algo = old_solve_algo;
|
||||
[junk,junk,multbar] = ramsey_static(oo_.steady_state,M_,options_,oo_,it_);
|
||||
[jacobia_,M_] = ramsey_dynamic(oo_.steady_state,multbar,M_,options_,oo_,it_);
|
||||
klen = M_.maximum_lag + M_.maximum_lead + 1;
|
||||
dr.ys = [oo_.steady_state;zeros(M_.exo_nbr,1);multbar];
|
||||
else
|
||||
klen = M_.maximum_lag + M_.maximum_lead + 1;
|
||||
iyv = M_.lead_lag_incidence';
|
||||
iyv = iyv(:);
|
||||
iyr0 = find(iyv) ;
|
||||
it_ = M_.maximum_lag + 1 ;
|
||||
|
||||
if M_.exo_nbr == 0
|
||||
oo_.exo_steady_state = [] ;
|
||||
end
|
||||
|
||||
z = repmat(dr.ys,1,klen);
|
||||
z = z(iyr0) ;
|
||||
[junk,jacobia_] = feval([M_.fname '_dynamic'],z,[oo_.exo_simul ...
|
||||
oo_.exo_det_simul], M_.params, it_);
|
||||
end
|
||||
|
||||
if options_.debug
|
||||
save([M_.fname '_debug.mat'],'jacobia_')
|
||||
end
|
||||
|
||||
dr=set_state_space(dr,M_);
|
||||
kstate = dr.kstate;
|
||||
kad = dr.kad;
|
||||
kae = dr.kae;
|
||||
nstatic = dr.nstatic;
|
||||
nfwrd = dr.nfwrd;
|
||||
npred = dr.npred;
|
||||
nboth = dr.nboth;
|
||||
order_var = dr.order_var;
|
||||
nd = size(kstate,1);
|
||||
nz = nnz(M_.lead_lag_incidence);
|
||||
|
||||
sdyn = M_.endo_nbr - nstatic;
|
||||
|
||||
k0 = M_.lead_lag_incidence(M_.maximum_endo_lag+1,order_var);
|
||||
k1 = M_.lead_lag_incidence(find([1:klen] ~= M_.maximum_endo_lag+1),:);
|
||||
b = jacobia_(:,k0);
|
||||
|
||||
if (options_.aim_solver == 1)
|
||||
error('Anderson and Moore AIM solver is not compatible with Partial Information models');
|
||||
end % end if useAIM and...
|
||||
|
||||
%forward--looking models
|
||||
if nstatic > 0
|
||||
[Q,R] = qr(b(:,1:nstatic));
|
||||
aa = Q'*jacobia_;
|
||||
else
|
||||
aa = jacobia_;
|
||||
end
|
||||
|
||||
% If required, try PCL86 solver, that is, if not the check being
|
||||
% performed only and if it is 1st order
|
||||
% create sparse, extended jacobia AA:
|
||||
nendo=M_.endo_nbr; % = size(aa,1)
|
||||
%%% OLD: aax=zeros(size(aa,1),size(aa,1)*klen);
|
||||
% partition jacobian:
|
||||
jlen=dr.nspred+dr.nsfwrd+M_.endo_nbr+M_.exo_nbr; % length of jacobian
|
||||
|
||||
|
||||
% find size xlen of the state vector Y and of A0, A1 and A2 transition matrices:
|
||||
% it is the sum the all i variables's lag/lead representations,
|
||||
% for each variable i representation being defined as:
|
||||
% Max (i_lags-1,0)+ Max (i_leads-1,0)+1
|
||||
% so that if variable x appears with 2 lags and 1 lead, and z
|
||||
% with 2 lags and 3 leads, the size of the state space is:
|
||||
% 1+0+1 + 1+2+1 =6
|
||||
% e.g. E_t Y(t+1)=
|
||||
% E_t x(t)
|
||||
% E_t x(t+1)
|
||||
% E_t z(t)
|
||||
% E_t z(t+1)
|
||||
% E_t z(t+2)
|
||||
% E_t z(t+3)
|
||||
|
||||
|
||||
% first transpose M_.lead_lag_incidence';
|
||||
lead_lag=M_.lead_lag_incidence';
|
||||
max_lead_lag=zeros(nendo,2); % lead/lag representation in Y for each endogenous variable i
|
||||
if ( M_.maximum_lag <= 1) && (M_.maximum_lead <= 1)
|
||||
xlen=nendo; %%=0;
|
||||
AA0=zeros(xlen,xlen); % empty A0
|
||||
AA2=AA0; % empty A2 and A3
|
||||
AA3=AA0;
|
||||
AA1=jacobia_(:,npred+1:npred+nendo);
|
||||
fnd = find(lead_lag(:,3));
|
||||
AA0(:, fnd)= jacobia_(:,nonzeros(lead_lag(:,3))); %forwd jacobian
|
||||
fnd = find(lead_lag(:,1));
|
||||
AA2(:, fnd)= jacobia_(:,nonzeros(lead_lag(:,1))); %backward
|
||||
|
||||
if M_.orig_endo_nbr<nendo
|
||||
exp_0= strmatch('AUX_EXPECT_LEAD_0_', M_.endo_names);
|
||||
num_exp_0=size(exp_0,1);
|
||||
if num_exp_0>0
|
||||
AA3(:,exp_0)=AA1(:,exp_0);
|
||||
XX0=zeros(nendo,num_exp_0);
|
||||
AA1(:,exp_0)=XX0(:,[1:num_exp_0])
|
||||
end
|
||||
end
|
||||
|
||||
else
|
||||
xlen=0;
|
||||
for i=1:nendo
|
||||
llmask=find(lead_lag(i,:)); % mask of leads and lags for var i
|
||||
nlag = max((M_.maximum_lag+1-min(llmask)), 0); % reduced no of lags and
|
||||
nlead = max((max(llmask)-(M_.maximum_lag+1)), 0); % reduced no of leads for var i
|
||||
max_lead_lag(i,:)=[nlag nlead]; % store for future reference
|
||||
%xlen=xlen+(nlag+nlead+1); % size as the sum over all the i variables
|
||||
xlen=xlen+(max(nlag-1,0)+max(nlead-1,0)+1); % size as the sum over all the i variables
|
||||
end
|
||||
AA0=zeros(xlen,xlen); % empty A0
|
||||
AA2=AA0; % empty A2 and A3
|
||||
AA3=AA0;
|
||||
end
|
||||
if (xlen>nendo )||( M_.maximum_lag >1) || (M_.maximum_lead >1) % we could not use shortcut above
|
||||
start=xlen -nendo+1;
|
||||
offset=0;
|
||||
for i=1:nendo
|
||||
llmask=find(lead_lag(i,:)); % mask of leads and lags for var i
|
||||
nlag=max_lead_lag(i,1); % size for the i'th variable lags
|
||||
nlead=max_lead_lag(i,2); % size for the i'th variable lead
|
||||
ilen=max(nlag-1,0)+max(nlead-1,0)+1; % size for the i'th variable
|
||||
if lead_lag(i,M_.maximum_lag-nlag+1) && nlag>0 %(j0==1 )&& lead_lag(i,j0) % !=0 %(ilen - iLagLen) %<=max(nlag,1)
|
||||
AA2( start:end, offset+1)=aa(:,lead_lag(i,M_.maximum_lag-nlag+1));
|
||||
else
|
||||
AA2( start+i-1, offset+1)=Inf;
|
||||
end
|
||||
if lead_lag(i,klen-M_.maximum_lead+nlead) && nlead>0 %(j0==1 )&& lead_lag(i,j0) % !=0 %(ilen - iLagLen) %<=max(nlag,1)
|
||||
AA0( start:end, offset+ilen)=aa(:,lead_lag(i,klen-M_.maximum_lead+nlead));
|
||||
else
|
||||
AA0( start+i-1, offset+ilen)=Inf;
|
||||
end
|
||||
for j0= 1:ilen % M_.maximum_lag
|
||||
if (j0<max(nlag-1,0)+1) && (ilen>nlag)&& lead_lag(i,j0+1) % !=0 %(ilen - iLagLen) %
|
||||
AA1( start:end, offset+j0)=aa(:,lead_lag(i,(j0+1)));
|
||||
elseif (j0==max(nlag-1,0)+1) && lead_lag(i,(M_.maximum_lag+1)) %&& (j0<=ilen-max(nlead-1,0) ) ...
|
||||
%AA1( start+i-1, offset+j0)=Inf; // jacobian at t
|
||||
AA1( start:end, offset+j0)=aa(:,lead_lag(i,(M_.maximum_lag+1)));
|
||||
elseif (j0>max(nlag-1,0)+1)&& (ilen>nlead)&& lead_lag(i,M_.maximum_lag+j0+1)
|
||||
AA1( start:end, offset+j0)=aa(:,lead_lag(i,(M_.maximum_lag+1+j0)));
|
||||
end
|
||||
end
|
||||
|
||||
offset=offset+ilen;
|
||||
if offset>xlen
|
||||
error(' dr1_PI: offset exceeds max xlen!');
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
PSI=-[[zeros(xlen-nendo,M_.exo_nbr)];[jacobia_(:, jlen-M_.exo_nbr+1:end)]]; % exog
|
||||
cc=0;
|
||||
NX=M_.exo_nbr; % no of exogenous varexo shock variables.
|
||||
NETA=nfwrd+nboth; % total no of exp. errors set to no of forward looking equations
|
||||
FL_RANK=rank(AA0); % nfwrd+nboth; % min total no of forward looking equations and vars
|
||||
|
||||
try
|
||||
|
||||
% call [G1pi,C,impact,nmat,TT1,TT2,gev,eu]=PI_gensys(a0,a1,a2,c,PSI,NX,NETA,NO_FL_EQS)
|
||||
% System given as
|
||||
% a0*E_t[y(t+1])+a1*y(t)=a2*y(t-1)+c+psi*eps(t)
|
||||
% with eps 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)']'.
|
||||
if(options_.ACES_solver==1)
|
||||
SSbar= diag(dr.ys);%(oo_.steady_state);
|
||||
AA0=AA0*SSbar;
|
||||
AA1=AA1*SSbar;
|
||||
AA2=AA2*SSbar;
|
||||
AA3=AA3*SSbar;
|
||||
end
|
||||
%%if any(AA3) % for expectational models when complete
|
||||
%% [G1pi,CC,impact,nmat,TT1,TT2,gev,eu, DD, E3,E5]=PI_gensysEXP(AA0,AA1,-AA2,AA3,cc,PSI,NX,NETA,FL_RANK, M_, options_);
|
||||
%%else
|
||||
[G1pi,CC,impact,nmat,TT1,TT2,gev,eu, DD, E3,E5]=PI_gensys(AA0,AA1,-AA2,AA3,cc,PSI,NX,NETA,FL_RANK, M_, options_);
|
||||
%%end
|
||||
|
||||
% reuse some of the bypassed code and tests that may be needed
|
||||
if eu ~=[1; 1]
|
||||
info(1) = abs(eu(1)+eu(2));
|
||||
info(2) = 1.0e+8;
|
||||
% return
|
||||
end
|
||||
|
||||
dr.PI_ghx=G1pi;
|
||||
dr.PI_ghu=impact;
|
||||
dr.PI_TT1=TT1;
|
||||
dr.PI_TT2=TT2;
|
||||
dr.PI_nmat=nmat;
|
||||
dr.PI_CC=CC;
|
||||
dr.PI_gev=gev;
|
||||
dr.PI_eu=eu;
|
||||
dr.PI_FL_RANK=FL_RANK;
|
||||
%dr.ys=zeros(nendo); % zero steady state
|
||||
dr.ghx=G1pi;
|
||||
dr.ghu=impact;
|
||||
dr.eigval = eig(G1pi);
|
||||
dr.rank=FL_RANK;
|
||||
|
||||
if options_.ACES_solver==1
|
||||
ACES.A=G1pi;
|
||||
ACES.C=impact; % (:,1);
|
||||
ACES.D=DD; %=impact (:,20);
|
||||
ACES.E3=E3;
|
||||
ACES.E5=E5;
|
||||
save ACES ;%ACES_A ACES_C ACES_D ACES_E2 ACES_E5
|
||||
%save([M_.fname '_ACES_IN'], 'ACES')
|
||||
save ([M_.fname '_ACES_A.txt'], 'G1pi', '-ascii', '-double', '-tabs');
|
||||
save ([M_.fname '_ACES_C.txt'], 'impact','-ascii', '-double', '-tabs');
|
||||
save ([M_.fname '_ACES_D.txt'], 'DD', '-ascii', '-double', '-tabs');
|
||||
save ([M_.fname '_ACES_E3.txt'], 'E3', '-ascii', '-double', '-tabs');
|
||||
save ([M_.fname '_ACES_E5.txt'], 'E5', '-ascii', '-double', '-tabs');
|
||||
end
|
||||
|
||||
catch
|
||||
disp('Problem with using Part Info solver - Using Dynare solver instead');
|
||||
lerror=lasterror;
|
||||
disp (lerror);
|
||||
options_.PartInfo = 0; % and then try mjdgges instead
|
||||
end
|
||||
|
||||
% TODO:
|
||||
% if options_.loglinear == 1
|
||||
% if exogenous deterministic variables
|
||||
|
||||
return;
|
|
@ -0,0 +1,146 @@
|
|||
function [dr,info]=resol(ys,check_flag)
|
||||
% function [dr,info]=resol(ys,check_flag)
|
||||
% Computes first and second order approximations
|
||||
%
|
||||
% INPUTS
|
||||
% ys: vector of variables in steady state
|
||||
% check_flag=0: all the approximation is computed
|
||||
% check_flag=1: computes only the eigenvalues
|
||||
%
|
||||
% OUTPUTS
|
||||
% dr: structure of decision rules for stochastic simulations
|
||||
% info=1: the model doesn't determine the current variables '...' uniquely
|
||||
% info=2: MJDGGES returns the following error code'
|
||||
% info=3: Blanchard Kahn conditions are not satisfied: no stable '...' equilibrium
|
||||
% info=4: Blanchard Kahn conditions are not satisfied:'...' indeterminacy
|
||||
% info=5: Blanchard Kahn conditions are not satisfied:'...' indeterminacy due to rank failure
|
||||
% info=6: The jacobian evaluated at the steady state is complex.
|
||||
% info=19: The steadystate file did not compute the steady state (inconsistent deep parameters).
|
||||
% info=20: can't find steady state info(2) contains sum of sqare residuals
|
||||
% info=21: steady state is complex
|
||||
% info(2) contains sum of sqare of
|
||||
% imaginary part of steady state
|
||||
% info=30: Variance can't be computed
|
||||
%
|
||||
% SPECIAL REQUIREMENTS
|
||||
% none
|
||||
|
||||
% Copyright (C) 2001-2009 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 <http://www.gnu.org/licenses/>.
|
||||
|
||||
global M_ options_ oo_
|
||||
global it_
|
||||
|
||||
jacobian_flag = 0;
|
||||
|
||||
options_ = set_default_option(options_,'jacobian_flag',1);
|
||||
info = 0;
|
||||
|
||||
it_ = M_.maximum_lag + 1 ;
|
||||
|
||||
if M_.exo_nbr == 0
|
||||
oo_.exo_steady_state = [] ;
|
||||
end
|
||||
|
||||
% check if ys is steady state
|
||||
tempex = oo_.exo_simul;
|
||||
oo_.exo_simul = repmat(oo_.exo_steady_state',M_.maximum_lag+M_.maximum_lead+1,1);
|
||||
if M_.exo_det_nbr > 0
|
||||
tempexdet = oo_.exo_det_simul;
|
||||
oo_.exo_det_simul = repmat(oo_.exo_det_steady_state',M_.maximum_lag+M_.maximum_lead+1,1);
|
||||
end
|
||||
dr.ys = ys;
|
||||
check1 = 0;
|
||||
% testing for steadystate file
|
||||
fh = str2func([M_.fname '_static']);
|
||||
if options_.steadystate_flag
|
||||
[dr.ys,check1] = feval([M_.fname '_steadystate'],dr.ys,...
|
||||
[oo_.exo_steady_state; ...
|
||||
oo_.exo_det_steady_state]);
|
||||
if size(dr.ys,1) < M_.endo_nbr
|
||||
if isfield(M_,'aux_vars')
|
||||
dr.ys = add_auxiliary_variables_to_steadystate(dr.ys,M_.aux_vars,...
|
||||
M_.fname,...
|
||||
oo_.exo_steady_state,...
|
||||
oo_.exo_det_steady_state,...
|
||||
M_.params);
|
||||
else
|
||||
error([M_.fname '_steadystate.m doesn''t match the model']);
|
||||
end
|
||||
end
|
||||
|
||||
else
|
||||
% testing if ys isn't a steady state or if we aren't computing Ramsey policy
|
||||
if options_.ramsey_policy == 0
|
||||
if options_.linear == 0
|
||||
% nonlinear models
|
||||
if max(abs(feval(fh,dr.ys,[oo_.exo_steady_state; ...
|
||||
oo_.exo_det_steady_state], M_.params))) > options_.dynatol
|
||||
[dr.ys,check1] = dynare_solve(fh,dr.ys,options_.jacobian_flag,...
|
||||
[oo_.exo_steady_state; ...
|
||||
oo_.exo_det_steady_state], M_.params);
|
||||
end
|
||||
else
|
||||
% linear models
|
||||
[fvec,jacob] = feval(fh,dr.ys,[oo_.exo_steady_state;...
|
||||
oo_.exo_det_steady_state], M_.params);
|
||||
if max(abs(fvec)) > 1e-12
|
||||
dr.ys = dr.ys-jacob\fvec;
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
% testing for problem
|
||||
if check1
|
||||
if options_.steadystate_flag
|
||||
info(1)= 19;
|
||||
resid = check1 ;
|
||||
else
|
||||
info(1)= 20;
|
||||
resid = feval(fh,ys,oo_.exo_steady_state, M_.params);
|
||||
end
|
||||
info(2) = resid'*resid ;
|
||||
return
|
||||
end
|
||||
|
||||
if ~isreal(dr.ys)
|
||||
info(1) = 21;
|
||||
info(2) = sum(imag(ys).^2);
|
||||
dr.ys = real(dr.ys);
|
||||
return
|
||||
end
|
||||
|
||||
dr.fbias = zeros(M_.endo_nbr,1);
|
||||
if( (options_.partial_information ==1) | (options_.ACES_solver==1))%&& (check_flag == 0)
|
||||
[dr,info,M_,options_,oo_] = dr1_PI(dr,check_flag,M_,options_,oo_);
|
||||
else
|
||||
[dr,info,M_,options_,oo_] = dr1(dr,check_flag,M_,options_,oo_);
|
||||
end
|
||||
if info(1)
|
||||
return
|
||||
end
|
||||
|
||||
if M_.exo_det_nbr > 0
|
||||
oo_.exo_det_simul = tempexdet;
|
||||
end
|
||||
oo_.exo_simul = tempex;
|
||||
tempex = [];
|
||||
|
||||
% 01/01/2003 MJ added dr_algo == 1
|
||||
% 08/24/2001 MJ uses Schmitt-Grohe and Uribe (2001) constant correction
|
||||
% in dr.ghs2
|
||||
% 05/26/2003 MJ added temporary values for oo_.exo_simul
|
|
@ -0,0 +1,327 @@
|
|||
function info=stoch_simul(var_list)
|
||||
|
||||
% 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 <http://www.gnu.org/licenses/>.
|
||||
|
||||
global M_ options_ oo_ it_
|
||||
|
||||
options_old = options_;
|
||||
if options_.linear
|
||||
options_.order = 1;
|
||||
end
|
||||
if (options_.order == 1)
|
||||
options_.replic = 1;
|
||||
elseif options_.order == 3
|
||||
options_.k_order_solver = 1;
|
||||
end
|
||||
|
||||
|
||||
TeX = options_.TeX;
|
||||
|
||||
iter_ = max(options_.periods,1);
|
||||
if M_.exo_nbr > 0
|
||||
oo_.exo_simul= ones(iter_ + M_.maximum_lag + M_.maximum_lead,1) * oo_.exo_steady_state';
|
||||
end
|
||||
|
||||
check_model;
|
||||
|
||||
[oo_.dr, info] = resol(oo_.steady_state,0);
|
||||
|
||||
if info(1)
|
||||
options_ = options_old;
|
||||
print_info(info, options_.noprint);
|
||||
return
|
||||
end
|
||||
|
||||
if ~options_.noprint
|
||||
disp(' ')
|
||||
disp('MODEL SUMMARY')
|
||||
disp(' ')
|
||||
disp([' Number of variables: ' int2str(M_.endo_nbr)])
|
||||
disp([' Number of stochastic shocks: ' int2str(M_.exo_nbr)])
|
||||
disp([' Number of state variables: ' ...
|
||||
int2str(length(find(oo_.dr.kstate(:,2) <= M_.maximum_lag+1)))])
|
||||
disp([' Number of jumpers: ' ...
|
||||
int2str(length(find(oo_.dr.kstate(:,2) == M_.maximum_lag+2)))])
|
||||
disp([' Number of static variables: ' int2str(oo_.dr.nstatic)])
|
||||
my_title='MATRIX OF COVARIANCE OF EXOGENOUS SHOCKS';
|
||||
labels = deblank(M_.exo_names);
|
||||
headers = strvcat('Variables',labels);
|
||||
lh = size(labels,2)+2;
|
||||
dyntable(my_title,headers,labels,M_.Sigma_e,lh,10,6);
|
||||
disp(' ')
|
||||
if (options_.order <= 2 & options_.partial_information ==0 & options_.ACES_solver==0)
|
||||
disp_dr(oo_.dr,options_.order,var_list);
|
||||
end
|
||||
end
|
||||
|
||||
if options_.irf
|
||||
if size(var_list,1) == 0
|
||||
var_list = M_.endo_names(1:M_.orig_endo_nbr, :);
|
||||
if TeX
|
||||
var_listTeX = M_.endo_names_tex(1:M_.orig_endo_nbr, :);
|
||||
end
|
||||
end;%else
|
||||
n = size(var_list,1);
|
||||
ivar=zeros(n,1);
|
||||
if TeX
|
||||
var_listTeX = [];
|
||||
end
|
||||
for i=1:n
|
||||
i_tmp = strmatch(var_list(i,:),M_.endo_names,'exact');
|
||||
if isempty(i_tmp)
|
||||
error (['One of the specified variables does not exist']) ;
|
||||
else
|
||||
ivar(i) = i_tmp;
|
||||
if TeX
|
||||
var_listTeX = strvcat(var_listTeX,deblank(M_.endo_names_tex(i_tmp,:)));
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
if ( options_.partial_information ==0 & options_.ACES_solver==0)
|
||||
if options_.periods == 0 && options_.nomoments == 0
|
||||
disp_th_moments(oo_.dr,var_list);
|
||||
elseif options_.periods ~= 0
|
||||
if options_.periods < options_.drop
|
||||
disp(['STOCH_SIMUL error: The horizon of simulation is shorter' ...
|
||||
' than the number of observations to be DROPed'])
|
||||
options_ =options_old;
|
||||
return
|
||||
end
|
||||
oo_.endo_simul = simult(repmat(oo_.dr.ys,1,M_.maximum_lag),oo_.dr);
|
||||
dyn2vec;
|
||||
if options_.nomoments == 0
|
||||
disp_moments(oo_.endo_simul,var_list);
|
||||
end
|
||||
end
|
||||
else % part info or ACES
|
||||
if ~isfield(options_,'varobs')
|
||||
nobs =0;
|
||||
else
|
||||
if isempty(options_.varobs)
|
||||
nobs =0;
|
||||
else
|
||||
varobs=options_.varobs;
|
||||
nobs = size(varobs,1);
|
||||
end
|
||||
end
|
||||
if nobs == 0
|
||||
ivarobs=ivar;
|
||||
else
|
||||
ivarobs=zeros(nobs,1);
|
||||
if TeX
|
||||
varobs_listTeX = [];
|
||||
end
|
||||
for i=1:nobs
|
||||
i_tmp = strmatch(varobs(i,:),M_.endo_names,'exact');
|
||||
if isempty(i_tmp)
|
||||
error (['One of the specified observed variables does not exist']) ;
|
||||
else
|
||||
ivarobs(i) = i_tmp;
|
||||
if TeX
|
||||
varobs_listTeX = strvcat(varobs_listTeX,deblank(M_.endo_names_tex(i_tmp,:)));
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
PCL_Part_info_moments (0, ivarobs', oo_.dr, ivar, var_list);
|
||||
end
|
||||
|
||||
if options_.irf
|
||||
|
||||
if TeX
|
||||
fidTeX = fopen([M_.fname '_IRF.TeX'],'w');
|
||||
fprintf(fidTeX,'%% TeX eps-loader file generated by stoch_simul.m (Dynare).\n');
|
||||
fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
|
||||
fprintf(fidTeX,' \n');
|
||||
end
|
||||
olditer = iter_;% Est-ce vraiment utile ? Il y a la même ligne dans irf...
|
||||
SS(M_.exo_names_orig_ord,M_.exo_names_orig_ord)=M_.Sigma_e+1e-14*eye(M_.exo_nbr);
|
||||
cs = transpose(chol(SS));
|
||||
tit(M_.exo_names_orig_ord,:) = M_.exo_names;
|
||||
if TeX
|
||||
titTeX(M_.exo_names_orig_ord,:) = M_.exo_names_tex;
|
||||
end
|
||||
for i=1:M_.exo_nbr
|
||||
if SS(i,i) > 1e-13
|
||||
if ( options_.partial_information ==1 | options_.ACES_solver==1)
|
||||
irfmat=PCL_Part_info_irf (0, ivarobs', M_, oo_.dr, options_.irf, i);
|
||||
y=zeros(M_.endo_nbr,options_.irf);
|
||||
y(ivar,:)=irfmat(ivar,1:options_.irf);
|
||||
else
|
||||
y=irf(oo_.dr,cs(M_.exo_names_orig_ord,i), options_.irf, options_.drop, ...
|
||||
options_.replic, options_.order);
|
||||
end
|
||||
if options_.relative_irf
|
||||
y = 100*y/cs(i,i);
|
||||
end
|
||||
irfs = [];
|
||||
mylist = [];
|
||||
if TeX
|
||||
mylistTeX = [];
|
||||
end
|
||||
for j = 1:n
|
||||
assignin('base',[deblank(M_.endo_names(ivar(j),:)) '_' deblank(M_.exo_names(i,:))],...
|
||||
y(ivar(j),:)');
|
||||
eval(['oo_.irfs.' deblank(M_.endo_names(ivar(j),:)) '_' ...
|
||||
deblank(M_.exo_names(i,:)) ' = y(ivar(j),:);']);
|
||||
if max(y(ivar(j),:)) - min(y(ivar(j),:)) > 1e-10
|
||||
irfs = cat(1,irfs,y(ivar(j),:));
|
||||
mylist = strvcat(mylist,deblank(var_list(j,:)));
|
||||
if TeX
|
||||
mylistTeX = strvcat(mylistTeX,deblank(var_listTeX(j,:)));
|
||||
end
|
||||
end
|
||||
end
|
||||
if options_.nograph == 0
|
||||
number_of_plots_to_draw = size(irfs,1);
|
||||
[nbplt,nr,nc,lr,lc,nstar] = pltorg(number_of_plots_to_draw);
|
||||
if nbplt == 0
|
||||
elseif nbplt == 1
|
||||
if options_.relative_irf
|
||||
hh = figure('Name',['Relative response to' ...
|
||||
' orthogonalized shock to ' tit(i,:)]);
|
||||
else
|
||||
hh = figure('Name',['Orthogonalized shock to' ...
|
||||
' ' tit(i,:)]);
|
||||
end
|
||||
for j = 1:number_of_plots_to_draw
|
||||
subplot(nr,nc,j);
|
||||
plot(1:options_.irf,transpose(irfs(j,:)),'-k','linewidth',1);
|
||||
hold on
|
||||
plot([1 options_.irf],[0 0],'-r','linewidth',0.5);
|
||||
hold off
|
||||
xlim([1 options_.irf]);
|
||||
title(deblank(mylist(j,:)),'Interpreter','none');
|
||||
end
|
||||
eval(['print -depsc2 ' M_.fname '_IRF_' deblank(tit(i,:)) '.eps']);
|
||||
if ~exist('OCTAVE_VERSION')
|
||||
eval(['print -dpdf ' M_.fname '_IRF_' deblank(tit(i,:))]);
|
||||
saveas(hh,[M_.fname '_IRF_' deblank(tit(i,:)) '.fig']);
|
||||
end
|
||||
if TeX
|
||||
fprintf(fidTeX,'\\begin{figure}[H]\n');
|
||||
for j = 1:number_of_plots_to_draw
|
||||
fprintf(fidTeX,['\\psfrag{%s}[1][][0.5][0]{$%s$}\n'],deblank(mylist(j,:)),deblank(mylistTeX(j,:)));
|
||||
end
|
||||
fprintf(fidTeX,'\\centering \n');
|
||||
fprintf(fidTeX,'\\includegraphics[scale=0.5]{%s_IRF_%s}\n',M_.fname,deblank(tit(i,:)));
|
||||
fprintf(fidTeX,'\\caption{Impulse response functions (orthogonalized shock to $%s$).}',titTeX(i,:));
|
||||
fprintf(fidTeX,'\\label{Fig:IRF:%s}\n',deblank(tit(i,:)));
|
||||
fprintf(fidTeX,'\\end{figure}\n');
|
||||
fprintf(fidTeX,' \n');
|
||||
end
|
||||
% close(hh)
|
||||
else
|
||||
for fig = 1:nbplt-1
|
||||
if options_.relative_irf == 1
|
||||
hh = figure('Name',['Relative response to orthogonalized shock' ...
|
||||
' to ' tit(i,:) ' figure ' int2str(fig)]);
|
||||
else
|
||||
hh = figure('Name',['Orthogonalized shock to ' tit(i,:) ...
|
||||
' figure ' int2str(fig)]);
|
||||
end
|
||||
for plt = 1:nstar
|
||||
subplot(nr,nc,plt);
|
||||
plot(1:options_.irf,transpose(irfs((fig-1)*nstar+plt,:)),'-k','linewidth',1);
|
||||
hold on
|
||||
plot([1 options_.irf],[0 0],'-r','linewidth',0.5);
|
||||
hold off
|
||||
xlim([1 options_.irf]);
|
||||
title(deblank(mylist((fig-1)*nstar+plt,:)),'Interpreter','none');
|
||||
end
|
||||
eval(['print -depsc2 ' M_.fname '_IRF_' deblank(tit(i,:)) int2str(fig) '.eps']);
|
||||
if ~exist('OCTAVE_VERSION')
|
||||
eval(['print -dpdf ' M_.fname '_IRF_' deblank(tit(i,:)) int2str(fig)]);
|
||||
saveas(hh,[M_.fname '_IRF_' deblank(tit(i,:)) int2str(fig) '.fig']);
|
||||
end
|
||||
if TeX
|
||||
fprintf(fidTeX,'\\begin{figure}[H]\n');
|
||||
for j = 1:nstar
|
||||
fprintf(fidTeX,['\\psfrag{%s}[1][][0.5][0]{$%s$}\n'],deblank(mylist((fig-1)*nstar+j,:)),deblank(mylistTeX((fig-1)*nstar+j,:)));
|
||||
end
|
||||
fprintf(fidTeX,'\\centering \n');
|
||||
fprintf(fidTeX,'\\includegraphics[scale=0.5]{%s_IRF_%s%s}\n',M_.fname,deblank(tit(i,:)),int2str(fig));
|
||||
if options_.relative_irf
|
||||
fprintf(fidTeX,['\\caption{Relative impulse response' ...
|
||||
' functions (orthogonalized shock to $%s$).}'],deblank(titTeX(i,:)));
|
||||
else
|
||||
fprintf(fidTeX,['\\caption{Impulse response functions' ...
|
||||
' (orthogonalized shock to $%s$).}'],deblank(titTeX(i,:)));
|
||||
end
|
||||
fprintf(fidTeX,'\\label{Fig:BayesianIRF:%s:%s}\n',deblank(tit(i,:)),int2str(fig));
|
||||
fprintf(fidTeX,'\\end{figure}\n');
|
||||
fprintf(fidTeX,' \n');
|
||||
end
|
||||
% close(hh);
|
||||
end
|
||||
hh = figure('Name',['Orthogonalized shock to ' tit(i,:) ' figure ' int2str(nbplt) '.']);
|
||||
m = 0;
|
||||
for plt = 1:number_of_plots_to_draw-(nbplt-1)*nstar;
|
||||
m = m+1;
|
||||
subplot(lr,lc,m);
|
||||
plot(1:options_.irf,transpose(irfs((nbplt-1)*nstar+plt,:)),'-k','linewidth',1);
|
||||
hold on
|
||||
plot([1 options_.irf],[0 0],'-r','linewidth',0.5);
|
||||
hold off
|
||||
xlim([1 options_.irf]);
|
||||
title(deblank(mylist((nbplt-1)*nstar+plt,:)),'Interpreter','none');
|
||||
end
|
||||
eval(['print -depsc2 ' M_.fname '_IRF_' deblank(tit(i,:)) int2str(nbplt) '.eps']);
|
||||
if ~exist('OCTAVE_VERSION')
|
||||
eval(['print -dpdf ' M_.fname '_IRF_' deblank(tit(i,:)) int2str(nbplt)]);
|
||||
saveas(hh,[M_.fname '_IRF_' deblank(tit(i,:)) int2str(nbplt) '.fig']);
|
||||
end
|
||||
if TeX
|
||||
fprintf(fidTeX,'\\begin{figure}[H]\n');
|
||||
for j = 1:m
|
||||
fprintf(fidTeX,['\\psfrag{%s}[1][][0.5][0]{$%s$}\n'],deblank(mylist((nbplt-1)*nstar+j,:)),deblank(mylistTeX((nbplt-1)*nstar+j,:)));
|
||||
end
|
||||
fprintf(fidTeX,'\\centering \n');
|
||||
fprintf(fidTeX,'\\includegraphics[scale=0.5]{%s_IRF_%s%s}\n',M_.fname,deblank(tit(i,:)),int2str(nbplt));
|
||||
if options_.relative_irf
|
||||
fprintf(fidTeX,['\\caption{Relative impulse response functions' ...
|
||||
' (orthogonalized shock to $%s$).}'],deblank(titTeX(i,:)));
|
||||
else
|
||||
fprintf(fidTeX,['\\caption{Impulse response functions' ...
|
||||
' (orthogonalized shock to $%s$).}'],deblank(titTeX(i,:)));
|
||||
end
|
||||
fprintf(fidTeX,'\\label{Fig:IRF:%s:%s}\n',deblank(tit(i,:)),int2str(nbplt));
|
||||
fprintf(fidTeX,'\\end{figure}\n');
|
||||
fprintf(fidTeX,' \n');
|
||||
end
|
||||
% close(hh);
|
||||
end
|
||||
end
|
||||
end
|
||||
iter_ = olditer;
|
||||
if TeX
|
||||
fprintf(fidTeX,' \n');
|
||||
fprintf(fidTeX,'%% End Of TeX file. \n');
|
||||
fclose(fidTeX);
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
if options_.SpectralDensity == 1
|
||||
[omega,f] = UnivariateSpectralDensity(oo_.dr,var_list);
|
||||
end
|
||||
|
||||
|
||||
options_ = options_old;
|
|
@ -0,0 +1,55 @@
|
|||
% price markup shock (eps_m) AND inflation target (pitarg) modelled as AR1
|
||||
% estimating varrho
|
||||
% cy = 0.614479/0.769365 - obtained from the computed steady state from the nonlinear counterpart
|
||||
% use usdata1.mat
|
||||
% this version - 13/07/09
|
||||
|
||||
|
||||
var pi mc mun muc c y n r g a pitarg PI Y RR;
|
||||
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;//0.5575;//0.64;
|
||||
hc = 0.0;//0.7984;//0.86;
|
||||
wd = 0.40;
|
||||
gamma = 0.5868;//0.5792;//0.00;
|
||||
sigma = 4.0897;//2.0034;//3.47
|
||||
cy = 0.614479/0.769365; // cy = Cbar/Ybar
|
||||
rho_g = 0.8325;//0.9132;//0.85;
|
||||
rho_a = 0.9827;//0.9787;//0.87;
|
||||
rho_r = 0.3529;//0.4945;//0.76;
|
||||
thetap = 2.2161;//2.2062;//2.30;
|
||||
varrho = 0.3853;//0.4432;//0.66;
|
||||
rho_targ = 0.6133;//0.5807;//0.9;
|
||||
|
||||
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;
|
||||
PI = pi;
|
||||
Y = y;
|
||||
RR = r;
|
||||
end;
|
||||
|
||||
shocks;
|
||||
var eps_g; stderr 3.8505;//4.1018;//0.325;
|
||||
var eps_a; stderr 0.7573;//0.9725;//0.598;
|
||||
var eps_e; stderr 0.2409;//0.2466;//0.150;
|
||||
var eps_m; stderr 0.8329;//0.6487;//0.2;
|
||||
var eps_targ; stderr 0.3978;//0.4577;//0.03;
|
||||
end;
|
||||
|
||||
//steady;
|
||||
//check;
|
||||
varobs pi mc mun muc c y n r g a pitarg;
|
||||
stoch_simul(linear,partial_information,periods=1000,irf=30)pi y r;//pi n c y r;
|
||||
|
|
@ -0,0 +1,55 @@
|
|||
% price markup shock (eps_m) AND inflation target (pitarg) modelled as AR1
|
||||
% estimating varrho
|
||||
% cy = 0.614479/0.769365 - obtained from the computed steady state from the nonlinear counterpart
|
||||
% use usdata1.mat
|
||||
% this version - 13/07/09
|
||||
|
||||
|
||||
var pi mc mun muc c y n r g a pitarg PI Y RR;
|
||||
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;//0.5575;//0.64;
|
||||
hc = 0.0;//0.7984;//0.86;
|
||||
wd = 0.40;
|
||||
gamma = 0.5868;//0.5792;//0.00;
|
||||
sigma = 4.0897;//2.0034;//3.47
|
||||
cy = 0.614479/0.769365; // cy = Cbar/Ybar
|
||||
rho_g = 0.8325;//0.9132;//0.85;
|
||||
rho_a = 0.9827;//0.9787;//0.87;
|
||||
rho_r = 0.3529;//0.4945;//0.76;
|
||||
thetap = 2.2161;//2.2062;//2.30;
|
||||
varrho = 0.3853;//0.4432;//0.66;
|
||||
rho_targ = 0.6133;//0.5807;//0.9;
|
||||
|
||||
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;
|
||||
PI = pi;
|
||||
Y = y;
|
||||
RR = r;
|
||||
end;
|
||||
|
||||
shocks;
|
||||
var eps_g; stderr 3.8505;//4.1018;//0.325;
|
||||
var eps_a; stderr 0.7573;//0.9725;//0.598;
|
||||
var eps_e; stderr 0.2409;//0.2466;//0.150;
|
||||
var eps_m; stderr 0.8329;//0.6487;//0.2;
|
||||
var eps_targ; stderr 0.3978;//0.4577;//0.03;
|
||||
end;
|
||||
|
||||
//steady;
|
||||
//check;
|
||||
varobs c n r ;
|
||||
stoch_simul(linear,partial_information,periods=1000,irf=30)pi y r;//pi n c y r;
|
||||
|
Loading…
Reference in New Issue