function res=steady_stateHH(X, parameters)
global Cr Cn C Nr Nn L W R Rk Lambda K MC Gamma_1 Gamma_2 Y Pstar_P Pi sp G B_P Pr;

Cr = X(1); 
Cn = X(2);
C = X(3);
Nr = X(4);
Nn = X(5);
L = X(6);
W = X(7);
R = X(8);
Rk = X(9);
Lambda = X(10);
K = X(11);
MC = X(12);
Gamma_1 = X(13);
Gamma_2 = X(14);
Y = X(15);
Pstar_P = X(16);
Pi = X(17);
sp = X(18);
G = X(19);
B_P = X(20);
Pr = X(21);



delta = parameters(1);
eta = parameters(2);
zeta = parameters(3);
lambda = parameters(4);
theta_p = parameters(5);
theta_w = parameters(6);
PHI_w = parameters(7);
alpha = parameters(8);
beta = parameters(9);
rho_R = parameters(10);
phi_pi = parameters(11);
phi_y = parameters(12);
rho_a = parameters(13);
rho_r = parameters(14);
rho_l = parameters(15);
rho_c = parameters(16);
sigma_a = parameters(17);
sigma_g = parameters(18);
sigma_r = parameters(19);
sigma_l = parameters(20);
sigma_c = parameters(21);
GY_ss = parameters(22);
psi = parameters(23);
k = parameters(24);
upsilon = parameters(25);
tau_w = parameters(26);
rho_b = parameters(27);
rho_g = parameters(28);
tau_c = parameters(29);
tau_k = parameters(30);
tau_p = parameters(31);
nu = parameters(32);



%% Steady State

res(1,1) = Nr - ( (1-tau_w)*W*(Cr*(1+tau_c))^(-1) )^(1/upsilon);   %                  
res(2,1) = R - ((1-tau_w)*beta)^(-1);   % differential in the expected real net returns
res(3,1) = Rk- ((1-tau_p)*beta)^(-1) - 1 + delta;
res(4,1) = Lambda - beta;
res(5,1) = Nn - ( (1-tau_w)*W*(Cn*(1+tau_c))^(-1) )^(1/upsilon);
res(6,1) = K - (alpha*W*L/((1-alpha)*Rk));
res(7,1) = MC - (zeta - 1)/zeta;
res(8,1) = Gamma_1 - MC*Y/(1 - theta_p*beta);
res(9,1) = Gamma_2 - Y/(1 - theta_p*beta);
res(10,1) = Pstar_P - 1;
res(11,1) = Pi - 1;
res(12,1) = sp - 1;
res(13,1) = G - GY_ss*Y;
res(14,1) = G + B_P*(1-1/R) - tau_w*W*L - tau_c*C - tau_k*Pr - tau_p*(1-delta+Rk)*K;
res(15,1) = C - Y + delta*K + G;
res(16,1) = Y - (K^alpha)*(L^(1-alpha));
res(17,1) = L - (((1-alpha)*Rk/(alpha*W))^alpha)*Y;
res(18,1) = L - nu*Nn - (1-nu)*Nr;
res(19,1) = C - nu*Cn - (1-nu)*Cr;
res(20,1) = Nn - (1+tau_w)*Cn/((1-tau_w)*W);
res(21,1) = Pr - Y + L*W + Rk*K;

% G + B_P*(1-1/R) - tau_w*W*L - tau_c*C - tau_k*Pr - tau_p*(1-delta+Rk)*K;
% (1+tau_c)*Cr + B_P*(1-R) - (1-tau_w)*W*Nr + (1 - (1-tau_p)*(Rk+1-delta))*K - (1-tau_k)*Pr
res=res'



% X = [Cr Cn C Nr Nn L W R Rk Lambda K MC Gamma_1 Gamma_2 Y Pstar_P Pi sp G B_P Pr]';
end