%% Model parameters

% Household sector
gamma=0.6;
rho=3;
sigmal=1;
beta=0.99;
hb=0.56;
phib=1;

%production & capital markets
mu_w=1.2;
thetaw=0.6;
comega=0.6;
alpha=0.6;
delta=0.3;
tau=1;
epsil=0.9;
chifa=0.5;
chifastar=0.4;
iota=0.8;
%foreihn economics
chi=1;
omega=0.5;
%capital production
gammai=0.6;
rho_i=10;
%homegoods retailer
gammah=0.5;
rhopi=0.4;
thetah=0.6;
mu_hss=1.2;
%imported goods prodecers
gammaf=0.5;
thetaf=0.6;
mu_fss=1.1;
%monetary policy
rho_m=0.3;
taupi=1.5;
rhopi=0.4;
tauy=0.5;
zss=1;

%% assumption

ratioleverage=5;
ratioadjustcost=0.2;
ratioc=1;

%% parameters for shocks
rhoc = 0.8;
rhol = 0.8;
rhophi= 0.8;
rhozeta= 0.8;
rhon = 0.8;
rhoi= 0.9;
rhomuh = 0.2;
rho1muf   = 0.6;
rhopi  = 0.8;
rhoystar   = 0.9;
rhopstar    = 0.8;
rhop_hstar=0.8;
rhoistar   = 0.8;
rhos=0.7;
rhou=0.6;

stdl_xi_c = 0.017;
std1_xi_l      = 0.083;
std1_xi_phi   = 0.003;
std1_zeta    = 0.032;
std1_xi_n     = 0.021;
std1_xi_i = 0.007;
std1_mu_h   = 0.005;
std1_mu_f   = 0.005;
std1_pi_f  = 0.050;
std1_ystar  = 0.018;
std1_epsil_m    = 0.133;
stdl_pstar=0.056;
stdl_istar=0.104;
stdl_s=0.043;
stdl_u=0.023;
stdl_p_hstar=0.256;


epsil_css=0;
epsil_lss=0;
epsil_phiss=0;
epsil_zetass=0;
epsil_nss=0;
epsil_iss=0;
epsil_muhss=0;
epsil_mufss=0;
epsil_piss=0;
epsil_ystarss=0;
epsil_mss=0;
epsil_pstarss=0;
epsil_istarss=0;
epsil_sss=0;
epsil_uss=0;
epsil_p_fstarss=0;

%% -------------- COMPUTE THE SS for some endogenous variables -------------------------------

lambdass=1-beta*hb;
niss=1/beta;
piss=1;
pi_hss=1;
pi_fss=1;

zetass=1;
xi_css=1;
xi_lss=1;
xi_phi=1;
xi_nss=1;
xi_iss=1;
pstarss=1;
ystarss=2;
istarss=0.02;
sss=6;
uss=0.7;
p_fstarss=1;

p_fss=sss*mu_fss*p_fstarss;
pw_fss=sss*p_fstarss;

kss= (32*epsil_iss^2*uss^(2*epsil) + 64*epsil*epsil_iss^2*uss^(2*epsil) + 32*epsil^2*epsil_iss^2*uss^(2*epsil))/(125*tau^2*uss^2 + 125*delta^2*uss^(2*epsil) + 250*delta^2*epsil*uss^(2*epsil) + 125*delta^2*epsil^2*uss^(2*epsil) + 250*delta*tau*uss*uss^epsil + 250*delta*epsil*tau*uss*uss^epsil);
iss=(0.4*kss)^(1/2);
ratioadjustcost=iss^2/2/kss;

p_hss=ratioc^(-1/rho)*(gamma/1-gamma)^(1/rho)*p_fss;
p_hstarss=p_hss/sss;
pss=(gamma*p_hss^(1-rho)+(1-gamma)*(p_fss)^(1-rho))^(1/(1-rho));
pwss=p_hss/mu_hss;
p_iss=(gammai*p_hss^(1-rhoi)+(1-gammai)*(p_fss)^(1-rhoi))^(1/(1-rhoi));
qss=(3*kss*p_iss+3*beta*kss*p_iss-3*beta*(delta+tau/(1+epsil)*uss^(1-epsil))*kss*p_iss)/2/(beta*(delta+tau/(1+epsil)*uss^(1-epsil))-beta+1)/(9*xi_iss*pss*p_iss^2/4-3*xi_iss*kss*pss/2);
nss=qss*kss/ratioleverage;
ratioleverage=qss*kss/nss;
yss=tau*(uss^(-epsil))*kss*p_iss/pwss*uss/alpha;
r_kss=(pwss/pss*alpha*yss/kss-p_iss/pss*(delta+tau/(1+epsil)*uss^(1-epsil)+qss))/qss;
gyss=pwss/pss*yss+(qss-p_iss/pss*(delta+tau/(1+epsil)*uss^(1-epsil)))*kss;
lss=((1-alpha)*(1-omega)*yss*lambdass*pwss/xi_lss/pss)^(1/(1+sigmal));
hss=(yss/(uss*kss)^alpha)^(1/(1-alpha))/zss;
less=(hss/lss^comega)^(1/(1-comega));
wss=(1-alpha)*(1-comega)*yss/lss*pwss;
w_ess=(1-alpha)*comega*yss/less*pwss;

x0=[0.5;0.5];
x= fsolve(@(x) functioni(x,iss,gammai,rhoi,p_hss,p_fss),x0, ...
  optimset('TolFun',1e-6,'MaxIter',1000,'MaxFunEvals',10000000));

i_hss=x(1);
i_fss=x(2);

c_hstarss=(p_hstarss/pstarss)^(-chi)*ystarss;
vss=(nss-w_ess/pss)/iota*xi_nss;
c_ess=(1-iota*xi_nss)*vss;


xx0=[1/3;1/3/1.2];
xx= fsolve(@(xx) functionce(xx,c_ess,gamma,rho,p_hss,p_fss),xx0, ...
  optimset('TolFun',1e-6,'MaxIter',1000,'MaxFunEvals',10000000));

c_ehss=xx(1);
c_efss=xx(2);
c_hss=yss-c_ehss-i_hss-c_hstarss;
c_fss=c_hss/(gamma/(1-gamma)*(p_hss/p_fss)^(-rho));

ratioc=c_hss/c_fss;

css=(gamma^(1/rho)*c_hss^((rho-1)/rho)+(1-gamma)^(1/rho)*c_fsss^((rho-1)/rho))^(rho/(rho-1));
d_fss=3*css;
bss=(r_kss*qss*kss-chifastar*ratioleverage*istarss*sss*d_fss/pss)*pss/(chifa*ratioleverage*iss);
y_fss=c_fss+c_efss+i_fss;


xxx0=1/2;
xxx= fsolve(@(xxx) funbstar(xxx, wss,pss,lss, yss ,p_hss, pwss, y_fss, p_fss, pw_fss, bss, niss, css, sss, istarss, phib, xi_phiss),xxx0, ...
  optimset('TolFun',1e-6,'MaxIter',1000,'MaxFunEvals',10000000));

bstarss=xxx;



save steadystate.m