/*
** SOLVLD.SRC -     Procedure that solves a labor demand model with
**                  with two types of labor
**                  The algorithm exploits a "smooth" Bellman's equation.
**                  The procedure returns a vector containing the
**                  "smooth" value function, and a vector with the
**                  conditional choice probabilities
**
** by Victor Aguirregabiria
**      January, 1999
**      Revised: 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:
**
** { svfun , procho } =
**          solvld(dpvec,dtvec,xpvec,xtvec,omevec,tranxp,tranxt,tranome)
**
**  INPUTS:
**          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 of transitions: xp' = xp+dp
**
**          tranxt - (numdt x numxt) vector of transitions: xt' = xt+dt
**
**          tranome- (numome x numome) matrix of transition probs. of omega
**
**
**  OUTPUTS:
**          svfun   - (numxp * numxt * numphi x 1) vector smooth value function
**
**          procho  - (numxp * numxt * numphi x numdp * numdt) matrix with the
**                    optimal choice probabilities.
**
*/

proc (2) = solvld(dpvec,dtvec,xpvec,xtvec,omevec,tranxp,tranxt,tranome);
  local alpha0, alphaL, lamda, wageper, wagetem, sigmaa, rho, beta, sigmaeps,
        Fhireper, Ffireper, Fhiretem, Ffiretem,
        Lhireper, Lfireper, Lhiretem, Lfiretem,
        Qhireper, Qfireper, Qhiretem, Qfiretem,
        sigmaphi, numdp, numdt, numxp, numxt, numome, numstate,
        smat, dmat, xplusdp, xplusdt, profit,
        ans1, ansome, ansxp, ansxt, myopic, vomexpxt,
        tranxpxt, V0, epsconv, criter, iter, V1, procho,
        grav0, grapcho ;

  "" ;
  "/*************************/" ;
  "/* 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 ;

  sigmaphi = sigmaa / sqrt(1-rho*rho) ;
  numdp = rows(dpvec) ;
  numdt = rows(dtvec) ;
  numxp = rows(xpvec) ;
  numxt = rows(xtvec) ;
  numome = rows(omevec) ;
  numstate = numxp * numxt * numome ;

  "" ;
  "/************************/" ;
  "/* 1. Matrix of profits */" ;
  "/************************/" ;
  "" ;
  smat = zeros(3,numstate) ;
  smat[1,.] = (xpvec') .*. ones(1,numxt*numome) ;
  smat[2,.] = ones(1,numxp) .*. ((xtvec') .*.  ones(1,numome)) ;
  smat[3,.] = ones(1,numxp*numxt) .*. (omevec') ;

  dmat = zeros(numdt*numdp,2) ;
  dmat[.,1] = dpvec .*. ones(numdt,1) ;
  dmat[.,2] = ones(numdp,1) .*. dtvec ;

  xplusdp = smat[1,.] + dmat[.,1] ;
  xplusdt = smat[2,.] + dmat[.,2] ;

  @ Profit function @

  profit = exp(alpha0 + smat[3,.] + alphaL.*xplusdp
               + alphaL*ln(1+lamda*(exp(xplusdt)-1)./exp(xplusdp)) )
         - wageper * exp(xplusdp)
         - wagetem * (exp(xplusdt)-1)
         - Fhireper * (dmat[.,1].>0)
         - Ffireper * (dmat[.,1].<0)
         - Fhiretem * (dmat[.,2].>0)
         - Ffiretem * (dmat[.,2].<0)
         - Lhireper * (dmat[.,1].>0) .* (exp(xplusdp) - exp(smat[1,.]))
         + Lfireper * (dmat[.,1].<0) .* (exp(xplusdp) - exp(smat[1,.]))
         - Lhiretem * (dmat[.,2].>0) .* (exp(xplusdt) - exp(smat[2,.]))
         + Lfiretem * (dmat[.,2].<0) .* (exp(xplusdt) - exp(smat[2,.]))
         - Qhireper * (dmat[.,1].>0) .* (exp(dmat[.,1])-1).*(exp(dmat[.,1])-1)
         - Qfireper * (dmat[.,1].<0) .* (exp(dmat[.,1])-1).*(exp(dmat[.,1])-1)
         - Qhiretem * (dmat[.,2].>0) .* (exp(dmat[.,2])-1).*(exp(dmat[.,2])-1)
         - Qfiretem * (dmat[.,2].<0) .* (exp(dmat[.,2])-1).*(exp(dmat[.,2])-1);

  "" ;
  "/*******************************/" ;
  "/* 2. Graph of profit function */" ;
  "/*******************************/" ;
  "" ;
  ans1 = "y" ;
  do while (ans1$=="Y")OR(ans1$=="y") ;
    ans1 = "" ;
    do while (ans1$/="Y")AND(ans1$/="y")
          AND(ans1$/="N")AND(ans1$/="n") ;
      "Do you want to graph the profit function ? (Answer Y or N)" ;
      ans1 = cons ;
    endo ;
    if (ans1$=="Y")OR(ans1$=="y") ;
      vomexpxt = int(numdp/2)*numdt + int(numdt/2) + 1 ;
      myopic = reshape(profit[vomexpxt,.],numxp*numxt,numome) ;
      myopic = reshape(myopic[.,int(numome/2)+1],numxp,numxt) ;
      title("Profit as function of xp");
      xlabel("xp");
      ylabel("Profit");
      xy(xpvec,myopic) ;
      title("Profit as function of xt");
      xlabel("xt");
      ylabel("Profit");
      xy(xtvec,myopic') ;
    endif ;
  endo ;

  "" ;
  "/*******************************/" ;
  "/* 3. Graph of myopic decision */" ;
  "/*******************************/" ;
  "" ;
  ans1 = "" ;
  do while (ans1$/="Y")AND(ans1$/="y")
        AND(ans1$/="N")AND(ans1$/="n") ;
    "Do you want to graph the choice probabilities of" ;
    "decision variables for the myopic model ? (Y/N)" ;
    ans1 = cons ;
  endo ;
  if (ans1$=="Y")OR(ans1$=="y") ;
    myopic = maxc(profit) ;
    myopic = ( profit - (myopic') )/sigmaeps ;
    myopic = exp(myopic)./(sumc(exp(myopic))') ;
    ansome = 0 ;
    do while (ansome<1)OR(ansome>numome) ;
      "Cell for omega:";; "between 1 and ";; numome ;;
      ansome = con(1,1) ;
    endo ;
    ansxp = 0 ;
    do while (ansxp<1)OR(ansxp>numxp) ;
      "Cell for xp:";; "between 1 and ";; numxp ;
      ansxp = con(1,1) ;
    endo ;
    ansxt = 0 ;
    do while (ansxt<1)OR(ansxt>numxt) ;
      "Cell for xt:";; "between 1 and ";; numxt ;
      ansxt = con(1,1) ;
    endo ;
    vomexpxt = (ansxp-1)*numxt*numome
             + (ansxt-1)*numome + ansome ;
    myopic = reshape(myopic[.,vomexpxt],numdp,numdt) ;

    title("Joint distribution of dp and dt");
    xlabel("dt");
    ylabel("dp");
    zlabel("Pr( dp,dt | xp,xt,omega )");
    surface(dtvec',dpvec,myopic) ;
  endif ;
  clear myopic ;

  /****************************/
  /* 3. Fixed point algorithm */
  /****************************/
  tranxp = tranxp .*. ones(numdt,numxt) ;
  tranxp = reshape(tranxp,numxp*numxt*numdp*numdt,1) ;
  tranxt = ones(numdp,numxp) .*. tranxt ;
  tranxt = reshape(tranxt,numxp*numxt*numdp*numdt,1) ;
  tranxpxt = (tranxp-1)*numxt+tranxt ;

  V0 = zeros(numxp*numxt*numome,1) ;
  epsconv = 1e-03 ;
  criter = 1e+06 ;

  iter=1 ;
  do while criter>epsconv ;
    "Iteration = " iter ;
    "Criterion = " criter ;
    V1 = reshape(V0,numxp*numxt,numome) ;
    V1 = V1 * tranome ;
    V1 = V1[tranxpxt,.] ;
    V1 = reshape(V1,numdp*numdt,numstate) ;
    V1 = (profit + beta * V1)./sigmaeps ;
    V1 = sigmaeps .* ( ln( sumc( exp(V1 - (maxc(V1)')) ) ) + maxc(V1) ) ;

    criter = maxc( abs(V1-V0) ) ;
    V0 = V1 ;
    iter=iter +1 ;
  endo ;

  V1 = reshape(V0,numxp*numxt,numome) ;
  V1 = V1 * tranome ;
  V1 = V1[tranxpxt,.] ;
  V1 = reshape(V1,numdp*numdt,numstate) ;
  V1 = (profit + beta * V1)./sigmaeps ;

  procho = maxc(V1) ;
  procho = V1 - (procho') ;
  procho = exp(procho)./(sumc(exp(procho))') ;

  /********************************/ ;
  /* 4. Graph of dynamic decision */ ;
  /********************************/ ;
  ans1 = "" ;
  do while (ans1$/="Y")AND(ans1$/="y")
        AND(ans1$/="N")AND(ans1$/="n") ;
    "Do you want to graph the value function ? (Y/N)" ;
    ans1 = cons ;
  endo ;
  if (ans1$=="Y")OR(ans1$=="y") ;
    ansome = 0 ;
    do while (ansome<1)OR(ansome>numome) ;
      "Cell for omega:";; "between 1 and ";; numome ;
      ansome = con(1,1) ;
    endo ;
    grav0 = reshape(V0,numxp*numxt,numome) ;
    grav0 = reshape(grav0[.,ansome],numxp,numxt) ;
    title("Smooth value function");
    xlabel("xt");
    xlabel("xp");
    zlabel("V( xp , xt, omega )");
    surface(xtvec',xpvec,grav0) ;
  endif ;

  ans1 = "" ;
  do while (ans1$/="Y")AND(ans1$/="y")
        AND(ans1$/="N")AND(ans1$/="n") ;
    "Do you want to graph the choice probabilities ? (Y/N)" ;
    ans1 = cons ;
  endo ;
  if (ans1$=="Y")OR(ans1$=="y") ;
    ansome = 0 ;
    do while (ansome<1)OR(ansome>numome) ;
      "Cell for omega:";; "between 1 and ";; numome ;
      ansome = con(1,1) ;
    endo ;
    ansxp = 0 ;
    do while (ansxp<1)OR(ansxp>numxp) ;
      "Cell for xp:";; "between 1 and ";; numxp ;;
      ansxp = con(1,1) ;
    endo ;
    ansxt = 0 ;
    do while (ansxt<1)OR(ansxt>numxt) ;
      "Cell for xt:";; "between 1 and ";; numxt ;;
      ansxt = con(1,1) ;
    endo ;
    vomexpxt = (ansxp-1)*numxt*numome
             + (ansxt-1)*numome + ansome ;
    grapcho = reshape(procho[.,vomexpxt],numdp,numdt) ;
    title("Joint distribution of dp and dt");
    xlabel("dt");
    xlabel("dp");
    zlabel("Pr( dp,dt | xp,xt,omega )");
    surface(dtvec',dpvec,grapcho) ;
  endif ;

  retp(V0,procho) ;
endp ;




