load param

Lams = beta;
Rds = Pis/beta;
RFss = (1+sp)^(1/4)*Rds;
RFps = Deltas1*(RFss - Rds) + Rds;
RBs = Deltas2*(RFss - Rds) + Rds;
QFss = (RFss - kappaf)^(-1);
QFps = (RFps - kappaf)^(-1);
QBs = (RBs - kappab)^(-1);
RLFss = QFss^(-1) + kappaf;
RLFps = QFps^(-1) + kappaf;
RLBs = QBs^(-1) + kappab;
QEHs = (Rds - kappab)^(-1);
REHs = QEHs^(-1) + kappab;
M2ss = Lams*Pis^(-1)/(QFss*(1-kappaf*Lams*Pis^(-1)));
M1ss = (M2ss-1)*psi+1;
M2ps = Lams*Pis^(-1)/(QFps*(1-kappaf*Lams*Pis^(-1)));
M1ps = (M2ps-1)*1+1;

% sectoral outputs: we can observe ps*ys/(pp*yp)
rhoY = fsolve(@(x)solverhoY(x,ratio,sigmaY,phips,phipp,Pis,epsips,epsipp,gammaps,gammapp, ...
    Lams,gamma,xis,Ldss,M1ss,delta,alpha,Ldps,M1ps),0.5);
pss = (1/((1+ratio)*(1-rhoY)))^(1/(1-sigmaY));
pps = pss*((1-rhoY)/rhoY)^(1/(1-sigmaY))*ratio^(1/(1-sigmaY));
prss = ( (1 - phips*Pis^((1-epsips)*(gammaps-1)))*pss^(1-epsips)      / (1-phips)         )^(1/(1-epsips)); % steady state reset price SOE
prps = ( (1 - phipp*Pis^((1-epsipp)*(gammapp-1)))*pps^(1-epsipp)      / (1-phipp)         )^(1/(1-epsipp)); % steady state reset price POE
pmss = ((epsips-1)/epsips)*prss*(1-phips*Lams*Pis^(epsips*(1-gammaps)))/(1-phips*Lams*Pis^((1-epsips)*(gammaps-1))); % steady state marginal cost SOE
pmps = ((epsipp-1)/epsipp)*prps*(1-phipp*Lams*Pis^(epsipp*(1-gammapp)))/(1-phipp*Lams*Pis^((1-epsipp)*(gammapp-1))); % steady state marginal cost SOE
Kss = (gamma*pmss*xis*Ldss^(1-gamma)/(pks*M1ss/Lams - (1-delta)*M1ss*pks))^(1/(1-gamma));
Ymss = xis*Kss^(gamma)*Ldss^(1-gamma);
wss = xis*pmss*(1-gamma)*Kss^(gamma)*Ldss^(-gamma);
Ihss = delta*Kss;
Kps = (alpha*pmps*Ldps^(1-alpha)/(pks*M1ps/Lams - (1-delta)*M1ps*pks))^(1/(1-alpha));
Ymps = Kps^(alpha)*Ldps^(1-alpha);
wps = pmps*(1-alpha)*Kps^(alpha)*Ldps^(-alpha);
Ihps = delta*Kps;
I_s = Ihss+Ihps;

vpss = ((1-phips)*(prss/pss)^(-epsips))/(1-phips*Pis^(epsips*(1-gammaps)));
vpps = ((1-phipp)*(prps/pps)^(-epsipp))/(1-phipp*Pis^(epsipp*(1-gammapp)));
Yss = Ymss/vpss;
Yps = Ymps/vpps;
Y_s = (rhoY^(1/sigmaY)*Yps^((sigmaY-1)/sigmaY)+(1-rhoY)^(1/sigmaY)*Yss^((sigmaY-1)/sigmaY))^(sigmaY/(sigmaY-1));
Gs = gy*Y_s;
Cs = Y_s - I_s - Gs;

wrss = (wss^(1-epsiws)*(1-phiws*Pis^((1-epsiws)*(gammaws-1)))/(1-phiws))^(1/(1-epsiws)); % reset wage SOE
wrps = (wps^(1-epsiwp)*(1-phiwp*Pis^((1-epsiwp)*(gammawp-1)))/(1-phiwp))^(1/(1-epsiwp)); % reset wage POE
mrsss = ((epsiws-1)/epsiws)*wrss*(1-phiws*Lams*Pis^(epsiws*(1-gammaws)))/(1-phiws*Lams*Pis^((1-epsiws)*(gammaws-1)));
mrsps = ((epsiwp-1)/epsiwp)*wrps*(1-phiwp*Lams*Pis^(epsiwp*(1-gammawp)))/(1-phiwp*Lams*Pis^((1-epsiwp)*(gammawp-1)));
vwss = (1-phiws)*(wrss/wss)^(-epsiws)/(1-phiws*Pis^(epsiws*(1-gammaws)));
vwps = (1-phiwp)*(wrps/wps)^(-epsiwp)/(1-phiwp*Pis^(epsiwp*(1-gammawp)));
Lss = Ldss*vwss;
Lps = Ldps*vwps;
rhoL = 1/(1+Lss/Lps*(mrsss/mrsps)^sigmaL);
L_s = (rhoL^(1/sigmaL)*Lps^((sigmaL-1)/sigmaL)+(1-rhoL)^(1/sigmaL)*Lss^((sigmaL-1)/sigmaL))^(sigmaL/(sigmaL-1));
mus = (1/Cs)*(1/(1-bh) - beta*bh/(1-bh));
chi = mus*mrsps*L_s^(-eta)*rhoL^(-1/sigmaL)*(Lps/L_s)^(1/sigmaL);

x1ss = pss^(epsips)*pmss*Yss/(1-phips*Lams*Pis^(epsips*(1-gammaps)));
x2ss = pss^(epsips)*Yss/(1-phips*Lams*Pis^((1-epsips)*(gammaps-1)));
x1ps = pps^(epsipp)*pmps*Yps/(1-phipp*Lams*Pis^(epsipp*(1-gammapp)));
x2ps = pps^(epsipp)*Yps/(1-phipp*Lams*Pis^((1-epsipp)*(gammapp-1)));
f1ss = mrsss*wss^(epsiws)/(1-phiws*Lams*Pis^(epsiws*(1-gammaws)));
f2ss = wss^(epsiws)/(1-phiws*Lams*Pis^((1-epsiws)*(gammaws-1)));
f1ps = mrsps*wps^(epsiwp)/(1-phiwp*Lams*Pis^(epsiwp*(1-gammawp)));
f2ps = wps^(epsiwp)/(1-phiwp*Lams*Pis^((1-epsiwp)*(gammawp-1)));

fprss = psi*Ihss/(QFss*(1-kappaf*Pis^(-1)));
fprps = Ihps/(QFps*(1-kappaf*Pis^(-1)));
bGs = bgy*Y_s/QBs;
bcbs = Y_s*fracb/QBs; 
fcbss = Y_s*fracfs/QFss;
fcbps = Y_s*fracfp/QFps;
bs = bGs - bcbs;
res = QFss*fcbss + QFps*fcbps + QBs*bcbs;
fss = fprss - fcbss;
fps = fprps - fcbps;
ns = (QFss*fss + QFps*fps + QBs*bs + res)/levs;
ds = QFss*fss + QFps*fps + QBs*bs + res - ns;
phis = (QFss*fss + Deltas1*QFps*fps + Deltas2*QBs*bs)/ns;
thetas = ( (1-sigma)*Rds*Pis^(-1)*Lams + (1-sigma)*(RFss-Rds)*phis*Pis^(-1)*Lams)/(phis - sigma*(RFss - Rds)*phis^2*Pis^(-1)*Lams - sigma*phis*Rds*Pis^(-1)*Lams);
X = ns*(1-sigma*((RFss-Rds)*phis + Rds)*Pis^(-1));
as = thetas*phis;
Omegas = 1-sigma + sigma*as;
lams = ((thetas/((RFss - Rds)*Lams*Omegas*Pis^(-1))) - 1)^(-1);
Res = Rds;
omegas = (Rds - Res)*(1+lams)/Pis*Lams*Omegas;
Tcbs = (1+kappaf*QFss)*Pis^(-1)*fcbss + (1+kappaf*QFps)*Pis^(-1)*fcbps + (1+kappab*QBs)*Pis^(-1)*bcbs - Res*Pis^(-1)*res;
Ts = Gs + bGs*Pis^(-1) - Tcbs - QBs*(bGs - kappab*Pis^(-1)*bGs);

function dis = solverhoY(rhoY,ratio,sigmaY,phips,phipp,Pis,epsips,epsipp,gammaps,gammapp, ...
    Lams,gamma,xis,Ldss,M1ss,delta,alpha,Ldps,M1ps)
    % ratio = (pp*yp)/(ps*ys)
    pss = (1/((1+ratio)*(1-rhoY)))^(1/(1-sigmaY));
    pps = pss*((1-rhoY)/rhoY)^(1/(1-sigmaY))*ratio^(1/(1-sigmaY));
    prss = ( (1 - phips*Pis^((1-epsips)*(gammaps-1)))*pss^(1-epsips)      / (1-phips)         )^(1/(1-epsips)); % steady state reset price SOE
    prps = ( (1 - phipp*Pis^((1-epsipp)*(gammapp-1)))*pps^(1-epsipp)      / (1-phipp)         )^(1/(1-epsipp)); % steady state reset price POE
    pmss = ((epsips-1)/epsips)*prss*(1-phips*Lams*Pis^(epsips*(1-gammaps)))/(1-phips*Lams*Pis^((1-epsips)*(gammaps-1))); % steady state marginal cost SOE
    pmps = ((epsipp-1)/epsipp)*prps*(1-phipp*Lams*Pis^(epsipp*(1-gammapp)))/(1-phipp*Lams*Pis^((1-epsipp)*(gammapp-1))); % steady state marginal cost SOE
    Kss = (gamma*pmss*xis*Ldss^(1-gamma)/(M1ss/Lams - (1-delta)*M1ss))^(1/(1-gamma));
    Ymss = xis*Kss^(gamma)*Ldss^(1-gamma);
    wss = pmss*xis*(1-gamma)*Kss^(gamma)*Ldss^(-gamma);
    Ihss = delta*Kss;
    Kps = (alpha*pmps*Ldps^(1-alpha)/(M1ps/Lams - (1-delta)*M1ps))^(1/(1-alpha));
    Ymps = Kps^(alpha)*Ldps^(1-alpha);
    wps = pmps*(1-alpha)*Kps^(alpha)*Ldps^(-alpha);
    Ihps = delta*Kps;
    Is = Ihss+Ihps;
    vpss = ((1-phips)*(prss/pss)^(-epsips))/(1-phips*Pis^(epsips*(1-gammaps)));
    vpps = ((1-phipp)*(prps/pps)^(-epsipp))/(1-phipp*Pis^(epsipp*(1-gammapp)));
    Yss = Ymss/vpss;
    Yps = Ymps/vpps;
    dis = ratio-pps*Yps/(pss*Yss);
end

