close all

var LAM R PI RK C W L Y A K MC B G U I RR TAX;
varexo eps_A eps_U eps_G;

parameters beta delta sigmaL omega alpha phiP Phi tau jotaR jotaPI jotaY rhoA rhoU rhoG sigmaA sigmaU sigmaG nu 
TAXss Gss Yss Rss sc AUX Kss Lss Css Wss Iss RKss LAMss Ass Uss Bss t;

beta = 0.90;
delta = 0.025; 
sigmaL = 5;
omega = 1;
alpha = 2/3;
phiP = 1.25;
Phi = 1.5;
tau = 5;
jotaR = 0;
jotaPI = 1.5;
jotaY = 0.125;
rhoA = 0.90;
rhoU = 0.90;
rhoG = 0.90;
sigmaA = 0.01;
sigmaU = 0.01;
sigmaG = 0.01;
nu = 0.8;
sc = 0.2;
t = 0.1;
AUX = ((1-alpha)*(1-1/Phi)/(1/beta - 1 + delta))^((1-alpha)/alpha);
Wss = AUX *alpha*(1-1/Phi);
Lss = ((1-sc-delta*AUX^(alpha/(1-alpha)))*sigmaL/(sigmaL-1)*omega*1/alpha*(1-1/Phi)^(-1))^(1/(-nu-1));
Kss = Lss*AUX^(1/(1-alpha));
Yss = Lss^(alpha)*Kss^(1-alpha);
Gss = sc*Yss;
Css = Yss - delta*Kss - Gss;
Rss = 1/beta; 
TAXss = t*Yss;
Iss = delta*Kss;
RKss = 1/beta - 1 + delta;
LAMss = 1/Css;
Ass = 1;
Uss = 1;
Bss = (TAXss - Gss)/(Rss-tau);

model;
LAM = beta*LAM(+1)*R/PI(+1);
LAM = beta*LAM(+1)*(RK(+1) + (1-delta));
LAM = 1/C;
W = sigmaL/(sigmaL-1)*omega*L^(nu)*1/LAM;
Y = A*L^(alpha)*K(-1)^(1-alpha);
RK = (1-alpha)*MC*A*L^(alpha)*K(-1)^(-alpha);
W = alpha*MC*A*L^(alpha-1)*K(-1)^(1-alpha);
1 - phiP*(PI-1)*PI + beta*phiP*LAM(+1)/LAM*(PI(+1)-1)*Y(+1)/Y*PI(+1) - (1-MC)*Phi = 0;
B = R(-1)*B(-1) + G - TAX;
TAX = TAXss + tau*B;  
R/Rss = (R(-1)/Rss)^(jotaR)*((1+PI)^(jotaPI)*(Y/Yss)^(jotaY))^(1-jotaR)*U;
K = (1-delta)*K(-1) + I;  
Y = C + I + G + phiP/2*(PI-1)^2*Y;
RR = R/PI(+1);
A = ((Ass)^(1-rhoA) )*( A(-1)^rhoA )*exp(eps_A);
U = ((Uss)^(1-rhoU) )*( U(-1)^rhoU )*exp(eps_U);
G = ((Gss)^(1-rhoG) )*( G(-1)^rhoG )*exp(eps_G);
end;


initval;
B = Bss;
PI = 1;
A = Ass;
U = Uss;
G = Gss^(1-rhoG);
MC = (Phi-1)/Phi;
TAX = TAXss;
R = Rss;
W = Wss; 
L = Lss;
K = Kss;
Y = Yss;
C = Css;
LAM = LAMss;
RK = RKss;
I = Iss;
RR = Rss;
end;

check;
steady;

shocks;
var eps_A = sigmaA;
var eps_U = sigmaU;
var eps_G = sigmaG;
end;

stoch_simul(periods=1000,irf=20,order=1);