var c_b c_d h_b h_d N_H N_C p lambda_b mc inflation 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 er;

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

set_param_value('phi',phi)
set_param_value('sigma',sigma)
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*p = (c_b)^(-1);
 
% 2) EULER EQUATION BORROWERS (FOC_B)

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_bh^(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*p + q*(h_b - (1-delta)*h_b(-1)) + B = w_c*n_bh + w_c*n_bc - r*B(-1) ;

% 8) HOUSING INVESTMENT BORROWERS

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

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

% 9) FOC c_d

lambda_d*p = (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*p + 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

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

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

% 17) PROFITS

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

% 18) WAGE IN CONSUMPTION SECTOR

%w_c = z;

% 19) WAGE IN HOUSING SECTOR

w_c = q*z*Lambda;

% 20) OPTIMAL PRICE SETTING

((1 - sigma) + sigma*mc) = theta*(inflation - 1)*inflation - theta*((Lambda(+1)/Lambda)*(q(+1)/q)*(Y_H(+1)/Y_H)*(inflation(+1) - 1)*inflation(+1));

% 21) SUPPLY H

 Y_H = z*N_H;

% 22) SUPPLY C

 Y_C = z*N_C;

%------------------------- AGGREGATE AND INTERESTING 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;

% 28) MARGINAL COST

 mc = w_h/q;

% 29) INFLATION IN HOUSING PRICES

 inflation = q/q(-1);

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

% 30) MARKET CLEARING DURABLE GOODS

Y_C = C;

% 31) MARKET CLEARING HOUSING

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

% 32) CLEARING CREDIT MARKET: LOANS = DEPOSITS

%D = B;

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

% 33) PRODUCTIVITY LAW OF MOTION

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

% 34) INTEREST RATE RULE

r = (H*inflation+C*p/p(-1))^(phi)*er*r(-1);

end;

%------------------------- SS INITIAL VALUES -------------------------%

initval;

p=1;
c_b = 1;
h_b = 5;
r =1/beta_d;
lambda_b = c_b^(-1); 
mu = (1-beta_b*r)*lambda_b ;
B = D;
mc=(sigma-1)/sigma;
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=0.5;
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;

steady(maxit=10000,solve_algo=4);
check;

shocks;
var ez; stderr sig_ez;
var er; stderr sig_ez;
end;

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