var c_b c_d h_b h_d N_H N_C lambda_b lambda_d mu C H IH i_hd i_hb B D r z q Lambda pi w_h w_c n_dh n_bh n_dc n_bc Y_H Y_C;
varexo ez;

parameters theta,beta_b,beta_d,rho,m,delta,sig_ez,mi,epsilon_b,epsilon_d,j;

set_param_value('theta',theta)
set_param_value('delta',delta)
set_param_value('beta_b',beta_b)
set_param_value('beta_d',beta_d)
set_param_value('rho',rho)
set_param_value('sig_ez',sig_ez)
set_param_value('m',m)
set_param_value('epsilon_d',epsilon_d)
set_param_value('epsilon_b',epsilon_b)
set_param_value('mi',mi)
set_param_value('j',j)

model;

%---------------------------- BORROWERS ----------------------------------%

% 1) FOC c_b
 
lambda_b = (c_b)^(-1);
 
% 2) EULER EQUATION BORROWERS (FOC_B): mu is bigger then 0 so the collateral constraint here is binding

mu = lambda_b - beta_b*r(+1)*lambda_b(+1);

% 3) BINDING COLLATERAL CONSTRAINT

B = h_b*m*q;

% 4) FOC h_b

q*lambda_b = j*h_b^(-1) + m*q*mu + beta_b*lambda_b(+1)*q(+1)*(1-delta);

% 5) FOC n_bc

w_c*lambda_b = (n_bc^(epsilon_b))*(n_bc^(1+epsilon_b) + n_bc^(1+epsilon_b))^(((1+mi)/(1+epsilon_b))-1);

% 6) FOC n_bh

w_c*lambda_b = (n_bh^(epsilon_b))*(n_bc^(1+epsilon_b) + n_bh^(1+epsilon_b))^(((1+mi)/(1+epsilon_b))-1);

% 7) BUDGET CONSTRAINT

c_b + q*(h_b - (1-delta)*h_b(-1)) + B = w_c*n_bh + w_c*n_bc - r*B(-1) ;

% 8) HOUSING INVESTMENT BORROWERS

(h_b - (1-delta)*h_b(-1)) = i_hb;

%---------------------------- DEPOSITORS ----------------------------------%

% 9) FOC c_d

lambda_d = (c_d)^(-1);
 
% 10) EULER EQUATION (FOC_D)

lambda_d = beta_d*r(+1)*lambda_d(+1);

% 11) FOC h_b

q*lambda_d = j*h_d^(-1) + beta_d*lambda_d(+1)*q(+1)*(1-delta);

% 12) FOC n_dc

w_c*lambda_d = (n_dc^(epsilon_d))*(n_dc^(1+epsilon_d) + n_dh^(1+epsilon_d))^(((1+mi)/(1+epsilon_d))-1);

% 13) FOC n_dh 

w_c*lambda_d = (n_dh^(epsilon_d))*(n_dc^(1+epsilon_d) + n_dh^(1+epsilon_d))^(((1+mi)/(1+epsilon_d))-1);

% 14) BUDGET CONSTRAINT

c_d + q*(h_d - (1-delta)*h_d(-1)) - D = w_c*n_dh + w_c*n_dc + r*B(-1) + pi;

% 15) STOCHASTIC DISCOUNT FACTOR

Lambda = beta_d*lambda_d(+1)*lambda_d^(-1);

% 16) HOUSING INVESTMENT DEPOSITORS

(h_d - (1-delta)*h_d(-1)) = i_hd;

%---------------------------- FIRMS --------------------------------------%

% 17) PROFITS

pi = q*Y_H - w_h*N_H - q*Y_H*theta^(-1)*(q*q(-1)^(-1)-1)^(2);

% 18) WAGE IN CONSUMPTION SECTOR

w_c = z;

% 19) WAGE IN HOUSING SECTOR

w_h = q*z*Lambda;

% 20) OPTIMAL PRICE SETTING

Lambda*(z*N_H - Y_H*theta*2^(-1)*(q*q(-1)^(-1)-1)^(2) - Y_H*q*q(-1)^(-1)*2*theta*(q*q(-1)^(-1)-1)) = Lambda(+1)*(((q(+1)/q)^(2))*Y_H(+1)*2*theta*(q(+1)*q^(-1)-1));

% 21) SUPPLY H

Y_H = z*N_H;

% 22) SUPPLY C

Y_C = z*N_C;

%------------------------- AGGREGATE VARIABLES ---------------------------%

% 23) AGGREGATE LABOR IN HOUSING SECTOR

N_H = n_bh + n_dh;

% 24) AGGREGATE LABOR IN CONSUMPTION SECTOR

N_C = n_bc + n_dc;

% 25) AGGRGATE CONSUMPTION

C = c_d + c_b;

% 26) AGGREGATE HOUSING

H = h_d + h_b;

% 27) AGGREGATE HOUSING INVESTMENT

IH = i_hb + i_hd;

%------------------------------- MARKET CLEARING -------------------------%

% 28) MARKET CLEARING DURABLE GOODS

Y_C = C;

% 29) MARKET CLEARING HOUSING

Y_H = IH + q*Y_H*theta*2^(-1)*(q*q(-1)^(-1)-1)^(2);

% 30) CLEARING CREDIT MARKET: LOANS = DEPOSITS

D = B;

%--------------------- EXOGENOUS PROCESS ---------------------------------%

% 31) PRODUCTIVITY LAW OF MOTION

z = (1-rho) + rho*z(-1) + sig_ez*ez;

end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% SS initial values %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

initval;

c_b = 1;
h_b = 5;
r =1/beta_d;
lambda_b = c_b^(-1); 
mu = (1-beta_b*r)*lambda_b ;
B = D;
w_c = 0.8;
w_h =w_c;
n_bc =1;
n_bh =1;
c_d=1;
lambda_d =c_d^(-1);
h_d=6;
H = h_d + h_b;
C =c_d+c_b;
D =0.5;
n_dc=0.5;
n_dh=0.5;
z=1;
q=2;
N_H=n_bh+n_dh;
N_C=n_bc+n_dc;
Y_C=N_C;
Y_H=N_H;
i_hb=delta*h_b;
i_hd =delta*h_d;
IH = delta*H;
pi=1;
Lambda=beta_d;

end;

resid(1);
steady(solve_algo=4);
check;

shocks;
var ez; stderr sig_ez;
end;

stoch_simul(periods=2100, order=1,nocorr, drop=400, nomoments,IRF=0,replic=1000);