function F = ss_fun(logit_varrho,M_)

%% Load model parameters

NumberOfParameters = M_.param_nbr;                            % Number of deep parameters.
for i = 1:NumberOfParameters                                  % Loop...
  paramname = deblank(M_.param_names(i,:));                  %    Get the name of parameter i.
  if strcmp(paramname,'xx')==0
  eval([ paramname ' = M_.params(' int2str(i) ');'])         %    Get the value of parameter i.
  end
end        
%%
varrho=exp(logit_varrho)/(1+exp(logit_varrho));
Pi = Pi_bar;
PiF = Pi_bar;
Pi_star = Pi_star_bar;
PiS = PiF/Pi_star;
phi = ((Pi/betta) / R_star_bar) / PiS;
J_by_JJ   = ( (1 - xi*Pi^(zzeta-1))/(1-xi) )^(1/(1-zzeta));
MC = J_by_JJ * ((zzeta-1)/zzeta) * ( 1 - xi*betta * Pi^zzeta ) / ( 1 - xi * betta * Pi^(zzeta-1) ); 
Disp = ( (1-xi)*J_by_JJ  ^(-zzeta) ) / ( 1- xi * Pi^zzeta );
h_by_K = ( ( 1/betta - (1-delta) ) / ( MC*(1-alp) ) )^(1/alp);
WP = alp * MC * h_by_K^(alp-1);
YW = H_bar*( h_by_K )^(alp-1);
Y = YW/Disp;
BF = -Y * log(phi) / phiB;
TB = BF * ( 1/(phi*R_star_bar) - PiS/Pi );
C = Y+YO -  delta*H_bar/h_by_K -  gy*Y - TB;
EX = Y - delta*wI*H_bar/h_by_K - C*wC - gy*Y;
Gov_bar = gy*Y;
HRoT = (1-varrho)/(1-varrho*chi);
C = ( YW/Disp - wI*delta*YW*( h_by_K )^(-alp) - EX - Gov_bar )/wC;
CR = C/(1-lambda) - HRoT*WP*lambda/(1-lambda);
HR = H_bar/(1-lambda) - HRoT*lambda/(1-lambda);

F = ( varrho / (1-varrho) )*(1-chi)*CR/(1-HR) - WP;
