% uf.m: utility function for retirement (employment) problem

  function [ufv,dufv]=uf(age,income,pension,work_state,work_decision)

  global work_disutility_age search_cost_age crra_coef u_scale unemp_prob;

  dufv=zeros(3,1);

  if (work_state > 0);  % a worker

   if (work_decision > 0);   % a worker decides to continue working

      a2=age*age;
      li=log(income);
      lp=log(pension);
      ufv=(1-unemp_prob)*(u_scale*exp(crra_coef*li)-work_disutility_age*a2)+...
           unemp_prob*u_scale*exp(crra_coef*lp);

      dufv(1)=-a2*(1-unemp_prob);
      dufv(2)=0.0;
      dufv(3)=(1-unemp_prob)*u_scale*li*exp(crra_coef*li)+unemp_prob*u_scale*lp*exp(crra_coef*lp);

   else;                     % a worker decides to stop working/retire, no disutility of work any more
    if (age > 55);
      lp=log(pension);
      ufv=u_scale*exp(crra_coef*lp);
      dufv(1)=0;
      dufv(2)=0;
      dufv(3)=lp*ufv;
     
    else;
      ufv=0;
      dufv(1)=0;
      durv(2)=0;
      dufv(3)=0;
    end;

   end;

  else; % work_state is 0 (not working at start of period)

   if (work_decision > 0);    % a retiree decides to go back to work and incurs a fixed
                              % "search cost" for doing so

     if (age > 55);
       a2=age*age;
       li=log(income);
       lp=log(pension);
       ufv=(1-unemp_prob)*(u_scale*exp(crra_coef*log(income+pension))-...
              (search_cost_age+work_disutility_age)*a2)+...
          unemp_prob*(u_scale*exp(crra_coef*log(pension))-search_cost_age*a2);
       dufv(1)=-a2*(1-unemp_prob);
       dufv(2)=-a2;
       dufv(3)=(1-unemp_prob)*u_scale*li*exp(crra_coef*li)+unemp_prob*u_scale*lp*exp(crra_coef*lp);
     else;
       a2=age*age; 
       li=log(income);
       lp=log(pension);
       ufv=(1-unemp_prob)*(u_scale*exp(crra_coef*log(income))-...
              (search_cost_age+work_disutility_age)*a2)+...
          unemp_prob*(u_scale*exp(crra_coef*log(pension))-search_cost_age*a2);
       dufv(1)=-a2*(1-unemp_prob);
       dufv(2)=-a2;
       dufv(3)=(1-unemp_prob)*u_scale*li*exp(crra_coef*li)+unemp_prob*u_scale*lp*exp(crra_coef*lp);
     end;

   else;                      % a retiree decides to stay retired/not working

    if (age > 55);
       lp=log(pension);
       ufv=u_scale*exp(crra_coef*lp);
       dufv(1)=0;
       dufv(2)=0;
       dufv(3)=lp*ufv;
    else;
       ufv=0;
       dufv(1)=0;
       dufv(2)=0;
       dufv(3)=0;
    end;

   end;

  end;
