var r P P_F P_I w_F w_I w K_F K_I K A_F A_I g H_F H_I C C_F C_I Y Y_F Y_I M;

varexo e_A_F e_A_I e_g;
parameters upsilon sigma beta theta delta psi iota 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.36;
delta = 0.25;
rho = 0.75;
psi = 11;
iota = 11;
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 = (((iota-1)*(1-alpha)^(1-alpha)*alpha^alpha)/(iota*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/iota)+(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 - ((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; //FOC of intermediate good firm
0 = A_F+(1-theta)*H_F+theta*K_F-Y_F; //aggregate production function


//Informal
0 = beta*P_I(+1)+((1-alpha)*(1-rho)*(1-beta*rho))/rho * w_I + (alpha*(1-rho)*(1-beta*rho))/rho * r - ((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; //FOC of intermediate good firm
0 = A_I+(1-alpha)*H_I+alpha*K_I-Y_I; //aggregate production function


0 = w+beta*steady_r*r(+1)-w(+1); //FOC households
0 = g+M(-1)-P-C;
0 = steady_P_F*gamma^(mu)*P_F+steady_P_I*(1-gamma)^(mu)*P_I -steady_P*P; //aggregated price
0 = steady_C_F*gamma^(mu/(mu-1))*C_F + steady_C_I*(1-gamma)^(mu/(mu-1))*C_I-steady_C*C; //aggregated consumption
0 = (1-delta)*steady_K*K(-1)-(steady_M/steady_P)*M + (steady_M/steady_P)*P +steady_Y*Y -steady_K*K; //budget restriction
0 = w+P-pi*g-M; // 2nd FOC and cash-in-advance combined
0 = upsilon^((1+theta_L)/theta_L)*steady_w_F*w_F+(1-upsilon)^((1+theta_L)/theta_L)*steady_w_I*w_I-steady_w*w; //agregated wage
0 = M-M(-1)-g; //aggregated money stock
0 = A_F*omega-A_I; //relationship technologies
0 = Y_F+Y_I-Y; // aggregated income (household)
0 = K_F+K_I-K; //aggregated capital
0 = pi*g-g(+1); // law of motion
0 = theta_L*H_F+w_F-theta_L*H_I-w_I; //demand formal and informal labor

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

initval;
   r = 0;
 P_F = 0;
 P_I = 0;
   P = 0;
 w_F = 0;
 w_I = 0;
   w = 0;
 K_F = 0;
 K_I = 0;
   K = 0;
 A_F = 0;
 A_I = 0;
   g = 0;
 H_F = 0;
 H_I = 0;
   C = 0;
 C_F = 0;
 C_I = 0;
 Y_F = 0; 
 Y_I = 0;
   Y = 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) ;