%----------------------------------------------------------
% Generate simulated data using a New Keynesian DSGE model
% with zero lower bound on the nominal interest rate
%----------------------------------------------------------
clear all;
randn('state', 777);

% Sample size and other global variables
nobs    = 3;		        % number of observed variables
simsmp  = 200;			% sample size in simulation
simdrop = 100;			% number of simulations to be dropped
simtt = simdrop + simsmp;	% total number of simulations
k = 3;				% number of structural shocks (excluding binding shocks)

nrep = 500;           % number of replication

%----------------------------------------
%   Assign parameter values
%----------------------------------------
sig    = 1;     % inverse of the intertemporal elasticity of substitution
chi    = 2;     % inverse of the elasticity of labor supply 
gam    = 0.5;   % habit persistence
xi     = 0.75;   % Calvo parameter
iota   = 0.5;   % inflation indexation
phi_pi = 1.5;   % policy response to inflation
phi_y  = 0.5; % policy response to output gap
phi_r  = 0.5;   % policy smoothing
piss   = 0.5;   % steady-state inflation
rss    = 0.5;  % steady-state real interest rate
rho_d  = 0.5;   % AR(1) coef on demand shock 
rho_p  = 0.5;   % AR(1) coef on cost shock 
rho_r  = 0.5;   % AR(1) coef on monetary policy shock 
sig_d  = 0.25;   % s.d. of demand shock
sig_p  = 0.25;   % s.d. of cost shock
sig_r  = 0.25;   % s.d. of policy shock

noeval = 0;

% auxiliary parameters
bet = 1/(1+rss/100);            % subjective discount factor
kap  = (1-xi)*(1-xi*bet)/xi;    % slope of the Phillips curve

%----------------------------------------
%  Matrices of canonical system (Sims' form)
%  GAM0*s(t) = GAM1*s(t-1) + PSI*e(t) + PPI*eta(t)
%----------------------------------------

% equation indices
eq_is    = 1;
eq_pc    = 2;
eq_rule  = 3;
eq_zd    = 4;
eq_zp    = 5;
eq_zr    = 6;
eq_Ey    = 7;
eq_Epi   = 8;
eq_m0    = 9;
eq_m1    = 10;
eq_m2    = 11;
eq_m3    = 12;
eq_m4    = 13;
eq_m5    = 14;
eq_m6    = 15;
eq_m7    = 16;
eq_m8    = 17;
eq_m9    = 18;
eq_m10   = 19;
eq_m11   = 20;
eq_m12   = 21;

% variable indices
v_y     = 1;  % output gap
v_pi    = 2;  % inflation rate
v_r     = 3;  % nominal interest rate
v_zd    = 4;  % demand shock
v_zp    = 5;  % cost shock
v_zr    = 6;  % monetary policy shock
v_Ey    = 7;  % E(t)y(t+1)
v_Epi   = 8;  % E(t)pi(t+1)
v_m0    = 9;  % binding shock with unanticipated component
v_m1    = 10; % 1-period-ahead anticipated binding shock
v_m2    = 11; % 2-period-ahead anticipated binding shock
v_m3    = 12; % 3-period-ahead anticipated binding shock
v_m4    = 13; % 4-period-ahead anticipated binding shock
v_m5    = 14; % 5-period-ahead anticipated binding shock
v_m6    = 15; % 6-period-ahead anticipated binding shock
v_m7    = 16; % 7-period-ahead anticipated binding shock
v_m8    = 17; % 8-period-ahead anticipated binding shock
v_m9    = 18; % 9-period-ahead anticipated binding shock
v_m10   = 19; % 10-period-ahead anticipated binding shock
v_m11   = 20; % 11-period-ahead anticipated binding shock
v_m12   = 21; % 12-period-ahead anticipated binding shock

% shock indices
e_d  = 1;
e_p  = 2;
e_r  = 3;
e_m0 = 4;
e_m1 = 5;
e_m2 = 6;
e_m3 = 7;
e_m4 = 8;
e_m5 = 9;
e_m6 = 10;
e_m7 = 11;
e_m8 = 12;
e_m9 = 13;
e_m10= 14;
e_m11= 15;
e_m12= 16;

% expectation error indices
n_y  = 1;
n_pi = 2;

% summary
neq  = 21;
neps = 16;
neta = 2;

% initialize matrices
GAM0 = zeros(neq,neq);
GAM1 = zeros(neq,neq);
C = zeros(neq,1);
PSI = zeros(neq,neps);
PPI = zeros(neq,neta);

%----------------------------------------
%      1. IS curve
%----------------------------------------
GAM0(eq_is,v_y)   =  1;
GAM0(eq_is,v_Ey)  = -1/(1+gam);
GAM0(eq_is,v_r)   =  (1-gam)/(1+gam)/sig;
GAM0(eq_is,v_Epi) = -(1-gam)/(1+gam)/sig;
GAM0(eq_is,v_zd)  = -1;
GAM1(eq_is,v_y)   =  gam/(1+gam);

%----------------------------------------
%      2. Phillips curve
%----------------------------------------
GAM0(eq_pc,v_pi)  = 1;
GAM0(eq_pc,v_Epi) = -bet/(1+bet*iota);
GAM0(eq_pc,v_y)   = -kap*(chi+sig/(1-gam))/(1+bet*iota);
GAM0(eq_pc,v_zp)  = -1;
GAM1(eq_pc,v_pi)  = iota/(1+bet*iota);
GAM1(eq_pc,v_y)   = -kap*(sig*gam/(1-gam))/(1+bet*iota);

%----------------------------------------
%      3. Monetary policy rule
%----------------------------------------
GAM0(eq_rule,v_r)  = 1;
GAM0(eq_rule,v_pi) = -(1-phi_r)*phi_pi;
GAM0(eq_rule,v_y)  = -(1-phi_r)*phi_y;
GAM0(eq_rule,v_zr) = -1;
GAM0(eq_rule,v_m0) = -1;
GAM1(eq_rule,v_r)  = phi_r;
%GAM1(eq_rule,v_y)  = -(1-phi_r)*phi_y;  %output growth rule

%----------------------------------------
%      4.- Shock Processes
%----------------------------------------
% (unanticipated) demand shock
GAM0(eq_zd,v_zd)  = 1;
GAM1(eq_zd,v_zd)  = rho_d;
PSI(eq_zd,e_d)    = 1;

% (unanticipated) cost shock
GAM0(eq_zp,v_zp)  = 1;
GAM1(eq_zp,v_zp)  = rho_p;
PSI(eq_zp,e_p)    = 1;

% (unanticipated) monetary policy shock
GAM0(eq_zr,v_zr)  = 1;
GAM1(eq_zr,v_zr)  = rho_r;
PSI(eq_zr,e_r)    = 1;

% binding shock with unanticipated component
GAM0(eq_m0,v_m0)  = 1;
PSI(eq_m0,e_m0)   = 1;
GAM1(eq_m0,v_m1) = 1;

% 1-period ahead anticipated binding shock
GAM0(eq_m1,v_m1) = 1;
PSI(eq_m1,e_m1)   = 1;
GAM1(eq_m1,v_m2) = 1;

% 2-period ahead anticipated binding shock
GAM0(eq_m2,v_m2) = 1;
PSI(eq_m2,e_m2)   = 1;
GAM1(eq_m2,v_m3) = 1;

% 3-period ahead anticipated binding shock
GAM0(eq_m3,v_m3) = 1;
PSI(eq_m3,e_m3)   = 1;
GAM1(eq_m3,v_m4) = 1;

% 4-period ahead anticipated binding shock
GAM0(eq_m4,v_m4) = 1;
PSI(eq_m4,e_m4)   = 1;
GAM1(eq_m4,v_m5) = 1;

% 5-period ahead anticipated binding shock
GAM0(eq_m5,v_m5) = 1;
PSI(eq_m5,e_m5)   = 1;
GAM1(eq_m5,v_m6) = 1;

% 6-period ahead anticipated binding shock
GAM0(eq_m6,v_m6) = 1;
PSI(eq_m6,e_m6)   = 1;
GAM1(eq_m6,v_m7) = 1;

% 7-period ahead anticipated binding shock
GAM0(eq_m7,v_m7) = 1;
PSI(eq_m7,e_m7)   = 1;
GAM1(eq_m7,v_m8) = 1;

% 8-period ahead anticipated binding shock
GAM0(eq_m8,v_m8) = 1;
PSI(eq_m8,e_m8)   = 1;
GAM1(eq_m8,v_m9) = 1;

% 9-period ahead anticipated binding shock
GAM0(eq_m9,v_m9) = 1;
PSI(eq_m9,e_m9)   = 1;
GAM1(eq_m9,v_m10) = 1;

% 10-period ahead anticipated binding shock
GAM0(eq_m10,v_m10) = 1;
PSI(eq_m10,e_m10)   = 1;
GAM1(eq_m10,v_m11) = 1;

% 11-period ahead anticipated binding shock
GAM0(eq_m11,v_m11) = 1;
PSI(eq_m11,e_m11)   = 1;
GAM1(eq_m11,v_m12) = 1;

% 10-period ahead anticipated binding shock
GAM0(eq_m12,v_m12) = 1;
PSI(eq_m12,e_m12)   = 1;

%----------------------------------------
%      6.- 7. Expectation error
%----------------------------------------
GAM0(eq_Ey,v_y)  = 1;
GAM1(eq_Ey,v_Ey) = 1;
PPI(eq_Ey,n_y)  = 1;

GAM0(eq_Epi,v_pi)  = 1;
GAM1(eq_Epi,v_Epi) = 1;
PPI(eq_Epi,n_pi)  = 1;

%----------------------------------------
%  Solve the model: QZ(generalized Schur) decomposition by GENSYS
%  -> s(t) = G1*s(t-1) + impact*e(t)
%----------------------------------------
[G1,const,impact,fmat,fwt,ywt,gev,eu,loose]=gensys(GAM0,GAM1,C,PSI,PPI);

if min(eu)<1
    noeval = 1;
end

% augmentation: v_ylag
%v_ylag = neq + 1;
%
%GG = [G1 zeros(neq,1); 1 zeros(1,neq)];
GG = G1;
%RR = [impact; zeros(1,neps)];
RR = impact;

%----------------------------------------
%  Measurement equations
%   Output growth rate = y;
%   Inflation rate     = piss + pi;
%   FF rate            = piss + rss + r;
%
%   y(t) = DD + ZZ*s(t)
%----------------------------------------

% observation equation indices
eq_yobs  = 1;
eq_piobs = 2;
eq_robs  = 3;

% observation variable indices
v_yobs  = 1;
v_piobs = 2;
v_robs  = 3;

DD = zeros(nobs,1);
DD(eq_piobs,1) = piss;
DD(eq_robs,1) = piss + rss;

ZZ = zeros(nobs,neq);
ZZ(eq_yobs,v_y) = 1;
ZZ(eq_piobs,v_pi) = 1;
ZZ(eq_robs,v_r) = 1;

%----------------------------------------
% Generate simulated series
%----------------------------------------
simyy = zeros(simsmp,nobs,nrep);
probzlbs = zeros(1,nrep);
distzlbs = [];

for rep = 1:nrep
    
sims = zeros(simtt,neq);
simy = zeros(simtt,nobs);

n = 15;    % number of period for IRFs
irfszlb = zeros(n,neq);

% shocks to be used in simulation
shocks = randn(k,simtt);

s = 1;
while s <= simtt
    
    st = irfszlb(1,:)';		% initially zero. for s>1, simulated values are used.
    stdevr = [sig_d sig_p sig_r];
    et = zeros(neps,1);
    et(1:size(stdevr,2),1) = stdevr'.*shocks(:,s);
    sts = zeros(neq,n);
    yts = zeros(nobs,n);
    
    % Compute IRFs without zero lower bound
    for t=1:n
        if t==1;
            st = GG*st+RR*et;
        else
            st = GG*st;
        end
        sts(:,t) = st;
        yts(:,t) = DD + ZZ*st;
    end
    
    irfs = sts';
    irfy = yts';
    
    % Find the periods when r drops below zero
    ind = find(yts(v_robs,:) < 0);
    
    % Compute IRFs with zero lower bound
    if size(ind,2)>0
        
        if max(ind)>13
            noeval = 1;
            break
        end
        
        % Add unanticipated and anticipated MP disturbances
        % until the negative r_obs gets close to zero
        itr = 1;
        while any(abs(yts(v_robs,ind)) > 0.001)
            
            % Shocks occur at t.
            % If r(t) is negative, add unanticipated binding shock e_m0(t).
            % If r(t+1) is negative, add 1-period-ahead anticipated binding shock e_m1(t).
            % If r(t+2) is negative, add 2-period-ahead anticipated binding shock e_m2(t).
            % And so on.
            
            for i=ind
                et(e_m0+i-1,1)=et(e_m0+i-1,1)-yts(v_robs,i);
            end
            
            for t=1:n
                if t==1;
                    st = irfszlb(1,:)';
                    st = GG*st+RR*et;
                else
                    st = GG*st;
                end
                sts(:,t)=st;
                yts(:,t)=DD+ZZ*st;
            end
            
            if abs(yts(eq_yobs,1))>10 || abs(yts(eq_piobs,1))>10
                noeval = 1;
                break
            end
            
            % The addition of MP disturbances sometimes generates more periods when r_obs becomes negative.
            % Such a possibility is considered below.
            ind2 = find(yts(v_robs,:) < 0);
            if size(ind2,2)>0
                if (min(ind2) < min(ind)) || (max(ind2) > max(ind))
                    ind2dif = setdiff(ind2,ind);
                    ind = [ind ind2dif];
                end
            end
            
            if max(ind)>13
                noeval = 1;
                break
            end
            
            itr = itr +1;
            if itr == 200
                noeval = 1;
                break
            end
        end            
    end
    
    irfszlb = sts';
    irfyzlb = yts';
    
    sims(s,:) = irfszlb(1,:);
    simy(s,:) = irfyzlb(1,:);
    
    if noeval==1
        if s==1
            irfszlb(1,:) = zeros(1,neq);
        else
            irfszlb(1,:) = sims(s-1,:);
        end
        shocks(:,s)=randn(k,1);   % Resample structural shocks.
        s = s-1;
        noeval=0;
    end
    
    s = s + 1;

end

if noeval==0
    simy(1:simdrop,:) = [];     % eliminate initial simulations
    %disp('"simy" contains simulated data.')
else
    disp('Failure in Simulation!')
end

simyy(:,:,rep) = simy;

tzlb = find(abs(simy(:,v_robs))<0.001);
nzlb = size(tzlb,1);
probzlb = nzlb/simsmp;
probzlbs(1,rep) = probzlb;

durzlb = 1;
distzlb = [];

if nzlb==0
    distzlb=0;
elseif nzlb ==1
    distzlb=1;
else
    for u = 2:nzlb
        if tzlb(u)-tzlb(u-1)==1
            durzlb = durzlb +1;
        else
            distzlb = [distzlb;durzlb];
            durzlb = 1;
        end
    end
    distzlb = [distzlb;durzlb];
end

distzlbs = [distzlbs;distzlb];

end

distzlbs
mean(probzlbs)
mean(distzlbs)
hist(distzlbs)

save('simyy','simyy');