var r_F r_I P_F P_I w_F w_I K_F K_I A_F A_I g H_F H_I C_F C_I Y_F Y_I M;

varexo e_A_F e_A_I e_g;
parameters upsilon sigma beta theta delta psi epsilon alpha rho pi B A h0 gamma mu phi omega theta_L steady_Y steady_r steady_K steady_w steady_w_F steady_w_I steady_g  steady_C steady_C_F steady_C_I steady_K_F steady_K_I steady_Y_F steady_Y_I steady_P steady_P_F steady_P_I steady_H_F steady_H_I steady_M;

sigma = 0.007;
beta = 0.99;
theta = 0.26;
delta = 0.011;
rho = 0.75;
psi = 11;
epsilon = 9;
alpha = 0.2;
gamma = 0.426;
phi = 0.95;
omega = 0.1;
A = 1.72;
upsilon = 0.29;
mu = 0.7;
theta_L = 2;
h0 = 0.583;
B = (A*log(1-h0))/h0;


steady_r = 1/beta-(1-delta);

steady_w_F = (((psi-1)*(1-theta)^(1-theta)*theta^theta)/(psi*steady_r^theta))^(1/(1-theta));
steady_w_I = (((epsilon-1)*(1-alpha)^(1-alpha)*alpha^alpha)/(epsilon*steady_r^alpha))^(1/(1-alpha));
steady_w = (upsilon*steady_w_F^(theta_L/(1+theta_L))+(1-upsilon)*steady_w_I^(theta_L/(1+theta_L)))^((1+theta_L)/theta_L);

steady_C_F = -1*(beta*steady_w_F)/B;
steady_C_I = -1*(beta*steady_w_I)/B;

steady_C = (gamma*steady_C_F^((mu-1)/mu)+(1-gamma)*steady_C_I^((mu-1)/mu))^((mu)/mu-1);
steady_M = 1;
steady_P = steady_M/steady_C;

steady_P_F = (steady_w_F^(1/(theta_L))*steady_P)/((gamma*steady_w)^(1/(theta_L)));
steady_P_I = (steady_w_I^(1/(theta_L))*steady_P)/(((1-gamma)*steady_w)^(1/(theta_L)));

steady_Y_F = (-1*beta*steady_w_F)/(B*(steady_w_F*((steady_r*(1-theta)/(steady_w_F*theta))^theta)+ (1/psi)+(steady_r-delta)*((steady_r*(1-theta))/(steady_w_F*theta))^(theta-1)));
steady_Y_I = (-1*beta*steady_w_I)/(B*(steady_w_I*((steady_r*(1-alpha)/(steady_w_I*alpha))^alpha)+ (1/epsilon)+(steady_r-delta)*((steady_r*(1-alpha))/(steady_w_I*alpha))^(alpha-1)));
steady_Y = steady_Y_F + steady_Y_I;

steady_K_F = ((steady_r*(1-theta)/(steady_w_F*theta))^(theta-1))*steady_Y_F;
steady_K_I = ((steady_r*(1-alpha)/(steady_w_I*alpha))^(alpha-1))*steady_Y_I;
steady_K = steady_K_F + steady_K_I;

steady_H_F = ((steady_r*(1-theta)/(steady_w_F*theta))^theta)*steady_Y_F;
steady_H_I = ((steady_r*(1-alpha)/(steady_w_I*alpha))^alpha)*steady_Y_I;

steady_g = 1;
pi = steady_g;

model(linear);
//Formal
0 = beta*P_F(+1)+((1-theta)*(1-rho)*(1-beta*rho))/rho * w_F + (theta*(1-rho)*(1-beta*rho))/rho * r_F - ((1-rho)*(1-beta*rho))/rho *A_F - (1+beta)*P_F+P_F(-1); //Pricing rule
0 = H_F+w_F-K_F-r_F; //FOC of intermediate good firm
0 = A_F+(1-theta)*H_F+theta*K_F-Y_F; //aggregate production function
0 = (1-delta)*steady_K_F*K_F(-1)-(steady_M/steady_P_F)*M + (steady_M/steady_P_F)*P_F +steady_Y_F*Y_F -steady_K_F*K_F; //budget restriction
0 = w_F+beta*steady_r*r_F(+1)-w_F(+1); // 1st FOC households 5
0 = w_F+P_F-C_F(+1)-P_F(+1); // 2nd FOC households oder w+P-M-pi*g
//0 = w_F+P_F-M-pi*g; // instead before
0 = g+M(-1)-P_F-C_F; // Cash-in-advance

//Informal
0 = beta*P_I(+1)+((1-alpha)*(1-rho)*(1-beta*rho))/rho * w_I + (alpha*(1-rho)*(1-beta*rho))/rho * r_I - ((1-rho)*(1-beta*rho))/rho *A_I - (1+beta)*P_I+P_I(-1); //Pricing rule
0 = H_I+w_I-K_I-r_I; //FOC of intermediate good firm instead of 0 = w+P-pi*g-M (2nd FOC and cash-in-advance combined) 
0 = A_I+(1-alpha)*H_I+alpha*K_I-Y_I; //aggregate production function
0 = (1-delta)*steady_K_I*K_I(-1)-(steady_M/steady_P_I)*M + (steady_M/steady_P_I)*P_I +steady_Y_I*Y_I -steady_K_I*K_I; //budget restriction 12
0 = w_I+beta*steady_r*r_I(+1)-w_I(+1); // 1st FOC households 12
0 = w_I+P_I-C_I(+1)-P_I(+1); // 2nd FOC households oder w+P-M-pi*g
//0 = w_I+P_I-M-pi*g; // instead before
0 = g+M(-1)-P_I-C_I; // Cash-in-advance

0 = M-M(-1)-g; //aggregated money stock 15

//Shocks
A_F = phi*A_F(-1)+e_A_F; //20
A_I = phi*A_I(-1)+e_A_I; //21
g = phi*g(-1)+e_g; //21

end;

initval;
   r_F = 0;
   r_I = 0;
 P_F = 0;
 P_I = 0;
 w_F = 0;
 w_I = 0;
 K_F = 0;
 K_I = 0;
 A_F = 0;
 A_I = 0;
   g = 0;
 H_F = 0;
 H_I = 0;
 C_F = 0;
 C_I = 0;
 Y_F = 0; 
 Y_I = 0;
   M = 0;
end;

steady;

shocks;

var e_A_F = sigma^2;
var e_A_I = sigma^2;
var e_g = sigma^2;

end;

stoch_simul(irf=100) ;