/*
**  MPEPROB.SRC     Procedure that computes players' choice probabilities
**                  of Markov Perfect Equilibrium in a dynamic
**                  game of firms' market entry/exit with incomplete information
**                  Gauss-Seidel iterations
**
**  by Victor Aguirregabiria
**  First version: August, 2003
**  Last version: November, 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 active:  U(0) = e0[it]
**
**  Active:     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:
**  { prob, psteady, mstate, dconv } = 
**     mpeprob(theta,xval,ptranx)
**
**  Inputs:
**      theta   - Vector of parameters
**                (alpha|beta|delta|thetaE|disfact|sigmaeps)
**                alpha = (nplayer x 1) vector with players' fixed effects
**                beta  = Parameter associated with x[t]
**                delta = Parameter of strategic interactions.
**                thetaE = Cost of entry
**                dfact  = Time discount factor
**                sigma  = Standard deviation of epsilons
**      zval    - Vector with values of z[t]
**      ptranz  - Transition probability of z[t]   
**
**  Outputs:
**      prob   - (numx x nplayer) vector with MPE probs of entry
**               The states in the rows of prob are ordered as follows.
**               Example: zval=(1|2) and 3 players:
**                          z[t]    a[1,t-1]    a[2,t-1]    a[3,t-1]
**               Row 1:     1           0           0           0
**               Row 2:     1           0           0           1
**               Row 3:     1           0           1           0
**               Row 4:     1           0           1           1
**               Row 5:     1           1           0           0
**               Row 6:     1           1           0           1
**               Row 7:     1           1           1           0
**               Row 8:     1           1           1           1
**               Row 9:     2           0           0           0
**               Row 10:    2           0           0           1
**               Row 11:    2           0           1           0
**               Row 12:    2           0           1           1
**               Row 13:    2           1           0           0
**               Row 14:    2           1           0           1
**               Row 15:    2           1           1           0
**               Row 16:    2           1           1           1
**
**      psteady - (numx x 1) vector with steady-state 
**                distribution of {z[t],a[t-1]}
**
**      mstate  - (numx x (nplayer+1)) matrix with values 
**                of state variables {z[t],a[t-1]}. 
**                The states in the rows are ordered as 
**                in the matrix "prob" decribed above.
**
**      dconv   - Indicator for convergence: 
**                  dconv = 1 ===> convergence
**                  dconv = 0 ===> no convergence
**
*/

proc (4) = mpeprob(theta,zval,ptranz);
  local alpha, beta, delta, thetae, dfact, 
        sigma, nplayer, numz, numa, numx, aval, mstate, 
        i, prob0, critconv, criter, dconv, iter, prob1,         mi, ppi, ptranai0, ptranai1, hi,                       profit1, profit0, ptrana, iptran0, iptran1,        
        j, buff, v0, cbell, v1, maxv1, 
        psteady0, psteady1, ptran ;
        
  numz = rows(zval) ;  
  nplayer = rows(theta)-5 ;
  numa = 2^nplayer ;
  numx = numz * numa ;

  alpha  = theta[1:nplayer] ;
  beta   = theta[nplayer+1] ;
  delta  = theta[nplayer+2] ;
  thetae = theta[nplayer+3] ;
  dfact  = theta[nplayer+4] ;
  sigma  = theta[nplayer+5] ;
  
  @ ---------------------------------------------------@
  @ 1. Matrices 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 ;

  @ --------------------------------------- @
  @ 2. Initializing vector of probabilities @
  @ --------------------------------------- @
  prob0 = 0.5 * ones(numx,nplayer) ;
  critconv = (1e-3)*(1/numx) ;
  criter = 1000 ;
  dconv = 1 ;

  @ ---------------------- @
  @ 3. Iterative algorithm @
  @ ---------------------- @
  iter=1 ;
  do while (criter>critconv) ;
    "Iteration =" ;; iter ;; "Criterion =" ;; criter ;
    prob1 = prob0 ; 
    @ ----------------------------------------------------@
    @ 3.1 Matrix of transition probs Pr(a[t]|z[t],a[t-1]) @
    @ ----------------------------------------------------@
    ptrana = (prob1[.,1].^(aval[.,1]'))
           .*((1-prob1[.,1]).^(1-aval[.,1]')) ;
    i=2 ;
    do while i<=nplayer ;
      ptrana = ptrana
            .*(prob1[.,i].^(aval[.,i]'))
            .*((1-prob1[.,i]).^(1-aval[.,i]')) ;
      i=i+1 ;
    endo ;        
        
    i=1 ;    
    do while i<=nplayer ;
      @ ----------------------------------------@
      @ 3.2 Matrices Pr(a[t]|z[t],a[t-1],ai[t]) @
      @ ----------------------------------------@
      mi = aval[.,i]' ;
      ppi = prob1[.,i] ;
      iptran0 = (1-mi)./((ppi.^mi).*((1-ppi).^(1-mi))) ;
      iptran0 = ptrana .* iptran0 ;
      iptran1 = mi./((ppi.^mi).*((1-ppi).^(1-mi))) ;
      iptran1 = ptrana .* iptran1 ;
      clear mi ;
               
      @ ------------------------------------@
      @ 3.3 Computing hi = E(ln(Sum(aj)+1)) @
      @ ------------------------------------@
      hi = aval ;
      hi[.,i] = ones(numa,1) ;
      hi= iptran1 * ln(sumc(hi')) ;
      
      @ ---------------------------------------- @
      @ 3.4 Matrices with expected profit firm i @
      @ ---------------------------------------- @
      profit1 = alpha[i] + beta * mstate[.,1]
              - thetaE*(1-mstate[.,i+1])
              - delta * hi ;              
      profit0 = zeros(numx,1) ;
                                
      @ ---------------------------------------------------- @      
      @ 3.5 Transition probabilities for firm i              @
      @         Pr(x[t+1],a[t] | x[t],a[t-1],ai[t]=0)        @
      @         Pr(x[t+1],a[t] | x[t],a[t-1],ai[t]=1)        @      
      @ ---------------------------------------------------- @      
      iptran0 = (ptranz.*.ones(numa,numa))
              .*(ones(1,numz).*.iptran0);
      iptran1 = (ptranz.*.ones(numa,numa))
              .*(ones(1,numz).*.iptran1);
      
      @ -------------------------------------- @
      @ 5. Computing value function for firm i @
      @ -------------------------------------- @      
      v0 = zeros(numx,1) ;
      cbell = 1000 ;
      do while cbell>critconv ;      
        v1 = (profit0 + dfact * iptran0 * v0)
           ~ (profit1 + dfact * iptran1 * v0) ;
        v1 = v1./sigma ;
        maxv1 = maxc(v1') ;
        v1 = v1 - maxv1 ;
        v1 = sigma * ( maxv1 + ln(exp(v1[.,1]) + exp(v1[.,2])) ) ;
        cbell = maxc( abs(v1-v0) ) ;
        v0 = v1 ;
      endo ;
      
      @ ------------------------------------ @
      @ 6. Updating probabilities for firm i @
      @ ------------------------------------ @      
      v1 = (profit0 + dfact * iptran0 * v0)
         ~ (profit1 + dfact * iptran1 * v0) ;      
      v1 = v1./sigma ;
      maxv1 = maxc(v1') ;
      v1 = v1 - maxv1 ;    
      prob1[.,i] = exp(v1[.,2])./(exp(v1[.,1])+exp(v1[.,2])) ;                   
      i=i+1 ;
    endo ;    

    criter = maxc(maxc(abs(prob1-prob0))) ;
    prob0 = prob1 ;
    iter = iter + 1 ;
  endo ;  
  clear iptran0, iptran1, v0, v1 ;

  if (criter>critconv) ;
    dconv = 0 ;
    psteady1 = 0 ;    
  else ;
    ptrana = (prob1[.,1].^(aval[.,1]'))
           .*((1-prob1[.,1]).^(1-aval[.,1]')) ;
    i=2 ;
    do while i<=nplayer ;
      ptrana = ptrana
            .*(prob1[.,i].^(aval[.,i]'))
            .*((1-prob1[.,i]).^(1-aval[.,i]')) ;
      i=i+1 ;
    endo ;
    ptrana = (ptranz.*.ones(numa,numa))
           .*(ones(1,numz).*.ptrana) ;   
    criter = 1000 ;
    psteady0 = (1/numx) * ones(numx,1) ;
    do while criter>critconv ;      
      psteady1 = ptrana' * psteady0 ;
      criter = maxc( abs(psteady1-psteady0) ) ;
      psteady0 = psteady1 ;
    endo ;      
  endif ;  
    
  retp(prob1,psteady1,mstate,dconv) ;
endp ;
