/*
** DISSOLLD.SRC -   Procedure that discretizes decisions and state variables
**                  and obtains matrices of transition probabilities
**                  in a dynamic labor demand model with two types of labor
**
** 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 of the procedure:
**
** { dpvec, dtvec, xpvec, xtvec, omevec, tranxp, tranxt, tranome } =
**                  dissolld(numdp,numdt,numxp,numxt,numome)
**
**  INPUTS:
**          numdp   -  number of cells for dp (should be a odd number)
**
**          numdt   -  number of cells for dt (should be a odd number)
**
**          numxp   -  number of cells for xp
**
**          numxt   -  number of cells for xt
**
**          numome  -  number of cells for omega
**
**  OUTPUTS:
**          dpvec   - (numdp x 1) vector with discretized values of dp.
**
**          dtvec   - (numdt x 1) vector with discretized values of dt.
**
**          xpvec   - (numxp x 1) vector with discretized values of xp.
**
**          xtvec   - (numxt x 1) vector with discretized values of xt.
**
**          omevec  - (numome x 1) vector with discretized values of omega
**
**          tranxp  - (numdp x numxp) vector with the tran rule xp'= xp+dp
**
**          tranxt  - (numdt x numxt) vector with the tran rule xt'= xt+dt
**
**          tranome - (numome x numome) matrix of transition probabilities omega
**
*/

proc (8) = dissolld(numdp,numdt,numxp,numxt,numome) ;
  local sigmaa, rho, sigmaom, minome, maxome, stepome, omevec,
        mindp, maxdp, splitdp, stepdp, dpvec,
        mindt, maxdt, splitdt, stepdt, dtvec,
        minxp, maxxp, stepxp, xpvec,
        minxt, maxxt, stepxt, xtvec,
        omabove, ombelow, tranome, buffx, buffd, tranxp, tranxt ;

  load sigmaa ; load rho ; load alpha0 ; load alphaL ;
  load wageper ; load wagetem ; load lamda ;
  sigmaom = sigmaa / sqrt(1-rho*rho) ;

  @ 1. Discretization of omega @
    minome = -2*sigmaom ;
    minome = round(10000 * minome)/10000  ;
    maxome = 2*sigmaom ;
    maxome = round(10000 * maxome)/10000  ;
    stepome = round( 10000 * (maxome-minome)/(numome-1) )/10000  ;
    omevec  = seqa(minome,stepome,numome) ;
    omevec = round( 10000 * omevec ) / 10000 ;

  @ 2. Discretization of ln(Permanent Employment) @
    minxp = (1/(1-alphaL))*(2*ln(alphaL)+alpha0+omevec[1]-ln(wageper)) ;
    minxp = 0 ;
    minxp = round( 10000 * minxp )/10000 ;
    maxxp = (1/(1-alphaL))*(ln(alphaL)+alpha0+omevec[numome]-ln(wageper));
    maxxp = round( 10000 * maxxp )/10000 ;
    stepxp = round( 10000 * (maxxp-minxp)/(numxp-1) )/10000  ;
    xpvec = seqa(minxp,stepxp,numxp) ;
    xpvec = round( 10000 * xpvec )/10000 ;

  @ 3. Discretization of Temporary Employmnet @
    minxt = (1/(1-alphaL))*(2*ln(alphaL)+alpha0+omevec[1]-ln(wagetem))
          - ln(lamda) ;
    minxt = 0 ;
    minxt = round( 10000 * minxt )/10000 ;
    maxxt = (1/(1-alphaL))*(ln(alphaL)+alpha0+omevec[numome]-ln(wagetem))
          - ln(lamda) ;
    maxxt = round( 10000 * maxxt )/10000 ;
    stepxt = round( 10000 * (maxxt-minxt)/(numxt-1) )/10000  ;
    xtvec = seqa(minxt,stepxt,numxt) ;
    xtvec = round( 10000 * xtvec )/10000 ;

  @ 4. Discretization of decision dp @
    splitdp = (numdp-1)/2 + 1 ;
    if splitdp/=int(splitdp) ;
    "Error: numdp is not an odd integer" ;
       end ;
    endif ;
    mindp = -(splitdp-1)*stepxp ;
    maxdp = (splitdp-1)*stepxp ;
    stepdp = stepxp ;
    dpvec  = seqa(mindp,stepdp,numdp) ;
    dpvec = round( 10000 * dpvec )/10000 ;

  @ 5. Discretization of decision dt @
    splitdt = (numdt-1)/2 + 1 ;
    if splitdt/=int(splitdt) ;
    "Error: numdt is not an odd integer" ;
       end ;
    endif ;
    mindt = -(splitdt-1)*stepxt ;
    maxdt = (splitdt-1)*stepxt ;
    stepdt = stepxt ;
    dtvec  = seqa(mindt,stepdt,numdt) ;
    dtvec = round( 10000 * dtvec )/10000 ;

  @ 6. Transition probability matrix of omega @
    omabove = omevec + stepome/2 ;
    ombelow = omevec - stepome/2 ;
    ombelow[1] = -1e12 ;
    omabove[numome] = 1e12 ;
    tranome = cdfn( (omabove - rho*(omevec'))/sigmaa )
            - cdfn( (ombelow - rho*(omevec'))/sigmaa ) ;
    tranome = tranome./sumc(tranome) ;

  @ 7. Transition rule xp' = xp + dp @
    buffx = seqa(1,1,numxp) ;
    buffd = seqa(1,1,numdp) - splitdp ;
    tranxp = (buffx') + buffd ;
    tranxp = (tranxp.>=1).*(tranxp.<=numxp).*tranxp
          + (tranxp.<1).*1
          + (tranxp.>numxp).*numxp ;

  @ 8. Transition rule xt' = xt + dt @
    buffx = seqa(1,1,numxt) ;
    buffd = seqa(1,1,numdt) - splitdt ;
    tranxt = (buffx') + buffd ;
    tranxt = (tranxt.>=1).*(tranxt.<=numxt).*tranxt
          + (tranxt.<1).*1
          + (tranxt.>numxt).*numxt ;

  retp(dpvec,dtvec,xpvec,xtvec,omevec,tranxp,tranxt,tranome) ;
endp ;
