%% ZLB code for Smets-Wouters model
% This code computes ZLB impulse response function based on the Smets-Wouters
% (2007) model solved in dynare. To compute the model write the model name
% into the cell models_cell. To implement the ZLB you will neeed to tell the 
% file what the name of the nominal interest rate is in the model (it will 
% find the location automatically) -- save this as the variable "iname".
% 
% Note: The original Smets-Wouters code file has been downloaded from AER
% website and modified to run with Dynare 4.2.2. The
% estimates are identical up to the second digit to those reported in
% Smets-Wouters (2007). You will need to download and install Dynare from
% dynare.org.

%% Running dynare code
models_cell =   {'usmodel_final'};
iname       =   'r';
irflen      =   40;     % length of IRFs from Dynare
irflen2     =   40;     % length of IRFs computed in this code
save model_info_final models_cell iname irflen irflen2

%% Set-up for model implementation
eval(['dynare ' models_cell{1} '.mod'])
save model_all_final


% loading output from dynare
clear all
close all
load model_info_final
load model_all_final


%% Equation set-up
% Equation setup is 
% 
% $$ AE_{t}x_{t+1}+Bx_{t}+Cx_{t-1}+D\epsilon_{t} $$ 
% 
% where $\epsilon_{t}$ are the shocks. 

% Obtaining linarized solution from dynare
nstate  =   oo_.dr.npred;           % number of states
nshock  =   M_.exo_nbr;             % number of shocks
nlead   =   oo_.dr.nsfwrd;          % number of jumps
ndyn    =   nstate+nlead;           
nvar    =   M_.endo_nbr;            % number of variables
sstate  =   oo_.steady_state;       % steady state
dsstate =   [];
for jj  =   1:3
    eval(['ssloc' num2str(jj)  ' =   find(M_.lead_lag_incidence(jj,:)>0);'])
    eval(['dsstate =   [dsstate;sstate(ssloc' num2str(jj) ')];'])
%     eval(['setzero' num2str(jj) '=  setdiff(1:1:nvar,ssloc)  ;']);
end

% dsstate =   [zeros(nlead+1:ndyn,1);sstate;zeros(1:nlead,1)];   % dynamic steady state (includes lags and leads)
eval(['[stuff,jacobia_2] =   '  models_cell{1}  '_dynamic(dsstate, zeros(1,nshock), M_.params, oo_.steady_state,1);']);  % arguments: dyn ss (lags,contemp,leads), shock values, parameters, shock value selection

% cashing matrices
A   =   zeros(nvar,nvar);
B   =   zeros(nvar,nvar);
C   =   zeros(nvar,nvar);
D   =   zeros(nvar,nshock);

% distributing jacobian matrix
C(:,ssloc1)     =   jacobia_2(:,1:nstate);
B(:,ssloc2)     =   jacobia_2(:,nstate+1:nstate+nvar);
A(:,ssloc3)     =   jacobia_2(:,nstate+nvar+1:end-nshock);
D(:,:)          =   jacobia_2(:,end-nshock+1:end);

%%  Solution to the system
% 
% $$ x_{t} = P x_{t-1}+Q \epsilon_{t} $$ 
% 
% Model is then solved using shooting as in Erceg et al (2009) and CGW
% (2012). The matrices P and Q are taking from Dynare and put into the
% correct format.

P       =   zeros(length(oo_.steady_state)); % cashed P matrix
P1      =   oo_.dr.ghx;                      % partial matrix of P with state variables in cols
P1ord   =   oo_.dr.inv_order_var;            % order of P1 matrix
P1col   =   oo_.dr.kstate(find(oo_.dr.kstate(:,4)>0),1);    % state variables



% creating P matrix from dynare objects
P(:,P1col)  =   P1;
P   =   P(P1ord,:);
P   =   P(:,P1ord);

Q   =   oo_.dr.ghu;             % full matrix Q
Q   =   Q(P1ord,:);

%% Variables of interest: positions in the matrices
cpos    =   find(strcmp(cellstr(M_.endo_names),'c')==1);
ypos    =   find(strcmp(cellstr(M_.endo_names),'y')==1);
rpos    =   find(strcmp(cellstr(M_.endo_names),'r')==1);
pinfpos	=   find(strcmp(cellstr(M_.endo_names),'pinf')==1);
mpos    =   find(strcmp(cellstr(M_.exo_names),'em')==1);
bpos    =   find(strcmp(cellstr(M_.exo_names),'eb')==1);

%% Test IRF (1): Monetary policy shock in linearized model
shock       =   zeros(nshock,1);
shock(mpos)	=   -sqrt(M_.Sigma_e(mpos,mpos))*0.01;
% shock(2)    =   -oo_.posterior_mode.shocks_std.eb*0.01;
irf         =   Q*shock;
for jj  =   1:irflen2
    irf =   [irf P*irf(:,end)];
end



shock       =   zeros(nshock,1);
shock(mpos)	=   -sqrt(M_.Sigma_e(mpos,mpos))*0.01;
shock(bpos)	=   -sqrt(M_.Sigma_e(bpos,bpos))*0.01;
Dirf        =   Q*shock;
for jj  =   1:irflen2
    Dirf=   [Dirf P*Dirf(:,end)];
end

irf2    =   irf+Dirf;
irf3    =   irf+2*Dirf;
irf4    =   irf+4*Dirf;

% these two should be identical: IRF from dynare and manual IRF from code
% without ZLB.
plotset={'y','c','r','pinf'};
figure(10)
for jj=1:length(plotset)
    eval(['subplot(2,2,' num2str(jj) '), plot(0:1:irflen2,irf(' plotset{jj} 'pos,:)*100,''b-'',0:1:irflen-1,-oo_.irfs.' plotset{jj} '_em(1:irflen),''r--'');'])
    title(plotset{jj});
    legend('If not identical = error!')
    hold on;
end
hold off;


% compute dY/dr;
% peak output response and initial output response;
dy1(1)	=   max(exp(irf(ypos,:))-1);
dy1(2)  =   exp(irf(ypos,1))-1;
dc1(1)	=   max(exp(irf(cpos,:))-1);
dc1(2)  =   exp(irf(cpos,1))-1;
dr1     =   prod(exp(irf(rpos,1:irflen2)-irf(pinfpos,1:irflen2)))^(1/((irflen2)/4))-1;

dydr1   =   dy1./dr1;
dcdr1   =   dc1./dr1;

display(['Normal Times: dy/dr is ' num2str(dydr1(1),2) ' at peak and ' num2str(dydr1(2),2) ' at inception.'])
display(['Normal Times: dc/dr is ' num2str(dcdr1(1),2) ' at peak and ' num2str(dcdr1(2),2) ' at inception.'])



%% ZLB Implementation
% When zero lower bound binds system is described by
% 
% $$ AS*E_t(s_t+1)+BS*s_t+CS*s_t-1+DS*e_t+d=0 $$
% 
% where the response of interest rates to any variable is zero

% determining location of interest rate and interest rate rule
bound       =   find(B(:,rpos)==1);
if length(bound)>1
    bound       =   intersect(bound,find(abs(B(:,find(strcmp(cellstr(M_.endo_names),'ms')==1)))==1)); % using monetary policy shock location if we cannot uniquely determine location yet.
end

% steady-state nominal rate;
cgammabar   =   1+oo_.posterior_mode.parameters.ctrend/100;
csigmacbar  =   oo_.posterior_mode.parameters.csigma;
cpistarbar  =   1+oo_.posterior_mode.parameters.constepinf/100;
cbettabar   =   (1+oo_.posterior_mode.parameters.constebeta/100)^-1;
rb          =   -(cbettabar^-1*cgammabar^csigmacbar*cpistarbar-1);
shocknames  =   cellstr(M_.exo_names);

% creating matrices for zero bound system
BS              =   B;                       
BS(bound,:)     =   zeros(1,nvar);       
BS(bound,rpos)  =   B(bound,rpos);
AS              =   A;                       
AS(bound,:)     =   zeros(1,nvar);
CS              =   C;                       
CS(bound,:)     =   zeros(1,nvar);
DS              =   D;                       
DS(bound,:)     =   zeros(1,nshock);

%and interest rates equal the lower bound
d               =   zeros(nvar,1);            
d(bound,1)      =   -rb*BS(bound,rpos);     %0.0145            

%calculating the necessary inverses
I=eye(nvar);
AB=-BS\AS;      
CB=-BS\CS;
db=-BS\d;
DB=-BS\DS;

ABb=-B\A;      
CBb=-B\C;
dbb=zeros(nvar,1);
DBb=-B\D;

% Calulcating impulse response functions taking into account ZLB.
% Algorithm is as described in Erceg et al (2009) and Coibion et al (2012).
shockval=-0.18;    % value for shock to get to ZLB.
for zz  =   1:2
    
    shock   =   zeros(nshock,1);
    
    if zz   ==   1
        shock(2)        =   shockval;
    elseif zz   ==  2
        shock(2)        =   shockval;
        shock(mpos)     =   shock(mpos)-600*oo_.posterior_mode.shocks_std.em;
    end

    i=1;
    checklb=0;
    checkub=0;
    checkblb=0;
    checkbub=0;
    T   =   0;

    while checklb==0 || checkub==0

        if T==0

            s=zeros(nvar,T+2);
            s(:,1)=Q*shock;  % P*x(:,i-1)+ %when the liquidity trap is not expected to bind, it follows the normal laws of motion 
            s(:,2)=P*s(:,1);

            checklb=1;                   %cannot move the liquidity trap length to T=-1
            checkub=(rb<s(rpos,1));   %checks that the zero lower bound does not bind today

        else

            G=zeros(nvar,nvar*T);       %Cashing variables
            h=zeros(nvar,T);            
            s=zeros(nvar,T+2);

            %Zero lower bound binds up to T. For T+1 onwards, system
            %follows s_t=P*s_t-1+Q*e_t, which is the AIM solution. Substituting the
            %solution for s_T+1 (noting e_T=0) into the equation
            %AS*E_t(s_T+1)+BS*s_T+CS*s_T-1+d=0
            %yields the solution s_T=G(1)*s_T-1+h(1) which give us the first set of
            %recursive coefficients [G(1),h(1)]:

            G(:,1:nvar)=-inv(AS*P+BS)*CS;
            h(:,1)=-inv(AS*P+BS)*d;

            %For all other periods t<T, the laws of motion are given by
            %s_t=BS^-1*AS*s_t+1+BS^-1*CS*s_t-1+BS^-1*DS*s_t-2+BS^-1*d
            %s_1=BS^-1*AS*s_2+BS^-1*CS*s_0+BS^-1*DS*s_-1+BS^-1*ES*e_1+BS^-1*d


            %Using the solution for x_T then the recursion is given by
            %s_t=G(T-t+1)*s_t-1+h(T-t+1), 3<=t<=T
            %s_1=G(T)*s_0+(I-BS^-1*A*G(T-1))^-1*BS^-1*DS*e_1
            %where
            %G(T-t+1)=(I-BS^-1*A*G(T-t))^-1*BS^-1*C
            %h(T-t+1)=(I-BS^-1*A*G(T-t))^-1*(BS^-1*h(T-t)+d)


            for t=2:T
               G(:,1+(t-1)*nvar:t*nvar)=(I-AB*G(:,1+(t-2)*nvar:(t-1)*nvar))\CB;
               h(:,t)=(I-AB*G(:,1+(t-2)*nvar:(t-1)*nvar))\(AB*h(:,t-1)+db);
            end

            %If s_0=0, then s_1=G(T)*s_0+(I-BS^-1*A*G(T-1))^-1*BS^-1*D*e_1+h(T)
            %The initial condition is given by x_t-1

            G=[P,G];
            
            
%             disp(['Determinant = ' num2str(det(I-AB*G(:,1+(T-1)*nvar:T*nvar)))])
            if det(I-AB*G(:,1+(T-1)*nvar:T*nvar))<=0
                % This determinant has to be positive for a unique bounded
                % solution to exist.
                
                disp('Warning: Conditions for unique bounded solution violated.')
                disp('Adjust shock size or persistence downward.')
                disp('Exiting program...')
                return
            end

            s(:,1)=(I-AB*G(:,1+(T-1)*nvar:T*nvar))\DB*shock+h(:,T);     %G(:,1+T*nvar:(T+1)*nvar)*x(:,i)+



            %Using recursive coefficients to solve forward until the lower bound ceases
            %to bind

            for t=2:T
                s(:,t)=G(:,1+(T-t+1)*nvar:(T-t+2)*nvar)*s(:,t-1)+h(:,T-t+1);
            end

            %Impulse response after the zero bound ceases to bind

            s(:,T+1)=P*s(:,T);
            s(:,T+2)=P*s(:,T+1);


            checklb=1;                   %cannot move the liquidity trap length to T=-1
            checkub=(rb<s(rpos,T+1));   %checks that the zero lower bound does not bind tomorrw

        end
        
        T   =   T+(1-checkub);
    
    end
    
    for kk  = T+2:irflen2+1
        s(:,kk) =   P*s(:,kk-1);
    end   
    
    if zz   ==   1
        irfset.base =   s;
    elseif zz   ==  2
        eval(['irfset.' shocknames{mpos} ' =   s - irfset.base;']);
    end
    
end

close all

%% Graph Output: Tech shock
figure(1)
plot(0:1:irflen2,irf(ypos,:),'b-',0:1:irflen2,irfset.em(ypos,:),'r--','Linewidth',3);
ylabel('Output in Smets and Wouters (2007)')
xlabel('Quarters')
line([min(0) max(irflen)],[0,0],'LineWidth',1,'Color','k','Linestyle',':')
legend('Normal Times','8 Quarter ZLB')
% eval(['export_fig -transparent irf' shocknames{1} 'y.pdf'])

figure(2)
plot(0:1:irflen2,irf(cpos,:),'b-',0:1:irflen2,irfset.em(cpos,:),'r--','Linewidth',3);
ylabel('Consumption in Smets and Wouters (2007)')
xlabel('Quarters')
line([min(0) max(irflen)],[0,0],'LineWidth',1,'Color','k','Linestyle',':')
legend('Normal Times','8 Quarter ZLB')
% eval(['export_fig -transparent irf' shocknames{1} 'c.pdf'])

figure(3)
plot(0:1:irflen2,irf(pinfpos,:),'b-',0:1:irflen2,irfset.em(pinfpos,:),'r--','Linewidth',3);
ylabel('Inflation in Smets and Wouters (2007)')
xlabel('Quarters')
line([min(0) max(irflen)],[0,0],'LineWidth',1,'Color','k','Linestyle',':')
legend('Normal Times','8 Quarter ZLB')
% eval(['export_fig -transparent irf' shocknames{1} 'pinf.pdf'])

figure(4)
plot(0:1:irflen2,irf(rpos,:),'b-',0:1:irflen2,irfset.em(rpos,:),'r--','Linewidth',3);
ylabel('Inflation in Smets and Wouters (2007)')
xlabel('Quarters')
line([min(0) max(irflen)],[0,0],'LineWidth',1,'Color','k','Linestyle',':')
legend('Normal Times','8 Quarter ZLB')
% eval(['export_fig -transparent irf' shocknames{1} 'r.pdf'])

dy2(1)	=   max(exp(irfset.em(ypos,:))-1);
dy2(2)  =   exp(irfset.em(ypos,1))-1;
dc2(1)	=   max(exp(irfset.em(cpos,:))-1);
dc2(2)  =   exp(irfset.em(cpos,1))-1;
dr2     =   prod(exp(irfset.em(rpos,1:irflen2)-irfset.em(pinfpos,1:irflen2)))^(1/((irflen2)/4))-1;

dydr2   =   dy2./dr2;
dcdr2   =   dc2./dr2;

display(['ZLB: dy/dr is ' num2str(dydr2(1),2) ' at peak and ' num2str(dydr2(2),2) ' at inception.'])
display(['ZLB: dc/dr is ' num2str(dcdr2(1),2) ' at peak and ' num2str(dcdr2(2),2) ' at inception.'])

