var
y
c_p
c_i
l
w
lambda_p
lambda_i
r
R_d
%R_star
npa
%s
pi
d
tau_w
tau_h
mu_dti
q
h
mu_ltv
D
L
phi_dti
phi_ltv
R_m
k_b
tr
q_h
q_l
mu_b
%Gamm
epz
mc
%sdf
B
%epsilon_h
tau_c
%tau_b
mu_n
lambda
delta_w
%b_f
epsm
%epsq
;

varexo epzss eta_ltv eta_dti epsilon_w epsilon_h epsilon_c epsmss eps_w eps_lambda eps_q;

%predetermined_variables  s;

parameters beta_p beta_i Gamm_bar kappa theta delta_h rho_ltv rho_dti phi_ltv_bar theta_tilde
phi_dti_bar delta_k f nu_b kk tauw_bar tauh_bar phi_wb phi_hb neta zeta epsilon q_bar h_bar y_bar
rhos rhom rhoz psi a_p a_i pi_bar r_bar alpha B_bar rho_w rho_h d_bar epz_bar rhoq rhol gamma phir epsm_bar D_bar ohm R_star_bar s_bar epsilon_h_bar rhoh
k_tilde k_b_bar A lambda_bar delta_w_bar omega L_bar s_tilde npa_bar tauc_bar taub_bar rho_c phi_cb rho_b phi_bb q_h_bar q_l_bar U_bar R_l_star_bar bf_bar rhod rholtv rhoqh epsq_bar;

beta_p = 0.984131; %0.986603;
beta_i = 0.975;
kappa = 0.033;
theta = 0.1;
delta_h = 0.01;
%theta_tilde = exp(2.6); %6.75;
rho_ltv = 0.95;
rho_dti = 0.95;
phi_ltv_bar = 0.78;
phi_dti_bar = 0.034754955; %2.7071;
delta_k = 0.1049;
f       = 0.016822909; 
nu_b    = 0.174483; %0.09;
kk      = 2;
tauw_bar = 0.19;
tauh_bar = 0.000624;
rho_w = 0.95;
rho_h = 0.95;
phi_wb = 0.0001;
phi_hb = 0.003;
neta = 1.2;
zeta = 0.5;
epsilon = 4/3;
q_bar = 1;
h_bar = 1;
y_bar = 1;
rhos = 0.99;
rhoz = 0.90;
rhom = 0.55;
psi  = 0.9;
a_p  = 0.1; %0.8;
a_i  = 0.3813; %0.6846; %0.7134; %0.713476;
pi_bar = 1;
r_bar = 0.016125; 
%R_bar = 1.016125;
alpha = 0.4;
gamma = 0.4;
B_bar = 0;
rho_w = 0.9;
rho_h = 0.9;
d_bar = 0.146009; %0.040647107;
epz_bar = 1;
rhoq = 0.90;
rhol = 0.90;
phir = 0.1;
epsm_bar = 1;
D_bar = 5.35969697;
R_star_bar= 1.2701561073;
R_l_star_bar = 1.276523405;
s_bar = 1;
ohm  = 1;
epsilon_h_bar = 1;
rhoh = 0.95;
A    = (1/0.33);
q_h_bar = 1.5;
q_l_bar = 0.5;
omega  = 0.018928005;
f      = 0.016822909;
delta_w_bar = 0.028784;
s_tilde = 1.863899317;
k_tilde = 0.005;
lambda_bar = 0.5;
U_bar = 0.030861;
L_bar = 0.176870;
npa_bar = 0.099961;
tauc_bar = 0.25;
taub_bar = 0.20;
rho_c = 0.9;
rho_b = 0.9;
phi_cb = 0.3;
phi_bb = 0.3;
k_b_bar = 0.030861;
Gamm_bar = 1;
theta_tilde = 9;
bf_bar = 0;
rhod = 0.9;
rholtv = 5;
rhoqh = 0.9;
epsq_bar = 1;


model;
//Savers (Equation 1 - 5)

lambda_p - (1-a_p)*(c_p - a_p*c_p(-1))^(-1) = 0;
lambda_p = beta_p*lambda_p(+1)*(1+r)*pi(+1)^(-1) ;
%lambda_p = beta_p*lambda_p(+1)*R_star*Gamm*(1-tau_b(+1))*pi(+1)^(-1);
lambda_p = beta_p*lambda_p(+1)*R_d*(1-tau_c(+1))*pi(+1)^(-1);
%Gamm =  exp(- k_tilde *(b_f*s - s_bar*bf_bar))*exp(eps_gamm);
%s = s_bar;

//Borrowers (Equation 6 - 16)

y=A*epz*l;
w = mc * A * epz;
lambda_i - (1-a_i)*(c_i - a_i*c_i(-1))^(-1)=0;
lambda_i*(1-tau_w)*w + mu_dti*phi_dti*(pi(+1))*w(+1)*(kappa + r)^(-1) + lambda_i*(-(1+tau_h)*h +(1-delta_h)*h(-1))*(q_h-q_l)*(lambda_bar*exp(-((1+r)-(1+r_bar)))*exp(y-y_bar)*epz*A*exp(eps_lambda))=l^theta;
beta_i*lambda_i(+1)*q(+1)*(1-delta_h) + mu_ltv*phi_ltv*q(+1)*pi(+1) = lambda_i*(1+tau_h)*q;
- exp(-theta_tilde*(L-L_bar))*(1-kappa)+ theta_tilde*exp(-theta_tilde*(L-L_bar))*(1-kappa)*L + lambda_i - mu_dti - mu_ltv*((R_m^4) -1)*100 - beta_i*(lambda_i(+1)*(R_m + kappa)*(pi(+1))^(-1)) = 0;
%((R_m(-1)^4) - 1)*100*L(-1) = phi_ltv(-1)*q*h(-1)*pi;
((R_m^4) - 1)*100*L = phi_ltv*q(+1)*h*pi(+1);
%L* (kappa + r) - phi_dti*pi(+1)*w(+1)*l = 0;
L(-1)* (kappa + r) - phi_dti(-1)*pi*w*l(-1) = 0;
%phi_ltv=(1-rho_ltv)*phi_ltv_bar+rho_ltv*phi_ltv(-1)- rholtv*(D-D_bar)- eta_ltv;
%phi_dti=(1-rho_dti)*phi_dti_bar+rho_dti*phi_dti(-1) - rhod*(D-D_bar) - eta_dti;
phi_ltv=(1-rho_ltv)*phi_ltv_bar+rho_ltv*phi_ltv(-1)- eta_ltv;
phi_dti=(1-rho_dti)*phi_dti_bar+rho_dti*phi_dti(-1) - eta_dti;
D =(1-kappa)*D(-1) + L;

//Banks (17-24)
L=d+k_b;
k_b=(1-delta_k)*k_b(-1) + omega*(1-kappa)*L;
npa = (1-delta_w)*npa(-1) + f * (1-kappa)*L + 0.5*(1-(npa/npa(-1)))^2*npa;
delta_w = rho_w*delta_w(-1) + (1-rho_w)*delta_w_bar + eps_w;

kk*(k_b/L - nu_b)*(k_b/L)^2 + beta_p*(R_m + kappa - R_d) + mu_b*omega*(1-kappa) + mu_n*(1-kappa)*f = 0;

beta_p  *(R_d - 1) + beta_p*mu_b(+1)*(1-delta_k) = kk*(k_b/L - nu_b)*(k_b/L) + 0.5*kk*(k_b/L - nu_b)^2 + mu_b;

-s_tilde*exp(-(npa/npa_bar)^2) + 2*s_tilde*exp(-(npa/npa_bar)^2)*(npa/npa_bar)^2 + mu_n*(-1 + 0.5*(1-(npa/npa(-1)))^2 + (1- (npa/npa(-1)))*(npa/npa(-1))) + beta_p*mu_n(+1)*((1-delta_w(+1)) + (1 - (npa(+1)/npa))*(npa(+1)/npa)^2) = 0;

%k_b = k_b_bar;

//Government (25 - 29)
%tau_h*q*h + tau_w*w*l + B + tau_c*R_d*d + tau_b*R_star*b_f = tr + (1+r)*B(-1);
tau_h*q*h + tau_w*w*l + B + tau_c*R_d*d = tr + (1+r)*B(-1);
tau_w - tauw_bar = rho_w*(tau_w(-1)-tauw_bar) + phi_wb*(B(-1)/y(-1) - B_bar/y_bar) +  epsilon_w; % + phi_wb*(B(-1)/y(-1) - B_bar/y_bar) + epsilon_w;
tau_h - tauh_bar = rho_h*(tau_h(-1)-tauh_bar) + phi_hb*(B(-1)/y(-1) - B_bar/y_bar) + epsilon_h; %  + epsilon_h;
tau_c - tauc_bar = rho_c*(tau_c(-1)-tauc_bar) + phi_cb*(B(-1)/y(-1) - B_bar/y_bar) +  epsilon_c;
%tau_b - taub_bar = rho_b*(tau_w(-1)-tauw_bar) + phi_bb*(B(-1)/y(-1) - B_bar/y_bar) +  epsilon_b;

//Housing (30-33)

lambda = lambda_bar*exp(-((1+r)-(1+r_bar)))*exp((y-y_bar))*exp(eps_lambda);
%lambda = lambda_bar*exp(0.005*(y-y_bar))*exp(eps_lambda);
q = lambda*q_h + (1-lambda)*q_l + eps_q;

q_h = rhoq*q_h(-1)+(1-rhoq)*q_h_bar + (epz-epz_bar) ;

q_l = rhol*q_l(-1)+(1-rhol)*q_l_bar - (epz-epz_bar);


//Central Bank (34)

%1+r = phir*(1+r_bar) + (1-phir)*(1+r(-1)) + (1-phir)*(neta*(pi-pi_bar) + zeta*(y-y_bar)) + epsmss;
%1+r = (1+r_bar)^(1-phir) *(1+r(-1))^(phir)*((pi/pi_bar)^neta)^(phir)*((y/y_bar)^zeta)^(phir);
1+r = (1 + r_bar)^phir * (1 + r(-1))^(1-phir) * (pi/pi_bar)^((1-phir)*neta)*epsm ;
%1+r = (1-phir)*(1 + r_bar) + phir*(1 + r(-1)) + (1-phir)*(neta)*(pi-pi_bar) + epsmss;

//prices (35)
%sdf = beta_p*((c_p-a_p*c_p(-1))/(c_p(-1)-a_p*c_p(-2)))/pi;
y - epsilon*y + epsilon*mc*y - psi*(pi-1)*pi*y + psi*(beta_p*((c_p(+1)-a_p*c_p)/(c_p-a_p*c_p(-1)))/pi)*(pi(+1)-1)*y(+1)*pi(+1)^2=0;

//Market closing equations (36-40)
h=h_bar;
%y = c + k_b - k_b(-1) + h - (1-delta_h)*h(-1) + R_l_star*s*k_b + s_tilde*exp(-(npa/npa_bar)^2);
%y = c_p + c_i + s_tilde*exp(-(npa/npa_bar)^2)*npa_bar + q*(h - (1-delta_h)*h(-1)) ;
y = c_p + c_i + s_tilde*exp(-(npa/npa_bar)^2)*npa_bar + q*(h - (1-delta_h)*h(-1)) ;
%y = c_p + c_i + s_tilde*exp(-(npa/npa_bar)^2) + q*(h- (1-delta_h)*h(-1));
B = 0;
%R_star = R_star_bar;
%R_l_star = R_l_star_bar;
%l = 0.33;

//Exogenous Stochastic Processes
%eps = rhos*eps(-1)+(1-rhos)*eps_bar + epss;
epsm = rhom*epsm(-1) + (1-rhom)*epsm_bar + epsmss;
epz = rhoz*epz(-1) + (1-rhoz)*epz_bar + epzss;
%epsq = rhoqh*epsq(-1) + (1-rhoqh)*epsq_bar + eps_q;
%epsilon_h = rhoh*epsilon_h(-1) + (1-rhoh)*epsilon_h_bar + epsh;
end;

initval;
y = 1;
B = 0;
%b_f = 0;
q = q_bar;
q_h = q_h_bar;
q_l = q_l_bar;
h=h_bar;
lambda_i = 1.223553996;
c_i = 0.817291271; 
c_p = 0.104166905;
lambda_p = 9.599978035;
%c   = c_p+c_i;
l   = 0.33;
w   = 0.25/0.33;
npa = 0.099961;
delta_w = 0.028784;
mu_n  = 15.51457458;
r        = 0.016125; %0.0136; %0.01357871;
R_d      = 1.35483333; 
%R_star   = 1.2701561073;
%R_l_star   = 1.276523405;
%s        = 1;
pi       = pi_bar;
tau_w    = 0.19;
tau_h    = 0.000624;
tau_c    = 0.25;
%tau_b    = 0.20; 
mu_ltv     = 0.055489742; 
mu_dti    = 0.305869394;
D        = 5.3596969697;  %1.35355;
L        = 0.176870; %0.04466715
phi_dti  = 0.034754955; %2.7071
phi_ltv  = 0.78;
R_m      = 1.010847272; %1.041058
k_b      = 0.030861;
d        = 0.146009;
tr       = 0.097578465;
mu_b     = 2.93190386;
%Gamm     = 1;
epsm     = epsm_bar;
epz      = 1;
mc       = 0.25;
y        = y_bar;
lambda = lambda_bar;
npa = 0.099961;
%epsq = 1;
end;

shocks;
var epzss; 
stderr 0.05;
var epsmss; 
stderr 0.05;
var eta_ltv; 
stderr 0.05;
var eta_dti; 
stderr 0.05;
var eps_w; 
stderr 0.001;
var eps_lambda; 
stderr 0.05;
var eps_q; 
stderr 0.01;
end;

steady;

check(qz_zero_threshold=1e-17);
model_diagnostics;
stoch_simul(order=1,irf=20) epsm y c_p c_i l w phi_ltv lambda q q_h q_l h r R_d R_m npa D L d k_b mc pi;
%stoch_simul(order=1,irf=20);
