%% steady state calculation

clear all
clc

% Deep parameters
betal = 0.995; % discount factor of lender
betab = 0.97; % discount factor of borrower
eta = 1; % frisch labor supply elasticity
phih = 0.05; % housing adjustment cost
deltah = 0.01; % housing depreciation rate
j = 0.12; % relative housing preference
chil = 7; % disutility parameter for lender
chib = 9; % disutility parameter for borrower
m = 0.65; % LTV ratio
alphal = 0.6; % labor income share of lender
theta = 0.83; % price stickiness parameter
epsilon = 11; % elasticity of substitution btw differentiated goods
gammar = 0.5; % interest rate smoothing parameter
gammapi = 1.5; % inflation responsiveness
gammay = 0.5; % output responsiveness
hbar = 1; % total housing stock(normalized)
rhou = 0.9; % persistence of costpush shock
rhoz = 0.9; % persistence of demand shock
su = 0.01; % stdev of costpush shock innovation
sz = 0.01; % stdev of demand shock innovation
uss = 1;
zss = 1;

% Steady state
piss = 1; %steady state inflation rate
Rss = piss/betal;
rss = Rss/piss;
rrbar = Rss;
pirss = ((1-theta*piss^(epsilon-1))/(1-theta))^(1/(1-epsilon));
vpss = ((1-theta)*pirss^(-epsilon))/(1-theta*piss^(epsilon));
mcss = pirss*((epsilon-1)/epsilon)*((1-theta*betal*piss^(epsilon))/(1-theta*betal*piss^(epsilon-1)));

% Auxiliary variables
sig1 = j/(1-betal*(1-deltah));
sig2 = j/(1-betab*(1-deltah)-m*(betal-betab));
cb_y = ((1-alphal)*mcss*vpss)/(1+deltah*sig2+(1-betal)*m*sig2);
cl_y = (alphal*mcss*vpss + (1-mcss*vpss) - (betal-1)*m*sig2*cb_y)/(1+deltah*sig1);

% Steady state
Llss = (alphal*mcss*vpss/(chil*cl_y))^(1/(1+eta));
Lbss = ((1-alphal)*mcss*vpss/(chib*cb_y))^(1/(1+eta));
yss = ((Llss^(alphal))*(Lbss^(1-alphal)))/vpss;
wlss = alphal*mcss*vpss*yss/Llss;
wbss = (1-alphal)*mcss*vpss*yss/Lbss;
clss = cl_y*yss;
cbss = cb_y*yss;
qss = (sig1*clss + sig2*cbss)/hbar;
hlss = sig1*clss/qss;
hbss = sig2*cbss/qss;
bbss = betal*m*qss*hbss;
mubss = (betal-betab)/(cbss*piss);
sdfss = betal;
z1ss = (mcss*yss)/(1-theta*sdfss*(piss^(epsilon)));
z2ss = yss/(1-theta*sdfss*(piss^(epsilon-1)));
Vlss = (log(clss)+j*log(hlss)-((Llss)^(1+eta))/(1+eta))/(1-betal);
Vbss = (log(cbss)+j*log(hbss)-((Lbss)^(1+eta))/(1+eta))/(1-betab);
omegabss = (qss*(1-deltah)*hbss - Rss*bbss/piss)/(qss*(1-deltah)*hbar);

pibar = piss;
ybar = ((Llss^(alphal))*(Lbss^(1-alphal)))/vpss;
ubar = 1;
zbar = 1;

save tank_costhousing_param_base betal betab alphal eta theta epsilon phih chil chib gammar gammapi gammay hbar rrbar deltah rhou rhoz su sz clss cbss hlss hbss Llss Lbss wlss wbss Rss rss bbss piss pirss vpss yss z1ss z2ss j m qss uss zss mubss mcss sdfss Vlss Vbss pibar ybar ubar zbar;

dynareOBC tank_consthousing.mod

% Wcl = (exp((1-betal)*(mean(Vl)-Vlss))-1)*100
% Wcb = (exp((1-betab)*(mean(Vb)-Vbss))-1)*100

