load param

spB= 0.01;
sp= 0.05;
Lams = beta;
Rds = Pis/beta;
RFs = (1+sp)^(1/4)*Rds;
%RBs = Deltas*(RFs - Rds) + Rds;
RBs = (1+spB)^(1/4)*Rds;

QFs = (RFs - kappaf)^(-1);
QBs = (RBs - kappab)^(-1);
RLFs = QFs^(-1) + kappaf;
RLBs = QBs^(-1) + kappab;
QEHs = (Rds - kappab)^(-1);
REHs = QEHs^(-1) + kappab;
M2s = Lams*Pis^(-1)/(QFs*(1-kappaf*Lams*Pis^(-1)));
Ms = (M2s-1)*psi+1;
ps = ( (1 - phip*Pis^((1-epsip)*(gammap-1)))      / (1-phip)         )^(1/(1-epsip)); % steady state reset price
pms = ((epsip-1)/epsip)*(1-phip*Lams*Pis^(epsip*(1-gammap)))/(1-phip*Lams*Pis^((1-epsip)*(gammap-1))); % steady state marginal cost
Ks = (alpha*pms/(Ms/Lams - (1-delta0)*Ms))^(1/(1-alpha));
Yms = Ks^(alpha);
ws = pms*(1-alpha)*Ks^(alpha);
Ihs = delta0*Ks;
Is = Ihs;
vps = ((1-phip)*ps^(-epsip))/(1-phip*Pis^(epsip*(1-gammap)));
Ys = Yms/vps;
Gs = gy*Ys;
Cs = Ys - Is - Gs;
wss = (ws^(1-epsiw)*(1-phiw*Pis^((1-epsiw)*(gammaw-1)))/(1-phiw))^(1/(1-epsiw)); % reset wage
mrss = ((epsiw-1)/epsiw)*wss*(1-phiw*Lams*Pis^(epsiw*(1-gammaw)))/(1-phiw*Lams*Pis^((1-epsiw)*(gammaw-1)));
vws = (1-phiw)*(wss/ws)^(-epsiw)/(1-phiw*Pis^(epsiw*(1-gammaw)));
Ls = vws;
mus = (1/Cs)*(1/(1-bh) - beta*bh/(1-bh));
chi = mus*mrss*Ls^(-eta);
delta1 = pms*alpha*Ks^(alpha-1)/Ms;
x1s = pms*Ys/(1-phip*Lams*Pis^(epsip*(1-gammap)));
x2s = Ys/(1-phip*Lams*Pis^((1-epsip)*(gammap-1)));
f1s = mrss*ws^(epsiw)/(1-phiw*Lams*Pis^(epsiw*(1-gammaw)));
f2s = ws^(epsiw)/(1-phiw*Lams*Pis^((1-epsiw)*(gammaw-1)));
fprs = psi*Ihs/(QFs*(1-kappaf*Pis^(-1)));
bGs = bgy*Ys/QBs;
bcbs = Ys*fracb/QBs; 
fcbs = Ys*fracf/QFs;
bs = bGs - bcbs;
res = QFs*fcbs + QBs*bcbs;
fs = fprs - fcbs;
ns = (QFs*fs + QBs*bs + res)/levs;
ds = QFs*fs + QBs*bs + res - ns;
phis = (QFs*fs + Deltas*QBs*bs)/ns;

% adding steady state for SLR based model
gammas = (QBs*bs)/(QFs*fs);
A = (QFs*fs)*(RFs-Rds+gammas*(RBs-Rds)) + Rds*ns;
thetas = (Pis^(-1)*Lams*(1-sigma)*A)/(QFs*fs*(1+Deltas*gammas) - A*Pis^(-1)*Lams*sigma*phis);
B = (A - Rds*ns)/(QFs*fs);
lam_tilda = (ns*B*thetas*phis)/(A*(1+Deltas*gammas));
lams = lam_tilda/(thetas - lam_tilda);
zeta_tilda = ((thetas*phis*ns*(RFs-Rds)/A) - lam_tilda)/gammas;
zetas = zeta_tilda*(1+lams);
X = ns - sigma*Pis^(-1)*(QFs*fs*(RFs-Rds+gammas*(RBs-Rds)) + Rds*ns);        % after eqn 149 in Sims Notes
%%check the value of X

%thetas = ( (1-sigma)*Rds*Pis^(-1)*Lams + (1-sigma)*(RFs-Rds)*phis*Pis^(-1)*Lams)/(phis - sigma*(RFs - Rds)*phis^2*Pis^(-1)*Lams - sigma*phis*Rds*Pis^(-1)*Lams);
%X = ns*(1-sigma*((RFs-Rds)*phis + Rds)*Pis^(-1));
as = thetas*phis;
Omegas = 1-sigma + sigma*as;
%lams = ((thetas/((RFs - Rds)*Lams*Omegas*Pis^(-1))) - 1)^(-1);
Res = Rds;
omegas = (Rds - Res)*(1+lams)/Pis*Lams*Omegas;
Tcbs = (1+kappaf*QFs)*Pis^(-1)*fcbs + (1+kappab*QBs)*Pis^(-1)*bcbs - Res*Pis^(-1)*res;
Ts = Gs + bGs*Pis^(-1) - Tcbs - QBs*(bGs - kappab*Pis^(-1)*bGs);