% simulate_data.m: program to simulate retirement/employment choices for a sample of individuals
%                 John Rust, Georgetown University, November 2012


  setup;

  firstage=20;
  lastage=80;

  income_const=2000;
  income_linear_term=60;
  income_quadratic_term=-0.7;
  income_lag=.7;
  income_std=100;
  
  retirement_age=55;                 % first age at which a worker can collect retirement benefits
  pension_replacement_rate=.6;       % pension replacement rate if over retirement age
  unemployment_replacement_rate=.4;  % unemployment benefit replacement rate

  Npeople=500;

  data=[];

  for i=1:Npeople;

   income0=income_const+income_std*randn;

   work_state=1;

   for t=firstage:lastage;

      income=income_linear_term*t+income_quadratic_term*(t^2)+income_lag*income0+income_std*randn; 

      if (t == firstage);
         avwage=income;
      else;
         if (working);
           avwage=((t-firstage)/(t+1-firstage))*avwage+income/(t+1-firstage);
         end;
      end;

      if (t < retirement_age);
        pension=unemployment_replacement_rate*avwage;
      else;
        pension=pension_replacement_rate*avwage;
      end;

      uw=uf(t,income,pension,work_state,1);
      unw=uf(t,income,pension,work_state,0);

      if (lsigma > 0);

        mu=max(uw,unw);
        euw=exp((uw-mu)/lsigma);
        eunw=exp((unw-mu)/lsigma);

        prob_work=euw/(euw+eunw);

        if (rand <= prob_work);
 
           working=1;

        else

           working=0;

        end;

      else;

       if (uw > unw);
         working=1;
       else;
         working=0;
       end;

      end;
   
      laid_off=0;

      if (working > 0);
        if (rand <=  unemp_prob);
         laid_off=1;
        end;
      end; 

      data=[data; [i working work_state t income pension laid_off]];

%      fprintf('age: %i income: %g income0 %g avwage=%g pension=%g uw=%g unw=%g work decision=%i work state=%i laid_off=%i\n',...
%               t,income,income0,avwage,pension,uw,unw,working,work_state,laid_off);

      if (laid_off);
        work_state=0;
      else;
        work_state=working;
      end;
      income0=income;

   end;
 
   casei=find(data(:,1) == i);
   retage=find(data(casei,2) == 0);
   fprintf('consumer %i retired at %i number retirement events %i\n',i,data(casei(retage(1)),4),size(retage,1));
  
  end;

  save('data','data');
