function y = F(x, b, bbeta, vvarphi, ssigma, oomega_T, oomega_N, oomega_C, cchi, ddelta, aalpha_T, aalpha_N, aalpha_C, ggamma_T, pphi_T, pphi_N, pphi_C, v_C, v_D, rho_AC, rho_AT, rho_AN,  Dss, r, mu_Tss, mu_Nss, mu_Css, TBss, P_Css, KL_C, W_Css, KL_T, W_Tss, eeta_N, eeta_T, eeta_C, A_Nss, A_Tss, A_Css, theta)

reer = x(1);
KL_N = x(2);
L_Nss = x(3);
Pss_N = x(4);

% Nontradable sector: {solve numerically for reer, Kss_N/L_Nss, L_Nss and Pss_N}
y(1) = reer - (cchi^vvarphi + ((1 - cchi)^vvarphi)*(Pss_N^(1 - vvarphi)))^(1/(1 - vvarphi));

y(2) = mu_Nss - ((aalpha_N*Pss_N*A_Nss*(KL_N^(aalpha_N - 1)))/(1 + eeta_N*(r/(1+r))));

y(3) = (1 + eeta_N*(r/(1+r)))*(reer/(1 - b*bbeta))*(L_Nss^(oomega_N - 1)) - (1 - aalpha_N)*Pss_N*A_Nss*(KL_N^aalpha_N);

y(4) = Pss_N - ((1 - cchi)/cchi)*((((theta*A_Tss*(KL_T^(aalpha_T + ggamma_T))*((((1 - b*bbeta)/reer)*W_Tss)^(1/(oomega_T - 1)))) - (ddelta*KL_T*((((1 - b*bbeta)/reer)*W_Tss)^(1/(oomega_T - 1)))) - (ddelta*KL_N*L_Nss) - (ddelta*KL_C*((((1 - b*bbeta)/reer)*W_Css)^(1/(oomega_C - 1)))) - (TBss - (P_Css*((A_Css*(KL_C^aalpha_C)*((((1 - b*bbeta)/reer)*W_Css)^(1/(oomega_C - 1)))) - ((ggamma_T*(theta*A_Tss*(KL_T^(aalpha_T + ggamma_T))*((((1 - b*bbeta)/reer)*W_Tss)^(1/(oomega_T - 1)))))/(P_Css*(1 + eeta_T*(r/(1 + r)))))))))/(A_Nss*(KL_N^aalpha_N)*L_Nss))^(1/vvarphi));


% C_T = (Y_Tss) - (I_Tss) - (I_Nss) - (I_Css) - (TB_Tss)

% C_N = Y_Nss

% Y_Tss = theta*A_Tss*(KL_T^(aalpha_T + ggamma_T))*(L_Tss)

% I_Tss = ddelta*KL_T*(L_Tss)

% I_Css = ddelta*KL_C*(L_Css)

% I_Nss = ddelta*KL_N*L_Nss

% L_Tss = (((1 - b*bbeta)/reer)*W_Tss)^(1/(oomega_T - 1))

% L_Css = (((1 - b*bbeta)/reer)*W_Css)^(1/(oomega_C - 1))

% TB_Tss = TBss - TB_Css

% TB_Css = Pss_C*((Y_Css) - (CM_Tss))

% Y_Css = A_Css*(KL_C^aalpha_C)*(L_Css)

% CM_Tss = (ggamma_T*(Y_Tss))/(Pss_C*(1 + eeta_T*(r/(1 + r))))


end