/*
** SIMLD.SRC -     Procedure that simulates data of state and decision
**                  variables in a labour demand model with
**                  linear hiring and firing costs.
**
** by Victor Aguirregabiria
** April, 1998
** Last revision: April, 2001
**
** For comments send e-mail to: vaguirre@bu.edu
**
**  ---------------------------------------------------------------------
**                              MODEL
**  ---------------------------------------------------------------------
**      (1) Variables:
**      --------------
**              xp = ln( permanent employment at t-1 )
**              xt = ln( 1 + temporary employment at t-1 )
**              dp = ln(per. employment at t) - ln(per. employment at t-1)
**              dt = ln(1 + tem. employment at t) - ln(1 + tem. emp at t-1)
**              omega = Productivity shock
**              State variables = { xp , xt , omega }
**              Decision variables = { dp , dt }
**
**      (2) One-period profit function:
**      -------------------------------
**          R(x,d,omega,eps) =
**  (Output)           y(xp+dp,xt+dt,omega)
**  (Wage bill)      - wageper * exp(xp+dp)
**                   - wagetem * exp(xt+dt)
**  (Fixed ACs)      - Fhireper * I(dp>0)
**                   - Ffireper * I(dp<0)
**                   - Fhiretem * I(dt>0)
**                   - Ffiretem * I(dt<0)
**  (Lin. ACs)       - Lhireper * I(dp>0) * abs(exp(xp+dp) - exp(xp))
**                   + Lfireper * I(dp<0) * abs(exp(xp+dp) - exp(xp))
**                   - Lhiretem * I(dt>0) * abs(exp(xt+dt) - exp(xt))
**                   + Lfiretem * I(dt<0) * abs(exp(xt+dt) - exp(xt))
**  (Quad. ACs)      - Qhireper * I(dp>0) * (exp(dp) - 1)*(exp(dp)-1)
**                   - Qfireper * I(dp<0) * (exp(dp) - 1)*(exp(dp)-1)
**                   - Qhiretem * I(dt>0) * (exp(dt) - 1)*(exp(dt)-1)
**                   - Qfiretem * I(dt<0) * (exp(dt) - 1)*(exp(dt)-1)
**  (Unobservable)   + eps(dp,dt)
**
**      Production Function
**      -------------------
**      y(xp+dp,xt+xt,omega) =  exp( alpha0 + omega + alphaL*(xp + dp)
**                           +       alphaL*lamda*exp(xt+dt)/exp(xp+dp) )
**
**      (3) Conditional transition probabilities of observable state variable :
**      -----------------------------------------------------------------------
**              xp[t+1] = xp[t] + dp[t]
**              xt[t+1] = xt[t] + dt[t]
**              omega[t+1] = rho * omega[t] + a[t], where a[t] iid N(0,sigmaa)
**
**      (4) Distribution of unobservables:
**      ----------------------------------
**              eps(dp,dt) is iid over (i,t,dp,dt) Extreme Value(sigmaeps)
**
**      (5) Parameters of the model:
**      ----------------------------
**      tetavec = { alpha0, alphaL, lamda, wageper, wagetem, rho, sigmaa, beta,
**                  Fhireper, Ffireper, Fhiretem, Ffiretem,
**                  Lhireper, Lfireper, Lhiretem, Lfiretem,
**                  Qhireper, Qfireper, Qhiretem, Qfiretem,
**                  sigmaeps }
**      where: beta = discount factor
**
**  ----------------------------------------------------------------------------
**
**  Format:
** { dpsim, dtsim, xpsim, xtsim, omesim, vfsim } =
** simld(probs,valfun,dpvec,dtvec,xpvec,xtvec,omevec,tranome,numN,numT,jump)
**
**  INPUTS:
**
**  probs  - (numdp*numdt x numxp*numxt*numome) matrix with indexes
**           optimal choice probabilities.
**
**  valfun - (numxp*numxt*numome x 1) vector with the value function
**
**  dpvec  - (numdp x 1) vector with the values of dp
**
**  dtvec  - (numdt x 1) vector with the values of dt
**
**  xpvec  - (numxp x 1) vector with the values of xp
**
**  xtvec  - (numxt x 1) vector with the values of xt
**
**  omevec - (numome x 1) vector with the values of omega
**
**  tranome- (numome x numome) matrix of transition probs of omega
**
**  numN   - Number of individuals in the simulations
**
**  numT   - Number of time periods in the simulations
**
**  jump   - Number of time periods which are used to obtain initial
**           simulations but which are not returned as simulations
**
**  OUTPUTS:
**
**  dpsim  - (numN x numT) with the simulations for dp
**
**  dtsim  - (numN x numT) with the simulations for dt
**
**  xpsim  - (numN x numT) with the simulations for xp
**
**  xtsim  - (numN x numT) with the simulations for xt
**
**  omesim - (numN x numT) with the simulations for omega
**
**  vfsim -  (numN x numT) with the simulations for value function
*/

proc (6) =
simld(probs,valfun,dpvec,dtvec,xpvec,xtvec,omevec,tranome,numN,numT,jump);
  local seed, numdp, numdt, numxp, numxt, numome, numstate,
        splitdp, splitdt, Pome_hi, Pome_lo, probs0, probs1,
        omesim, dpsim, dtsim, xpsim, xtsim, statesim,
        pbuff0, pbuff1, psid, t, psiome, vfsim ;


  seed = 7387905 ;
  numT = numT + jump ;
  numdp = rows(dpvec) ;
  numdt = rows(dtvec) ;
  numxp = rows(xpvec) ;
  numxt = rows(xtvec) ;
  numome = rows(omevec) ;
  numstate = numxp*numxt*numome ;
  splitdp = int(numdp/2) + 1 ;
  splitdt = int(numdt/2) + 1 ;
  Pome_hi = cumsumc(tranome) ;
  Pome_lo = ( zeros(1,numome) | Pome_hi[1:numome-1,.] ) ;

  probs = cumsumc(probs) ;
  probs1 = probs ;
  clear probs ;
  probs0 = zeros(1,numstate) | probs1[1:numdp*numdt-1,.] ;

  @ Simulations @
  omesim = zeros(numN,numT) ;
  dpsim  = zeros(numN,numT) ;
  dtsim  = zeros(numN,numT) ;
  xpsim  = zeros(numN,numT) ;
  xtsim  = zeros(numN,numT) ;

  @ Iteration 1 @
  xpsim[.,1]  = int(numxp*rndus(numN,1,seed)) + 1 ;
  xtsim[.,1]  = int(numxt*rndus(numN,1,seed)) + 1 ;
  omesim[.,1] = int(numome*rndus(numN,1,seed)) + 1 ;
  statesim = (xpsim[.,1]-1)*numxt*numome
           + (xtsim[.,1]-1)*numome + omesim[.,1] ;
  pbuff0 = probs0[.,statesim]';
  pbuff1 = probs1[.,statesim]';
  psid = rndus(numN,1,seed) ;
  psid = maxindc( ((psid.>pbuff0).*(psid.<=pbuff1))' ) ;
  dpsim[.,1] = int((psid-1)./numdt) + 1 ;
  dtsim[.,1] = psid - (dpsim[.,1]-1).*numdt ;
  xpsim[.,2] = xpsim[.,1] + dpsim[.,1] - splitdp ;
  xpsim[.,2] = (xpsim[.,2].>=1).*(xpsim[.,2].<=numxp).*xpsim[.,2]
             + (xpsim[.,2].<1).*1
             + (xpsim[.,2].>numxp).*numxp ;
  xtsim[.,2] = xtsim[.,1] + dtsim[.,1] - splitdt ;
  xtsim[.,2] = (xtsim[.,2].>=1).*(xtsim[.,2].<=numxt).*xtsim[.,2]
             + (xtsim[.,2].<1).*1
             + (xtsim[.,2].>numxt).*numxt ;
  t=2;
  do while t<=(numT-1) ;
    psiome = rndus(numN,1,seed) ;
    omesim[.,t] = maxindc( ((psiome') .> Pome_lo[.,omesim[.,t-1]])
                      .*   ((psiome') .<=Pome_hi[.,omesim[.,t-1]]) );
    statesim = (xpsim[.,t]-1)*numxt*numome
             + (xtsim[.,t]-1)*numome + omesim[.,t] ;
    pbuff0 = probs0[.,statesim]';
    pbuff1 = probs1[.,statesim]';
    psid = rndus(numN,1,seed) ;
    psid = maxindc( ((psid.>pbuff0).*(psid.<=pbuff1))' ) ;
    dpsim[.,t] = int((psid-1)./numdt) + 1 ;
    dtsim[.,t] = psid - (dpsim[.,t]-1).*numdt ;
    xpsim[.,t+1] = xpsim[.,t] + dpsim[.,t] - splitdp ;
    xpsim[.,t+1] = (xpsim[.,t+1].>=1).*(xpsim[.,t+1].<=numxp)
                   .*xpsim[.,t+1]
                 + (xpsim[.,t+1].<1).*1
                 + (xpsim[.,t+1].>numxp).*numxp ;
    xtsim[.,t+1] = xtsim[.,t] + dtsim[.,t] - splitdt ;
    xtsim[.,t+1] = (xtsim[.,t+1].>=1).*(xtsim[.,t+1].<=numxt)
                   .*xtsim[.,t+1]
               + (xtsim[.,t+1].<1).*1
               + (xtsim[.,t+1].>numxt).*numxt ;
    t=t+1;
  endo ;

  @ Iteration numT @
  psiome = rndus(numN,1,seed) ;
  omesim[.,t] = maxindc( ((psiome') .> Pome_lo[.,omesim[.,t-1]])
                    .*   ((psiome') .<=Pome_hi[.,omesim[.,t-1]]) );
  statesim = (xpsim[.,t]-1)*numxt*numome
           + (xtsim[.,t]-1)*numome + omesim[.,t] ;
  pbuff0 = probs0[.,statesim]';
  pbuff1 = probs1[.,statesim]';
  psid = rndus(numN,1,seed) ;
  psid = maxindc( ((psid.>pbuff0).*(psid.<=pbuff1))' ) ;
  dpsim[.,t] = int((psid-1)./numdt) + 1 ;
  dtsim[.,t] = psid - (dpsim[.,t]-1).*numdt ;

  vfsim = (xpsim-1)*numxt*numome + (xtsim-1)*numome + omesim ;
  vfsim = reshape(valfun[vfsim],numN,numT) ;
  vfsim = vfsim[.,jump+1:numT] ;
  dpsim = reshape(dpvec[dpsim],numN,numT) ;
  dpsim = dpsim[.,jump+1:numT] ;
  dtsim = reshape(dtvec[dtsim],numN,numT) ;
  dtsim = dtsim[.,jump+1:numT] ;
  xpsim = reshape(xpvec[xpsim],numN,numT) ;
  xpsim = xpsim[.,jump+1:numT] ;
  xtsim = reshape(xtvec[xtsim],numN,numT) ;
  xtsim = xtsim[.,jump+1:numT] ;
  omesim = reshape(omevec[omesim],numN,numT) ;
  omesim = omesim[.,jump+1:numT] ;

  retp(dpsim,dtsim,xpsim,xtsim,omesim,vfsim) ;
endp ;


