function out = fun_exp(c_star,s, nu, sigma);

delta = 0.9998;

%  Inputs: c*, s, nu, sigma
%  c* = s + theta
%  This function's output is: s - c* + theta;
%  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*

vector = [];
for i=1:10000
   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)));

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