Math test

[John Shannon]

My likelihood function is of the form:
\[ L_t = N_t \ln\theta + \ln\Gamma\left(\frac{1}\theta+N_t\right)-\ln\Gamma\left(\frac{1}\theta\right)+ \sum_{i=1}^{N _t}\ln\left(\rho\lambda\left(\frac{b_{it}}\lambda \right)^{\rho-1}\right) +\left(\frac{1}\theta+N_t\right)\ln\left( 1+\theta\sum_{i=1}^{N_t}\left(\frac{b_{it}}\lambda \right)^\rho\right) \]
The problem is that N_t, the number of bidders in auction t appears in one of the parameters I want to estimate via maximum likelihood:
\[ \rho=\gamma_1N_t+\gamma_0 \]
This means that N_t shows up in my ml program as both a dependent variable and an independent variable

program gammaweibull3
     args lnf lambda rho theta
     quietly replace `lnf' = 3 * ln(`theta') + lngamma(1/`theta' + 3 ) - lngamma(1/`theta') + ln(`rho'*`lambda'*($ML_y1 /`lambda')^(`rho'-1)) + ln(`rho'*`lambda'*($ML_y2 /`lambda')^(`rho'-1)) + ln(`rho'*`lambda'*($ML_y3 /`lambda')^(`rho'-1)) + (1/`theta' + 3 )*ln(1+`theta'*(($ML_y1 /`lambda')^(`rho') + ($ML_y2 /`lambda')^(`rho') + ($ML_y3 /`lambda')^(`rho')))
end

ml model lf gammaweibull3 (lambda: bid1 bid2 bid3 = nbidders ...) (rho:) (theta:)
ml maximize

which gives me this output:

initial:       log likelihood =  1796499.8
rescale:       log likelihood =  1796499.8
rescale eq:    log likelihood =  3.88e+298
could not calculate numerical derivatives -- discontinuous region with missing
values encountered
r(430);

[John Shannon]

This is my first time using Stata for maximum likelihood estimation, and I'm encountering a problem where my ml evaluator will not converge.

I am using a parametric specification for bid distributions based on Athey, Levin, and Seira (2011). My dataset has one observation per auction with variables for the value of the highest bid, second highest bid, etc. These variables have missing values if a bid of that rank was not submitted in that auction.

The log-likelihood function for auction t is:
\[ \ln L_t\left(\rho,\lambda,\theta\right) = N_t \ln\theta + \ln\Gamma\left(\frac{1}\theta+N_t\right)-\ln\Gamma\left(\frac{1}\theta\right)+\sum_{i=1}^{N _t}\ln\left(\rho\lambda\left(\frac{b _{it}}\lambda\right)^{\rho-1}\right)+\left(\frac{1}\theta+N_t\right)\ln\left( 1+\theta\sum_{i=1}^{N_t}\left(\frac{b_{it}}\lambda \right)^\rho\right) \]

N_t is the number of bidders in auction t and gamma is the gamma function.

Here's the ml program I wrote. $ML_y1 through $ML_y12 are the top 12 bids, and $ML_y13 is the number of bidders:

program gammaweibull
 args lnf lambda rho theta
 tempvar ncons nmult lsumpdf sumcdf lsumcdf ///
   ld1 ld2 ld3 ld4 ld5 ld6 ld7 ld8 ld9 ld10 ld11 ld12
 quietly {
  gen double `ncons' = $ML_y13 * ln(`theta') + lngamma(1/`theta' + $ML_y13 ) - lngamma(1/`theta')
  gen double `nmult' = 1/`theta' + $ML_y13
  
  *If bid values are missing from an auction, then the program observes ln(0)...
  *To correct for this, manually replace additional term with 0 for missing observations
  gen double `ld1'  = ln(`rho'*`lambda'*($ML_y1 /`lambda')^(`rho'-1))
  gen double `ld2'  = `ld1'  + ln(`rho'*`lambda'*($ML_y2  /`lambda')^(`rho'-1))
  replace    `ld2'  = `ld1'  if $ML_y2  ==.
  gen double `ld3'  = `ld2'  + ln(`rho'*`lambda'*($ML_y3  /`lambda')^(`rho'-1))
  replace    `ld3'  = `ld2'  if $ML_y3  ==.
  gen double `ld4'  = `ld3'  + ln(`rho'*`lambda'*($ML_y4  /`lambda')^(`rho'-1))
  replace    `ld4'  = `ld3'  if $ML_y4  ==.
  gen double `ld5'  = `ld4'  + ln(`rho'*`lambda'*($ML_y5  /`lambda')^(`rho'-1))
  replace    `ld5'  = `ld4'  if $ML_y5  ==.
  gen double `ld6'  = `ld5'  + ln(`rho'*`lambda'*($ML_y6  /`lambda')^(`rho'-1))
  replace    `ld6'  = `ld5'  if $ML_y6  ==.
  gen double `ld7'  = `ld6'  + ln(`rho'*`lambda'*($ML_y7  /`lambda')^(`rho'-1))
  replace    `ld7'  = `ld6'  if $ML_y7  ==.
  gen double `ld8'  = `ld7'  + ln(`rho'*`lambda'*($ML_y8  /`lambda')^(`rho'-1))
  replace    `ld8'  = `ld7'  if $ML_y8  ==.
  gen double `ld9'  = `ld8'  + ln(`rho'*`lambda'*($ML_y9  /`lambda')^(`rho'-1))
  replace    `ld9'  = `ld8'  if $ML_y9  ==.
  gen double `ld10' = `ld9'  + ln(`rho'*`lambda'*($ML_y10 /`lambda')^(`rho'-1))
  replace    `ld10' = `ld9'  if $ML_y10 ==.
  gen double `ld11' = `ld10' + ln(`rho'*`lambda'*($ML_y11 /`lambda')^(`rho'-1))
  replace    `ld11' = `ld10' if $ML_y11 ==.
  gen double `ld12' = `ld11' + ln(`rho'*`lambda'*($ML_y12 /`lambda')^(`rho'-1))
  replace    `ld12' = `ld11' if $ML_y12 ==.
  gen `lsumpdf' = `ld12'
  
  *Now if the bid is 0 then the additional terms will be 0, so we only need one variable
  gen double `sumcdf'  = ($ML_y1  / `lambda')^(`rho') ///
           + ($ML_y2  / `lambda')^(`rho') ///
           + ($ML_y3  / `lambda')^(`rho') ///
           + ($ML_y4  / `lambda')^(`rho') ///
           + ($ML_y5  / `lambda')^(`rho') ///
           + ($ML_y6  / `lambda')^(`rho') ///
           + ($ML_y7  / `lambda')^(`rho') ///
           + ($ML_y8  / `lambda')^(`rho') ///
           + ($ML_y9  / `lambda')^(`rho') ///
           + ($ML_y10 / `lambda')^(`rho') ///
           + ($ML_y11 / `lambda')^(`rho') ///
           + ($ML_y12 / `lambda')^(`rho')  
  gen double `lsumcdf'  = ln(1+ `theta'*`sumcdf' )
  
  replace `lnf' = `ncons'+`lsumpdf'+(`nmult'*`lsumcdf')
  }
end

When I use ml maximize I get this output:

. ml maximize

initial:       log likelihood =  134279.26
rescale:       log likelihood =  134279.26
rescale eq:    log likelihood =  141848.78
Iteration 0:   log likelihood =  141848.78  (not concave)
Iteration 1:   log likelihood =  141869.98  (not concave)
Iteration 2:   log likelihood =   141881.5  (not concave)
Iteration 3:   log likelihood =  141882.65  (not concave)
Iteration 4:   log likelihood =  141882.76  (not concave)
Iteration 5:   log likelihood =  141882.81  (not concave)
cannot compute an improvement -- discontinuous region encountered
r(430);

I get similar results when I use the "difficult" option. I tried writing a smaller version of the program for auctions with only three bidders and ran it on that sample and got the same results. I used ml query to look for the problem, and found this:

Current status
    Coefficient values
        1:                                   -1.5495e-08
        2:                                   -.000023924
        3:                                   265.716277
        4:                                   .000022557
        5:                                   108.852769
        6:                                   3.3868e+17
    Function value:                          .
    Converged:                               no

Coefficient 6 corresponds to theta. I've tried everything I can think of to achieve convergence, but haven't had any luck.

One other problem I have is that I want to use the number of bidders as an independent variable in one of my parameter equations, but Stata automatically omits it because of multicollinearity. Given the form of the likelihood function I don't see how multicollinearity would be a problem, but it could be that I'm missing something. Is there a way to get around this?