function [Y_b,X_b,C_b,Ny_b,Nx_b,K_b,Iy_b,Ix_b,Z_b,Px_b,W_b,m_b,e_b,Rz_b,B_b,Rf_b] = ss_AI(alpha,sigma,A,deltak,M,gamma,deltaz,eta,emin,E,sigmaz,sigmab,deltae,thetas,beta,psi,phi,B,C,tauss)

%% solve for the bubbly steady state
options = optimset('TolFun',1e-9);
e_b = fsolve(@(x)solve_estar_bubbly(x,eta,beta,deltae,sigmab,thetas),2);
L_b = (1/(eta-1)*(e_b/emin)^(-eta));
omega_b = (eta*e_b/(eta-1));
Rz_b = ((1-beta*(1-deltae)*(1-deltaz+sigmaz*L_b))/(beta*(1-deltae)*e_b*(1+L_b)));
Rf_b = (1/(beta*(1-deltae)*(1+L_b)));

% solve for cutoff efficiency m

m_b = fsolve(@(x)solve_m(x,C,gamma,A,sigma,Rz_b,tauss),1,options);
inte_b = 1/(A*(sigma-1))*(exp(A*(sigma-1)*m_b)-1);
W_b = B*(1-m_b+inte_b)^(1/(sigma-1))*exp(A*(M-m_b));
Px_b = 1/(1+tauss)*W_b/exp(A*(M-m_b));

% solve for Ny and then the whole steady state

Ny_b = fsolve(@(x)solve_ss(x,M,A,sigma,alpha,beta,deltak,deltaz,Rz_b,omega_b,gamma,psi,phi,E,W_b,Px_b,m_b),0.5,options);
inte_b = 1/(A*(sigma-1))*(exp(A*(sigma-1)*m_b)-1);
X_b = exp(A*(M-m_b))*(1-m_b)/inte_b*Ny_b;
K_b = (alpha*beta/(1-beta*(1-deltak)))^(1/(1-alpha))*(1-m_b+inte_b)^(sigma/(sigma-1))/inte_b*Ny_b*exp(A*(M-m_b));
Y_b = K_b^(alpha)*((1-m_b)^(1/sigma)*X_b^((sigma-1)/sigma)+inte_b^(1/sigma)*(exp(A*(M-m_b))*Ny_b)^((sigma-1)/sigma))^(sigma*(1-alpha)/(sigma-1));
Iy_b = deltak*K_b;
Ix_b = deltaz/omega_b*gamma/Rz_b*Px_b*X_b;
C_b = Y_b-Ix_b-Iy_b;
Z_b = gamma*Px_b*X_b/Rz_b;
N_b = (W_b/(psi*C_b))^(1/phi);
Nx_b = N_b-Ny_b;
B_b = (Ix_b/((e_b/emin)^(-eta))-gamma*Px_b*X_b-sigmaz/e_b*Z_b)/sigmab;


%% functions
function dis = solve_estar_bubbly(x,eta,beta,deltae,sigmab,thetas)
    emin = (eta-1)/eta;
    L = (1/(eta-1)*(x/emin)^(-eta));
    LHS = 1;
    RHS = beta*(1-deltae)*thetas*(1+sigmab*L);
    dis = LHS-RHS;
end

function dis = solve_m(m,C,gamma,A,sigma,Rz,tauss)
    inte = 1/(A*(sigma-1))*(exp(A*(sigma-1)*m)-1);
    RHS = C*exp((1-gamma)/gamma*A*m)*(1-m+inte)^(1/(sigma-1))*(1+tauss)^(-1/gamma);
    dis = Rz-RHS;
end

function dis = solve_ss(Ny,M,A,sigma,alpha,beta,deltak,deltaz,Rz,omega,gamma,psi,phi,E,W,Px,m)
    inte = 1/(A*(sigma-1))*(exp(A*(sigma-1)*m)-1);
    X = exp(A*(M-m))*(1-m)/inte*Ny;
    K = (alpha*beta/(1-beta*(1-deltak)))^(1/(1-alpha))*(1-m+inte)^(sigma/(sigma-1))/inte*Ny*exp(A*(M-m));
    Y = K^(alpha)*((1-m)^(1/sigma)*X^((sigma-1)/sigma)+inte^(1/sigma)*(exp(A*(M-m))*Ny)^((sigma-1)/sigma))^(sigma*(1-alpha)/(sigma-1));
    Iy = deltak*K;
    Ix = deltaz/omega*gamma/Rz*Px*X;
    C = Y-Ix-Iy;
    Z = gamma*Px*X/Rz;
    N = (W/(psi*C))^(1/phi);
    Nx = N-Ny;
    dis = X-E*Z^(gamma)*Nx^(1-gamma);
end
end
