/*
** NPLLD.SRC -  Procedure that estimates a dynamic labor demand model with
**              lump-sum and linear hiring and firing costs, using a nested
**              fixed point algorithm.
**              Features:
**                  - Cobb-Douglas production function
**                  - AR(1) productivity shock
**                  - No fixed effects in productivity
**                  - Constant wages: over time and individuals
**                  - Adjustment costs: Lump-sum + linear + quadratic
**
** by Victor Aguirregabiria
**      October, 1997
**      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
**  ---------------------------------------------------------------------
**
** { tetaest , varest , proest } =
**            nplld(dpobs,dtobs,xpobs,xtobs,omobs,numdp,numdt,numome,fixed)
**
**  INPUTS:
**          dobs   - (nobs x 1) vector with the observations of
**                   the dependent variable (not discretized): ln(1+DL/L_1)
**
**          xobs   - (nobs x 1 ) vector with the observations of
**                   ln(employment) at period t-1 (not discretized)
**
**          omobs  - (nobs x 1 ) vector with the observations of
**                   productivity shock (not discretized)
**
**          numdom - (2 x 1) with the number of cells in the discretization
**                   of dobs and omobs: (numd | numom)
**
**          fixed  - (K x 1) vector with the binary indicator of whether
**                   a parameter is fixed (=0) or it will be estimated (=1).
**
**          q0     - (K x 1) vector with initial values for struc. parameters
**     tetavec = {alfa0,alfaL,wage,FH,FF,LH,LF,QH,QF,rho,sigmaa,beta,sigeps}
**
**  OUTPUTS:
**          tetaest - (K x 1) vector with parameter estimates
**
**          varest  - (K x K) covariance matrix of the estimates.
**
**          proest  - (numxom x numd) matrix of estimated choice probabilities.
**                    where numxom = numx * numom
**
**  ---------------------------------------------------------------------
**                      PROCEDURES CALLED
**  ---------------------------------------------------------------------
**          diskld.src  :   Discretization
**
**
**  ---------------------------------------------------------------------
*/

proc (3) = nplld(dpobs,dtobs,xpobs,xtobs,omobs,numdp,numdt,numome,fixed) ;
  local nameteta, myzero,
        alpha0, alphaL, lamda, sigmaa, rho, beta, wageper, wagetem,
        Fhireper, Ffireper, Lhireper, Lfireper, Qhireper, Qfireper,
        Fhiretem, Ffiretem, Lhiretem, Lfiretem, Qhiretem, Qfiretem,
        sigmaeps,
        dpvec, dtvec, xpvec, xtvec, omevec, tranome, tranxp, tranxt,
        numxp, numxt, inddobs, inddpobs, inddtobs, indxpobs, indxtobs,
        indomobs,
        explan, bp0 , vbp0,
        tetaest, varest, pest ;

  "" ;
  "/*************************/" ;
  "/* 0. Loading parameters */" ;
  "/*************************/" ;
  "" ;
  load alpha0   ; load alphaL   ; load lamda    ; load sigmaa   ;
  load rho      ; load beta     ; load wageper  ; load wagetem  ;
  load Fhireper ; load Ffireper ; load Lhireper ; load Lfireper ;
  load Qhireper ; load Qfireper ; load Fhiretem ; load Ffiretem ;
  load Lhiretem ; load Lfiretem ; load Qhiretem ; load Qfiretem ;
  load sigmaeps ;


  nameteta = ("alfa0"    | "alfaL"    | "lamda"    | "sigmaa"   |
              "rho"      | "beta"     | "wageper"  | "wagetem"  |
              "Fhireper" | "Ffireper" | "Fhiretem" | "Ffiretem" |
              "Lhireper" | "Lfireper" | "Lhiretem" | "Lfiretem" |
              "Qhireper" | "Qfireper" | "Qhiretem" | "Qfiretem" |
              "sigmaeps" ) ;

  myzero = 1e-12 ;

  /*****************************************************/
  /* 1. Discretization of state and decision variables */
  /*****************************************************/

  { dpvec, dtvec, xpvec, xtvec, omevec,
    inddpobs, inddtobs, indxpobs, indxtobs, indomobs,
    tranxp, tranxt, tranome }
          = disestld(dpobs,dtobs,xpobs,xtobs,omobs,numdp,numdt,numome) ;
  numxp = rows(xpvec) ;
  numxt = rows(xtvec) ;

  inddpobs = discthre(dpobs,dpvec[1:numdp-1]+((dpvec[2]-dpvec[1])/2)) ;
  inddtobs = discthre(dtobs,dtvec[1:numdt-1]+((dtvec[2]-dtvec[1])/2)) ;
  indxpobs = discthre(xpobs,xpvec[1:numxp-1]+((xpvec[2]-xpvec[1])/2)) ;
  indxtobs = discthre(xtobs,xtvec[1:numxt-1]+((xtvec[2]-xtvec[1])/2)) ;
  indomobs = discthre(omobs,omevec[1:numome-1]+((omevec[2]-omevec[1])/2)) ;


  title("Rate of change in permanent employment \L Discretized distribution") ;
  xlabel("Rate of change in permanent employment (%)") ;
  ylabel("Frequency") ;
  histp(dpvec[inddpobs],200) ;

  title("Logarithm of permanent employment \L Discretized distribution") ;
  xlabel("Logarithm of permanent employment") ;
  ylabel("Frequency") ;
  histp(xpvec[indxpobs],200) ;

  title("Rate of change in temporary employment \L Discretized distribution") ;
  xlabel("Rate of change in temporary employment (%)") ;
  ylabel("Frequency") ;
  histp(dtvec[inddtobs],200) ;

  title("Logarithm of 1 + temporary employment \L Discretized distribution") ;
  xlabel("Logarithm of 1 + temporary employment") ;
  ylabel("Frequency") ;
  histp(xtvec[indxtobs],200) ;

  title("Productivity shock \L Discretized distribution") ;
  xlabel("Productivity shock") ;
  ylabel("Frequency") ;
  histp(omevec[indomobs],200) ;



  /********************************************************/ ;
  /* 2. Initial Choice Probabilities (Estimates from MNL) */ ;
  /********************************************************/ ;
  inddobs = (inddpobs-1)*numdt + inddtobs ;

  explan = ones(rows(xpobs),1)~xpobs~xtobs~omobs ;
  { bp0 , vbp0 } = multilog(inddobs,explan) ;

  tetaest = 0 ;
  varest = 0 ;
  pest = 0 ;

  retp(tetaest,varest,pest) ;
endp ;











