/*
**  NPLDYGAM.SRC   Procedure that estimates the structural parameters 
**                 of dynamic game of firms' entry/exit
**                 using a Nested Pseudo-Likelihood (NPL) algorithm
**
**  by Victor Aguirregabiria
**  First version: February, 2003
**  Last revision: December, 2003
**
**  ---------------------------------------------------------------------
**                              MODEL
**  ---------------------------------------------------------------------
**  i = Firm index; t = Time index
**      a[it] = Indicator of the event 'firm i operates 
**               in the market at period t'
**      z[t] = Characteristics of the market at period t.
**      e0[it] and e1[it] = Firm i's private information.
**
**  Profits:
**  --------
**  No operative: U(0) = e0[it]
**
**  Operative:    U(1) = alphai - thetaE*(1-a[i,t-1]) 
**                     + beta * z[t]
**                     - delta * log(N[-it]) + 1) + e1[it]
**  where:
**      alphai = Firm i's fixed effect
**      N[-it] = Number of firms, other than i, 
**               operating in the market at period t. 
**      thetaE, beta and delta are parameters
**
**  ----------------------------------------------------------------------------
**
**
**  FORMAT:
** {best,varb} = 
**      npldygam(aobs,zobs,aobs_1,zval,ptranz,pchoice,bdisc,kiter)
**
**  Inputs:
**      aobs    - (nobs x nplayer) matrix with observations of
**                firms' activity decisions (1=active ; 0=no active)
**      zobs    - (nobs x 1) vector with observations of market
**                exogenous characteristics.
**      aobs_1  - (nobs x nplayer) matrix with observations of
**                firms' initial states (1=incumbent; 0=potential entrant)
**      zval    - (numz x 1) vector with values of market characteristics
**      ptranz  - (numz x numz) matrix of transition probabilities
**                of market characteristics.
**      pchoice - ((numz*2^nplayer) x nplayer) matrix of players'
**                choice probabilities used to initialize the procedure.
**      bdisc   - Discount factor
**      kiter   - Number of NPL iterations
**
**  Outputs:
**      best    - (kparam x kiter) matrix with estimates of parameters
**                (alpha_1 | alpha_2 | ...    | alpha_n | 
**                    beta | delta   | thetaE )
**                First column is the 1st-stage NPL, second column is
**                the 2nd-stage NPL, and so on.
**      varb    - (kparam x kaparam*kiter) matrix with estimated 
**                covariance matrices of the NPL estimators. 
**
*/

proc (2) = npldygam(aobs,zobs,aobs_1,zval,ptranz,pchoice,bdisc,kiter);
  local eulerc, myzero, nobs, nplayer, numa, numz, numx, 
        kparam, best, varb, mstate, twop, indzobs, indobs,
        i, namesb, aval, iter, ptrana, invi_bf, j, buff, 
        uobs0, uobs1, eobs0, eobs1, mi, ppi, ptranai0, 
        ptranai1, hi, umat0, umat1, sumu, sume, ww, 
        utilda0, utilda1, etilda0, etilda1, tetaest, 
        varest, u0, u1, e0, e1 ;        
        
  @ ---------------@
  @ Some constants @
  @ ---------------@                
  eulerc = 0.5772 ;
  myzero = 1e-16 ;
  nobs = rows(aobs) ;
  nplayer = cols(aobs) ;
  numa = 2^nplayer ;
  numz = rows(zval) ;  
  numx = numz*numa ;
  kparam = nplayer + 3 ;
  best = zeros(kparam,kiter) ;
  varb = zeros(kparam,kparam*kiter) ;  
  namesb = (0 $+ "alpha" $+ ftocv(seqa(1,1,nplayer),1,0))
         | "betaz" | "delta" | "thetaE" ; 
   
  @ ----------------------------------------------@
  @ Matrix with values of states of (z[t],a[t-1]) @
  @ ----------------------------------------------@
  aval = zeros(numa,nplayer) ; 
  i=1 ;
  do while i<=nplayer ;
    aval[.,i] = ones(2^(i-1),1)
            .*.((0|1).*.ones(2^(nplayer-i),1)) ;               
    i=i+1 ;
  endo ; 
  mstate = zeros(numx,nplayer+1) ;
  mstate[.,1] = zval.*.ones(numa,1) ;
  mstate[.,2:nplayer+1] = ones(numz,1).*.aval ;
  
  @ ------------------------------------------------@  
  @ Matrix with observed indexes of state variables @
  @ ------------------------------------------------@ 
  indzobs = (zobs.==(zval'))*seqa(1,1,numz) ;  
  twop = 2.^(nplayer-seqa(1,1,nplayer)') ;
  indobs = sumc((aobs_1.*twop)') ;  
  indobs = (indzobs-1).*(2^nplayer) + indobs + 1 ;
   
  @ ------------- @
  @ NPL algorithm @
  @ ------------- @ 
  aobs = 1 + reshape(aobs',nobs*nplayer,1) ;
  u0 = zeros(numx*nplayer,kparam) ;
  u1 = zeros(numx*nplayer,kparam) ;
  e0 = zeros(numx*nplayer,1) ;
  e1 = zeros(numx*nplayer,1) ;    
  iter=1 ;
  do while iter<=kiter ;
    "" ;
    " -----------------------------------------------------" ;
    " POLICY ITERATION ESTIMATOR: STAGE =" ;; iter ;
    " -----------------------------------------------------" ;
    "" ;    
    @ -----------------------------------------------------------@
    @ 1. Matrix of transition probabilities Pr(a[t]|z[t],a[t-1]) @
    @ -----------------------------------------------------------@    
    ptrana = (pchoice[.,1].^(aval[.,1]'))
           .*((1-pchoice[.,1]).^(1-aval[.,1]')) ;
    i=2 ;
    do while i<=nplayer ;
      ptrana = ptrana
            .*(pchoice[.,i].^(aval[.,i]'))
            .*((1-pchoice[.,i]).^(1-aval[.,i]')) ;
      i=i+1 ;
    endo ;          
    
    @ --------------------@
    @ 2. Inverse of I-b*F @
    @ --------------------@    
    invi_bf = (ptranz.*.ones(numa,numa)).*(ones(1,numz).*.ptrana);
    invi_bf = crout(eye(numx)-bdisc*invi_bf) ;        
    
    @ -----------------------------------------@
    @ 3. Construction of explanatory variables @
    @ -----------------------------------------@       
    uobs0 = zeros(nobs*nplayer,kparam) ;
    uobs1 = zeros(nobs*nplayer,kparam) ;
    eobs0 = zeros(nobs*nplayer,1) ;
    eobs1 = zeros(nobs*nplayer,1) ;   
    i=1 ;
    do while i<=nplayer ;
      @ --------------------------------------------@
      @ 3.1. Matrices Pr(a[t] | z[t],a[t-1], ai[t]) @
      @ --------------------------------------------@
      mi = aval[.,i]' ;
      ppi = pchoice[.,i] ;
      ppi= (ppi.>=myzero).*(ppi.<=(1-myzero)).*ppi
         + (ppi.<myzero).*myzero
         + (ppi.>(1-myzero)).*(1-myzero) ;      
      ptranai0 = (1-mi)./((ppi.^mi).*((1-ppi).^(1-mi))) ;
      ptranai0 = ptrana .* ptranai0 ;
      ptranai1 = mi./((ppi.^mi).*((1-ppi).^(1-mi))) ;
      ptranai1 = ptrana .* ptranai1 ;
      clear mi ;
               
      @ ------------------------------------@
      @ 3.2 Computing hi = E(ln(Sum(aj)+1)) @
      @ ------------------------------------@
      hi = aval ;
      hi[.,i] = ones(numa,1) ;
      hi = ptranai1 * ln(sumc(hi')) ;
      
      @ -----------------------@
      @ 3.3 Creating U0 and U1 @
      @ -----------------------@
      umat0 = zeros(numx,nplayer+3) ;
      umat1 = eye(nplayer) ;
      umat1 = umat1[i,.] ;
      umat1 = (ones(numx,nplayer).*umat1)
             ~mstate[.,1]~(-hi)~(mstate[.,i+1]-1) ;
      clear hi ;
      
      @ ---------------------------@
      @ 3.4 Creating sumu and sume @
      @ ---------------------------@
      sumu = (1-ppi).*umat0 + ppi.*umat1 ;
      sume = (1-ppi).*(eulerc-ln(1-ppi)) + ppi.*(eulerc-ln(ppi));
      clear ppi ;

      @ ----------------@
      @ 3.5 Creating ww @
      @ ----------------@
      ww = qrtsol(sumu~sume,lowmat(invi_bf)) ;
      ww = qrsol(ww,upmat1(invi_bf)) ;
      clear sumu, sume ;
            
      @ -------------------------------@
      @ 3.6 Creating utilda and etilda @
      @ -------------------------------@
      ptranai0 = (ptranz.*.ones(numa,numa))
               .*(ones(1,numz).*.ptranai0) ;
      utilda0 = umat0 + bdisc*(ptranai0*ww[.,1:kparam]) ;
      etilda0 = bdisc*(ptranai0*ww[.,kparam+1]) ;
      clear umat0, ptranai0 ;
      
      ptranai1 = (ptranz.*.ones(numa,numa))
               .*(ones(1,numz).*.ptranai1) ;
      utilda1 = umat1 + bdisc*(ptranai1*ww[.,1:kparam]) ;
      etilda1 = bdisc*(ptranai1*ww[.,kparam+1]) ;
      clear umat1, ptranai1 ;
            
      @ ----------------------------------------@
      @ 3.7 Creating observations uobs and eobs @
      @ ----------------------------------------@
      uobs0[(i-1)*nobs+1:i*nobs,.] = utilda0[indobs,.] ;
      uobs1[(i-1)*nobs+1:i*nobs,.] = utilda1[indobs,.] ;
      eobs0[(i-1)*nobs+1:i*nobs,.] = etilda0[indobs,.] ;
      eobs1[(i-1)*nobs+1:i*nobs,.] = etilda1[indobs,.] ;      
      u0[(i-1)*numx+1:i*numx,.] = utilda0 ;
      u1[(i-1)*numx+1:i*numx,.] = utilda1 ;
      e0[(i-1)*numx+1:i*numx,.] = etilda0 ;
      e1[(i-1)*numx+1:i*numx,.] = etilda1 ;
      clear utilda0, utilda1, etilda0, etilda1 ;            
      i=i+1 ;
    endo ;
      
    @ ----------------------------------------@
    @ 4. Pseudo Maximum Likelihood Estimation @
    @ ----------------------------------------@                   
    {tetaest,varest} = 
              clogit(aobs,uobs0~uobs1,eobs0~eobs1,namesb) ;                  
    best[.,iter] = tetaest ;
    varb[.,(iter-1)*kparam+1:iter*kparam] = varest ;    
                   
    @ ------------------------- @
    @ 5. Updating probabilities @
    @ ------------------------- @
    i=1 ;
    do while i<=nplayer ;
      buff = (u1[(i-1)*numx+1:i*numx,.]-u0[(i-1)*numx+1:i*numx,.]) ;
      buff = buff*tetaest 
           + (e1[(i-1)*numx+1:i*numx,.]-e0[(i-1)*numx+1:i*numx,.]) ;
      pchoice[.,i] = exp(buff)./(1+exp(buff)) ;
      i=i+1 ;
    endo ;
    
    iter=iter+1 ;
  endo ;
    
  retp(best,varb) ;
endp ;
