/*
** DISCKLD.SRC -    Procedure that discretizes decision and state variables
**                  and obtains matrices of transition probabilities
**                  in a labor demand model with linear hiring and
**                  firing costs
**
** by Victor Aguirregabiria
**      April, 1998
**      Revised: March, 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:
**
** { dvec, xvec, omvec, tranom, tranx } =
**      = disckld(dobs,xobs,omobs,numd,numom,sigmaa,rho)
**
**  INPUTS:
**          dobs    -  (nobs x 1) vector with the observations of d
**
**          xobs    -  (nobs x 1) vector with the observations of x
**
**          omobs   -  (nobs x 1) vector with the observations of omega
**
**          numd    -  number of cells for d (should be a odd number)
**
**          numom   -  number of cells for omega
**
**          sigmaa  -  Std. dev. iid shock in stochastic process of omega
**
**          rho     -  Autoregressive parameter in stochastic process of omega
**
**  OUTPUTS:
**          dvec    - (numd x 1) vector with discretized values of d.
**
**          xvec    - (numx x 1) vector with discretized values of x.
**
**          omvec   - (numom x 1) vector with discretized values of omega
**
**          tranom  - (numom x numom) matrix of transition probabilities omega
**
**          tranx - (numx * numd x 1) vector with the transition rule x'= x + d.
**
*/

proc (5) = disckld(dobs,xobs,omobs,numd,numom,sigmaa,rho) ;
  local sigmaom, minom, maxom, stepom, omvec,
        mind, maxd, splitd, stepd, dvec,
        minx, maxx, stepx, numx, xvec,
        omabove, ombelow, tranom, buffx, buffd, tranx ;

  sigmaom = sigmaa / sqrt(1-rho*rho) ;

  "/* 1. Discretization of phi */" ;
    minom = pctiles(omobs,2) ;
    minom = round(10000 * minom)/10000  ;
    maxom = pctiles(omobs,98) ;
    maxom = round(10000 * maxom)/10000  ;
    stepom = round( 10000 * (maxom-minom)/(numom-1) )/10000  ;
    omvec  = seqa(minom,stepom,numom) ;
    omvec = round( 10000 * omvec ) / 10000 ;

  "/* 2. Discretization of decision */" ;
    maxd = minc ( abs(pctiles(dobs,2)) | abs(pctiles(dobs,98)) ) ;
    maxd = round(10000 * maxd)/10000  ;
    mind = -maxd ;
    splitd = (numd-1)/2 + 1 ;
    if splitd/=int(splitd) ;
    "Error: numd is not an odd integer" ;
       end ;
    endif ;
    stepd = round( 10000 * (maxd-mind)/(numd-1) )/10000  ;
    dvec  = seqa(mind,stepd,numd) ;
    dvec = round( 10000 * dvec )/10000 ;

  "/* 3. Discretization of Permanent Employmnet */" ;
    minx = pctiles(xobs,2) ;
    minx = round(10000 * minx)/10000  ;
    maxx = pctiles(xobs,98) ;
    maxx = round(10000 * maxx)/10000  ;
    stepx = stepd ;
    numx = int( (maxx-minx)/stepx ) + 1 ;
    xvec = seqa(minx,stepx,numx) ;
    xvec = round( 10000 * xvec )/10000 ;

  "/* 4. Transition probability matrix of omega */" ;
    omabove = omvec + stepom/2 ;
    ombelow = omvec - stepom/2 ;
    ombelow[1] = -1e12 ;
    omabove[numom] = 1e12 ;
    tranom = cdfn( (omabove - rho*(omvec'))/sigmaa )
            - cdfn( (ombelow - rho*(omvec'))/sigmaa ) ;
    tranom = tranom./sumc(tranom) ;

  "/* 5. Transition rule x' = x + d */" ;
    buffx = seqa(1,1,numx) ;
    buffd = seqa(1,1,numd) - splitd ;
    tranx = (buffx') + buffd ;
    tranx = (tranx.>=1).*(tranx.<=numx).*tranx
          + (tranx.<1).*1
          + (tranx.>numx).*numx ;
    tranx = reshape(tranx,numx*numd,1) ;

  retp(dvec, xvec, omvec, tranom, tranx) ;
endp ;


