function out = fun_sop(c_star,s, nu, sigma,beta);

delta = 0.9998;

%  Inputs: c*, s, nu, sigma,beta
%  c* = beta*( s + theta)
%  This function's output is: beta*(s + theta) - c*;
%  To solve for c*, we will set this function value to zero in Matlab's Fzero function

% Ec = expected value of c conditional on c<c*, where c comes from a normal distribution ~ N(nu, sigma)
vector = [];
for i=1:5000
   r=normrnd(nu, sigma);
   if r<c_star
      vector = [vector; r];
   end   
end
Ec = mean(vector);


theta = (delta * (1 - normcdf(c_star, nu, sigma))*s  + normcdf(c_star, nu, sigma)*Ec ) / (1 - delta*(1-normcdf(c_star, nu, sigma)));

%%%  This is equation (8) in Appendix of Laibson.  
%%%   out should be equal to zero when c = C*

out = beta * delta * (s  + theta)  - c_star;