var y n c w kk ku r i pi mc pk u invv psii rkk zr a rn ygap;

varexo ez er;

parameters pi_s LAMBDA KAPPA PSI_YA RHO RHOPIT RHOPI BETTA DELTA XI ALFA YBAR IBAR CBAR CY IY KY THETA EPS FI_Y VARPHI SIG FI_PI RHO_A RHO_V SIGMAa Rk MCBAR NBAR KBAR;

% Calibration of parameters
SIG = 1; 
EPS = 6; 
BETTA = 0.99;
VARPHI = 1;
THETA = 2/3;
RHO_V = 0.50;
RHO_A = 0.50; 
FI_Y = 0.5/4;
ALFA= 0.33;
DELTA= 0.025;
XI = 2.48; 
SIGMAa = 0.01;              
FI_PI = 1.5;
RHOPIT = 0.50;
RHOPI = 0.50; 
RHO     =  1/BETTA-1;     // steady state interest rate
pi_s = 1;                //Zero steady state inflation

Rk = (1 - BETTA * (1 - DELTA)) / BETTA;
MCBAR = (EPS - 1)/EPS;
NBAR = 1/3;
KBAR = NBAR * (Rk / ALFA * (1/MCBAR) )^( 1/(ALFA-1));
YBAR = KBAR^(ALFA) * NBAR^(1-ALFA);
IBAR = DELTA * KBAR;
CBAR = YBAR - IBAR;
CY = CBAR/YBAR;
IY = IBAR/YBAR;
KY = KBAR/YBAR;
PSI_YA  =  (1+VARPHI)/(VARPHI+ALFA+SIG*(1-ALFA)); 
LAMBDA = ((1 - THETA*pi_s^(EPS -1))*(1 - BETTA*THETA*(pi_s^(EPS -1)))/(THETA*(pi_s^(EPS -1))))*(1-ALFA)/(1-ALFA+ALFA*EPS);
KAPPA   =  LAMBDA*(SIG + (VARPHI+ALFA)/(1-ALFA));


model(linear);

% BLOCK 1 HOUSEHOLDS

(-SIG) * c = psii; % Household Consumption FOC
VARPHI * n = (- SIG) * c + w; %Household Labor supply
kk = (1 - DELTA ) * kk(-1) + DELTA * invv; % Capital Accumulation
psii + pk = psii(+1) + BETTA * (1 - DELTA) * pk(+1) + (1 - BETTA * (1 - DELTA)) * (rkk(+1)); %Household capital FOC
rkk = SIGMAa * u; % Household FOC with respect to capital utilization
-i + psii(+1) - psii = pi(+1) - RHO; % The Euler equation
pk = XI * (invv - invv(-1) - BETTA * (invv(+1) - invv)); % Household FOC with respect to Investment

%BLOCK 2 FIRMS

w = mc + y - n; % wage equation
rkk = mc + y - ku; % capital rent
ku =  u + kk(-1); % capital utilization
pi = BETTA * pi(+1) +  KAPPA*ygap; %The New-Keynsian Phillips cruve

% BLOCK 3 MONETARY AUTHORITY
i      = rn + FI_PI*pi + FI_Y*(ygap) + zr; %Taylor Rule

%BLOCK 4 MARKET CLEARING
y = CY * c + IY * invv + Rk * KY * u; %goods market
y = a + ALFA * ku + (1 - ALFA) * n; %Production function
rn     = RHO + SIG*PSI_YA*(a(+1) - a);
ygap   = ygap(+1) - 1/SIG*(i-pi(+1)-rn);
r      = i - pi(+1);                            //Fischer equation

% BLOCK 5 STOCHASTIC SHOCKS

zr    = RHO_V*zr(-1) + er;
a     = RHO_A*a(-1) + ez;

end;

initval;

er=0.000;
ez=0.000;
i     = BETTA;
rn    = BETTA;
r   = RHO;
ygap  = 0;
y     = 0;
n     = 0;
pi    = 0.00;
end;

steady;
check;

shocks;
var ez; stderr 0.25; 
var er; stderr 0.25;        
end;

stoch_simul(irf=12,  order=1);