Edit on GitHub

Simulations

Here we perform some monte-carlo simulations to illustrute package usage and estimator performance. See Model for an overview of the model. The tests directory contains some additional examples.

Generate Data

using BLPDemand, Statistics, PrettyTables, Printf, Random, LinearAlgebra
K = 2             # number of characteristics
J = 6             # number of products
S = 10            # draws of nu
T = 100           # number of markets
β = ones(K)*2
β[1] = -1.5       # important for equilibrium that higher prices lower sales
σ = ones(K)
σ[1] = 0.2
γ = ones(1)*0.1
Random.seed!(98426)

(sim, ξ, ω) = simulateBLP(J,T, β, σ, γ, S, varξ=0.2, varω=0.2);
@show quantile(vcat((d->d.s[:]).(sim)...), [0, 0.05, 0.5, 0.95, 1])
5-element Vector{Float64}:
 0.00014225499708970805
 0.008799502157570303
 0.07110141519276794
 0.23451158881708067
 0.43509329194476337

Estimation will encounter numerical difficulties if market shares are very close (within 1e-4) to 0 or 1. The above range should be fine.

Instruments

In the simulated data sim[].x[2:end,:] and sim[].w are uncorrelated with ξ and ω, but price, sim[].x[1,:], is endogenous. The price of the jth firm will depend on the characteristics and costs of all other goods, so these are available as instruments. makeivblp is a convenience function that constructs the sum of all other firms' exogenous variables to use as instruments. This is similar to what Berry, Levinsohn, and Pakes (1995) do. We will see below though that much more accurate estimates can be obtained by using the optimal instruments. makeivblp is used by simulateBLP to generate sim[].zd and sim[].zs.

Estimation

We can now estimate the model. Three methods are available: GMM with nested-fixed point (:NFXP), constrained GMM (:MPEC), and constrained GEL (:GEL). See estimateBLP. 

julia> @time nfxp = estimateBLP(sim, method=:NFXP, verbose=false) 56.283514 seconds (109.55 M allocations: 11.298 GiB, 6.01% gc time, 44.85% compilation time)
(β = [-1.249463056038996, 1.9630170845229709], σ = [0.7293030059188311, 0.9143005273895057], γ = [0.1648937304835401], ξ = [[0.23788902011850926, -0.062136337773911726, -0.15067909415155034, 0.14119138247963983, -0.7163793355333339, 0.49881215713832056], [-0.41676959586275647, -0.07433944112616508, -0.1257453903856729, 0.635857359067554, -0.1468921312569782, 0.33102022240193385], [0.37082589353074935, -0.1286539517880545, 0.26321398982898075, 0.4971619664350553, 0.008917924170646074, 0.6376429329612301], [-0.10113883697064452, 0.2064769579346244, -0.5501141316546871, -0.6357933343253537, -0.7715652491041144, 0.4471982923919048], [-1.5858459196431771, 0.37728339849944725, -0.2749046737318408, 0.10077983403525637, -0.4490182086635328, -0.5706464471525794], [0.07469969067444748, -0.2628446459817759, 0.47671742974634634, -0.13730259710824844, -0.4734139879644619, 0.2908014736657875], [0.4374697323130202, 0.6156491350836647, 0.29280145945222236, 0.5055355852931416, 0.7821730713145327, 0.1818656198366646], [-1.2771106796177518, -0.23137697598365747, -0.6064079506077049, 1.134469875158584, -1.029727773074149, -0.7858216505282098], [1.1390410376764846, -0.5136655334584572, 0.29259480630964046, 0.15229775468636275, -0.17380864966339302, 0.655284454950303], [0.27228328990582074, 0.09003529244290587, -0.5613545101746029, -0.23406952974096606, -0.07308156958396106, -0.7357553553101908]  …  [0.08494758118448753, 0.19422157686259445, -0.23002698006203204, 0.43395764046585905, 0.21965551180319887, -0.034550632409312554], [-0.052497427223710424, -0.2590028495386809, -0.09469120427925387, 0.09909283459982048, 0.2905956670472074, 0.2378400159247156], [0.0534959013090987, -0.2078053366601892, 0.4840445620503159, 0.46127934197637166, -0.3360940732189602, -0.6570519175449968], [-0.3914195370867204, 0.21817161229442172, 0.19049075043357894, -0.12889755737856756, -0.4081764987030777, -0.1524334779698322], [-0.610166745715226, -0.47730937126405837, -0.023803346906427048, -0.01476366662501638, -0.9232487039611283, 0.47601225491112886], [0.6981924070468238, -0.3389938653228821, 0.5784059942862644, 0.4552401991145746, 0.3069974864064311, 0.05627424414618698], [0.417655211008765, 0.22274124405332607, 0.12284061430681237, -0.15371576896818562, 0.15838208746614035, 0.1460961051776526], [-1.0098422264142672, -0.06322262785787114, 0.1771545430805963, -0.020522729760088065, 0.03050509103887089, 0.016132294282947712], [-0.3700542321012643, 0.293834123781042, 0.011797832591629698, -0.34019666416025807, -0.7886833465461791, 0.8978988878064149], [0.08582710498774215, -0.46602182862906005, 0.11028250772238724, 0.383240190180024, -0.9225794629815709, -0.5425146795598648]], ω = [[-0.7014042378448027, -0.4878008695779514, -0.7303246863622261, -0.7762201526051153, -0.5951460127844037, -0.44872994572793606], [-0.11887811219701198, 0.6142462081586622, -0.17283844125408948, -0.8709413418084426, -0.08724558010770539, -0.4318665620637805], [0.42492291260038256, 0.4408319256179237, -0.1798786956022574, -0.792680840765003, 0.09006828932757677, -0.5423612491310829], [1.1719578128349868, 0.5353837385049148, -0.08004376657943979, 0.569155293455716, -0.2054178660270042, -0.8599046015937191], [-0.5084927407742841, 0.20737396586631623, 0.642558737153649, 0.4792404229660143, 0.42604370454052903, 0.16018613563583323], [0.19883733164540132, 0.4998309379263953, 0.43661165618900527, 0.3023091532285942, -0.11385686672928275, -0.36853377354353334], [-0.4763673210803104, 0.04245148569922874, 0.8598343893287282, -0.5330777712920816, -0.618742995373333, -0.2959607146802188], [-0.6425556433573792, -0.12819544303944344, 0.14325974261973184, -0.6696438965800281, 0.4455300902879131, -0.3717334236320141], [0.48698458059541516, 0.34432058040498564, 0.02433405559298376, -0.4710375308059086, -0.6782439484413988, 0.3729044848052865], [-0.10481147864597938, -0.46897455250927317, -0.11775451470413996, 0.4622513465020994, 0.6356837547245937, 0.20988356434547806]  …  [0.6851449545075188, -0.16586068356690634, 0.8174394453709385, -0.3765906282339739, 0.0035040079854535, 0.05204275549758433], [-0.3387865155271307, -0.011936297920372735, 0.013493100159258795, -1.0101078946398012, 0.14154578158503925, -0.14902964330030374], [-0.4379675540752608, -0.2530050882470719, 0.022354604453730206, -0.6767621812657975, -0.07970811784173008, -0.04492391071297893], [-0.3417698163427129, -0.5441857868535341, -0.9345003836770915, -0.47483083269584536, -0.2206161866858677, -0.028413222459862983], [0.23094336665155032, 0.11224898231213093, -1.2715607825903197, -0.35537423503752547, -0.0021101646035253108, 0.02867639226383778], [-0.1416410527017636, 0.777167783722683, -0.39867666020056336, -0.36018852721873446, 0.2184861624393425, 0.35541034280997413], [-0.16763690539126208, 0.5600133413157935, -0.19754473202600997, -0.32708790487410244, -0.06538187234766661, 0.4882913733048263], [-1.0839733355935814, 0.6906349023689455, 0.5562911205774486, -0.30783641959896707, 0.764920058832687, -0.17353440735331727], [-0.14988698217979463, -0.07931749019677659, 0.5639838631403117, 1.1209995532528716, 0.83529016910402, 0.08443031277320959], [1.2919989153621974, 0.27981808185955676, -0.24016801954571146, 0.45324008001990196, 0.10761727582979169, 0.17650558689696894]], opt =  * Status: success

 * Candidate solution
    Final objective value:     1.212794e+00

 * Found with
    Algorithm:     L-BFGS

 * Convergence measures
    |x - x'|               = 0.00e+00 ≤ 0.0e+00
    |x - x'|/|x'|          = 0.00e+00 ≤ 0.0e+00
    |f(x) - f(x')|         = 0.00e+00 ≤ 0.0e+00
    |f(x) - f(x')|/|f(x')| = 0.00e+00 ≤ 0.0e+00
    |g(x)|                 = 1.58e-07 ≰ 1.0e-08

 * Work counters
    Seconds run:   31  (vs limit Inf)
    Iterations:    747
    f(x) calls:    3545
    ∇f(x) calls:   3545
)
julia> @time mpec = estimateBLP(sim, method=:MPEC,verbose=true);(β0, σ0, γ0) = ([-1.638020917308739, 2.1619552402599465], [0.05, 0.05], [0.1950352404764217])

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit https://github.com/coin-or/Ipopt
******************************************************************************

This is Ipopt version 3.14.10, running with linear solver MUMPS 5.5.1.

Number of nonzeros in equality constraint Jacobian...:   184207
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:    76807

Total number of variables............................:    21012
                     variables with only lower bounds:     6002
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:    21007
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  0.0000000e+00 4.05e+00 1.85e-03  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  1.9776651e+00 2.50e+02 7.23e+02  -1.0 1.04e+01    -  6.75e-02 8.12e-01h  1
   2  2.1222130e+00 9.25e+01 3.78e+01  -1.0 1.32e+00    -  3.45e-01 9.90e-01h  1
   3  2.3105932e+00 3.38e+01 2.05e+00  -1.0 9.91e-01    -  1.00e+00 1.00e+00h  1
   4  2.3715049e+00 2.21e+01 1.10e+03  -1.0 1.91e+00    -  7.45e-01 4.36e-01h  1
   5  2.3766046e+00 1.13e+01 3.93e+02  -1.0 1.02e+00    -  1.00e+00 6.61e-01h  1
   6  2.4832064e+00 5.71e+00 2.36e+02  -1.0 1.01e+00    -  1.00e+00 6.95e-01h  1
   7  2.6701738e+00 2.43e+00 1.44e+02  -1.0 1.01e+00    -  1.00e+00 8.45e-01h  1
   8  3.2170302e+00 9.23e-01 5.56e+02  -1.0 1.12e+00    -  4.10e-01 6.85e-01H  1
   9  2.7010234e+00 6.07e-01 2.80e+02  -1.0 2.03e+00    -  1.52e-01 3.44e-01F  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10  1.8556488e+00 1.06e+00 1.89e+02  -1.0 2.73e+00    -  4.01e-01 8.29e-01F  1
  11  1.9667907e+00 3.26e-02 3.64e+02  -1.0 1.80e+00    -  2.36e-01 1.00e+00f  1
  12  2.0278899e+00 3.40e+00 3.67e+01  -1.0 2.28e+01  -2.0 1.00e+00 1.00e+00F  1
  13  2.4808682e+00 3.48e+00 5.32e+00  -1.0 6.41e+00    -  1.00e+00 1.00e+00h  1
  14  3.4370802e+00 3.08e+00 8.42e+00  -1.0 1.62e+01    -  1.00e+00 1.00e+00h  1
  15  4.4008926e+00 2.47e+00 4.65e+00  -1.0 1.74e+01    -  1.00e+00 1.00e+00h  1
  16  4.5070605e+00 1.36e-02 5.50e-02  -1.0 3.78e+00    -  1.00e+00 1.00e+00h  1
  17  2.0557480e+00 2.35e+01 6.57e+02  -2.5 7.33e+01    -  4.02e-01 8.72e-01f  1
  18  1.9337156e+00 1.96e+00 5.81e+02  -2.5 1.80e+01    -  1.10e-02 1.00e+00h  1
  19  1.8557082e+00 9.64e-02 7.10e+01  -2.5 1.29e+00    -  8.76e-01 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  20  1.6405131e+00 1.98e-01 3.68e-01  -2.5 2.57e+00    -  1.00e+00 1.00e+00h  1
  21  1.2187716e+00 8.63e-01 1.08e+01  -2.5 5.11e+00    -  5.69e-01 1.00e+00h  1
  22  1.2604002e+00 1.12e-01 2.54e-01  -2.5 6.79e-01    -  1.00e+00 1.00e+00h  1
  23  1.2640633e+00 6.93e-04 3.24e-03  -2.5 2.34e-02    -  1.00e+00 1.00e+00h  1
  24  1.1852316e+00 2.06e-01 1.03e+01  -3.8 1.94e+00    -  7.89e-01 1.00e+00h  1
  25  1.2124428e+00 4.92e-03 4.54e-02  -3.8 1.01e-01    -  1.00e+00 1.00e+00h  1
  26  1.2129470e+00 3.28e-04 8.01e-04  -3.8 6.08e-02    -  1.00e+00 1.00e+00h  1
  27  1.2126202e+00 1.00e-03 1.11e-01  -5.7 1.14e-01    -  9.89e-01 1.00e+00h  1
  28  1.2127941e+00 1.63e-08 5.87e-07  -5.7 5.48e-04    -  1.00e+00 1.00e+00h  1
  29  1.2127941e+00 1.54e-07 1.26e-06  -8.6 1.40e-03    -  1.00e+00 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  30  1.2127941e+00 2.66e-15 1.14e-13  -8.6 9.31e-08    -  1.00e+00 1.00e+00h  1

Number of Iterations....: 30

                                   (scaled)                 (unscaled)
Objective...............:   1.2127941213904363e+00    1.2127941213904363e+00
Dual infeasibility......:   1.1368683772161603e-13    1.1368683772161603e-13
Constraint violation....:   2.6645352591003757e-15    2.6645352591003757e-15
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   2.5059035847836149e-09    2.5059035847836149e-09
Overall NLP error.......:   2.5059035847836149e-09    2.5059035847836149e-09


Number of objective function evaluations             = 37
Number of objective gradient evaluations             = 31
Number of equality constraint evaluations            = 37
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 31
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 30
Total seconds in IPOPT                               = 5.458

EXIT: Optimal Solution Found.
 13.470499 seconds (8.84 M allocations: 613.678 MiB, 3.51% gc time, 61.68% compilation time: 8% of which was recompilation)
julia> @time gel = estimateBLP(sim, method=:GEL,verbose=true);(β0, σ0, γ0) = ([-1.638020917308739, 2.1619552402599465], [0.05, 0.05], [0.1950352404764217])
This is Ipopt version 3.14.10, running with linear solver MUMPS 5.5.1.

Number of nonzeros in equality constraint Jacobian...:   192600
Number of nonzeros in inequality constraint Jacobian.:      600
Number of nonzeros in Lagrangian Hessian.............:    85800

Total number of variables............................:    21605
                     variables with only lower bounds:     6602
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:    21007
Total number of inequality constraints...............:        1
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        1

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0 -2.7631021e+03 5.00e+00 9.98e+01  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1 -2.7652634e+03 4.98e+00 8.44e+01  -1.0 1.01e+01    -  9.49e-02 6.28e-03h  1
   2 -3.6756804e+03 3.67e+02 9.84e+03  -1.0 9.36e+00  -2.0 5.79e-01 9.33e-01h  1
   3 -3.7034617e+03 2.98e+02 8.83e+03  -1.0 1.05e+00    -  8.71e-01 2.07e-01h  1
   4 -3.8408662e+03 1.11e+02 3.65e+03  -1.0 1.04e+00    -  1.00e+00 9.91e-01h  1
   5 -3.8415930e+03 4.06e+01 2.08e+03  -1.0 1.00e+00    -  1.00e+00 1.00e+00h  1
   6 -3.8412275e+03 1.49e+01 7.21e+02  -1.0 9.95e-01    -  8.54e-01 1.00e+00h  1
   7 -3.8411302e+03 5.52e+00 2.83e+02  -1.0 9.95e-01    -  1.00e+00 1.00e+00h  1
   8 -3.8413967e+03 2.02e+00 2.72e+01  -1.0 9.89e-01    -  1.00e+00 1.00e+00h  1
   9 -3.8417957e+03 9.64e-01 6.20e+01  -1.0 9.63e-01    -  4.84e-01 5.85e-01H  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10 -3.8417307e+03 8.72e-01 1.30e+02  -1.0 1.94e+00    -  3.07e-01 9.53e-02f  2
  11 -3.8416210e+03 6.75e-01 2.32e+02  -1.0 1.61e+00    -  7.52e-01 2.29e-01h  2
  12 -3.8414946e+03 5.41e-01 3.05e+02  -1.0 2.09e+00    -  6.95e-01 2.00e-01h  3
  13 -3.8413759e+03 4.38e-01 3.49e+02  -1.0 2.69e+00    -  1.00e+00 1.92e-01h  3
  14 -3.8412670e+03 3.59e-01 2.91e+02  -1.0 3.39e+00    -  9.20e-01 1.82e-01h  3
  15 -3.8411683e+03 2.97e-01 2.22e+02  -1.0 4.24e+00    -  1.00e+00 1.73e-01h  3
  16 -3.8410112e+03 4.78e-01 1.51e+02  -1.0 5.28e+00    -  1.00e+00 3.26e-01f  2
  17 -3.8409154e+03 7.23e-01 1.56e+02  -1.0 8.29e+00    -  9.20e-01 2.68e-01h  2
  18 -3.8409032e+03 9.30e-01 1.64e+02  -1.0 1.17e+01    -  1.00e+00 2.25e-01h  2
  19 -3.8413516e+03 2.33e+00 6.75e+01  -1.0 1.17e+01    -  3.65e-01 6.30e-01F  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  20 -3.8426333e+03 6.46e+00 3.20e+01  -1.0 1.79e+01    -  1.00e+00 1.00e+00f  1
  21 -3.8439386e+03 4.22e+00 1.64e+01  -1.0 2.48e+01    -  1.00e+00 1.00e+00h  1
  22 -3.8435395e+03 3.89e-01 1.64e+00  -1.0 3.20e+00    -  1.00e+00 1.00e+00h  1
  23 -3.8434770e+03 1.18e-02 4.36e-02  -1.0 8.17e-01    -  1.00e+00 1.00e+00h  1
  24 -3.8413086e+03 2.20e+01 5.89e+02  -2.5 7.47e+01    -  3.81e-01 8.03e-01f  1
  25 -3.8409312e+03 1.93e+00 5.81e+02  -2.5 1.65e+01    -  1.12e-02 1.00e+00h  1
  26 -3.8409065e+03 8.60e-02 7.67e+01  -2.5 1.08e+00    -  8.62e-01 1.00e+00h  1
  27 -3.8407588e+03 1.67e-01 9.55e-01  -2.5 2.26e+00    -  1.00e+00 1.00e+00h  1
  28 -3.8401996e+03 9.12e-01 2.44e+01  -2.5 9.03e+00    -  2.19e-01 7.61e-01H  1
  29 -3.8401899e+03 5.64e-02 3.94e+00  -2.5 7.48e-01    -  1.00e+00 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  30 -3.8400490e+03 1.22e-01 4.32e+02  -2.5 2.46e+00    -  1.00e+00 4.42e-01h  2
  31 -3.8398692e+03 1.02e-01 8.98e+00  -2.5 1.16e+00    -  1.00e+00 1.00e+00H  1
  32 -3.8398920e+03 5.95e-03 1.56e+03  -2.5 3.17e-01    -  7.20e-01 1.00e+00h  1
  33 -3.8399039e+03 9.38e-04 6.95e-02  -2.5 9.41e-02    -  1.00e+00 1.00e+00h  1
  34 -3.8399052e+03 1.19e-05 6.67e-04  -2.5 1.10e-02    -  1.00e+00 1.00e+00h  1
  35 -3.8397689e+03 2.96e-01 1.56e+01  -3.8 2.02e+00    -  6.52e-01 8.40e-01f  1
  36 -3.8398315e+03 1.06e-01 5.00e+04  -3.8 9.95e-01    -  7.80e-03 1.00e+00h  1
  37 -3.8398000e+03 2.46e-02 9.61e+01  -3.8 4.37e-01    -  6.53e-01 1.00e+00h  1
  38 -3.8397920e+03 1.35e-02 9.64e-01  -3.8 3.16e-01    -  1.00e+00 1.00e+00h  1
  39 -3.8397930e+03 1.04e-03 2.59e-02  -3.8 8.26e-02    -  1.00e+00 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  40 -3.8397931e+03 2.45e-07 3.61e-06  -3.8 2.23e-03    -  1.00e+00 1.00e+00h  1
  41 -3.8397926e+03 5.06e-04 4.96e-01  -5.7 6.15e-02    -  9.90e-01 1.00e+00h  1
  42 -3.8397927e+03 9.32e-08 9.15e-07  -5.7 1.14e-03    -  1.00e+00 1.00e+00h  1
  43 -3.8397927e+03 7.64e-08 2.89e-06  -8.6 7.64e-04    -  1.00e+00 1.00e+00h  1
  44 -3.8397927e+03 3.55e-15 4.71e-13  -8.6 1.68e-07    -  1.00e+00 1.00e+00h  1

Number of Iterations....: 44

                                   (scaled)                 (unscaled)
Objective...............:   3.8397926919649194e+03   -3.8397926919649194e+03
Dual infeasibility......:   4.7107360050721559e-13    4.7107360050721559e-13
Constraint violation....:   3.5527136788005009e-15    3.5527136788005009e-15
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   2.5059035656602840e-09    2.5059035656602840e-09
Overall NLP error.......:   2.5059035656602840e-09    2.5059035656602840e-09


Number of objective function evaluations             = 97
Number of objective gradient evaluations             = 45
Number of equality constraint evaluations            = 97
Number of inequality constraint evaluations          = 97
Number of equality constraint Jacobian evaluations   = 45
Number of inequality constraint Jacobian evaluations = 45
Number of Lagrangian Hessian evaluations             = 44
Total seconds in IPOPT                               = 7.616

EXIT: Optimal Solution Found.
  8.519893 seconds (2.87 M allocations: 211.235 MiB, 2.16% gc time, 5.53% compilation time)
tbl = vcat(["" "True" "NFXP" "MPEC" "GEL"],
           [(i->"β[$i]").(1:length(β)) β nfxp.β mpec.β gel.β],
           [(i->"σ[$i]").(1:length(σ)) σ nfxp.σ mpec.σ gel.σ],
           [(i->"γ[$i]").(1:length(γ)) γ nfxp.γ mpec.γ gel.γ])
pretty_table(tbl, noheader=true, formatter=ft_printf("%6.3f",3:5))

In the absence of optimization error and other numeric problems, :NFXP and :MPEC should produce identical results. In finite sample, we should expect the GEL estimates to differ. All three methods are consistent, but GEL is also asymptotically efficient (for a fixed choice of z; different z can have large effects on efficiency).

Inference

varianceBLP computes the variance of the estimates produced by either of GMM estimation methods.[simdraws]

[^simdraws:] The variance calculation ignores uncertainty from Monte-Carlo integration. For this to be valid, we must have the number of simulation draws grow faster than the sample size. In our notation, S is the number of simulation draws per observation, so it is sufficient at that $S \to \infty$ as $T \to \infty$.

vnfxp = varianceBLP(nfxp.β, nfxp.σ, nfxp.γ, sim);
vmpec = varianceBLP(mpec.β, mpec.σ, mpec.γ, sim);
nothing

Inference for GEL has not been directly implemented. However, GEL is first order asymptotically equivalent to efficiently weigthed GMM. In other words, GEL estimates have the same asymptotic variance as efficiently weighted GMM. 

v = varianceBLP(gel.β, gel.σ, gel.γ, sim)
vgel = varianceBLP(gel.β, gel.σ, gel.γ, sim, W=inv(v.varm))

f(v) = @sprintf("(%.2f)", norm(sqrt(Complex(v))))
vtbl = permutedims(hcat(tbl[1,:],
                        [hcat(tbl[i+1,:],
                              ["", "", f.([vnfxp.Σ[i,i], vmpec.Σ[i,i], vgel.Σ[i,i]])...])
                         for i in 1:size(vgel.Σ,1)]...
                        ))
pretty_table(vtbl, noheader=true, formatter=ft_printf("%6.3f",3:5))

These reusults look pretty good. In additional experiments (not shown) the results were somewhat sensitive to the simulation design and choice of instruments (see makeivblp).

Optimal Instruments

The above estimators use the unconditional moment restriction

\[E[(\xi, \omega) z] =0.\]

If we assume the conditional moment restriction,

\[E[(\xi, \omega) | z] =0.\]

then, we can potentially use any function of $z$ to form unconditional moments.

\[E[(\xi, \omega) f(z)] =0.\]

The optimal (minimal asymptotic variance) choice of f(z) is

\[\frac{\partial}{\partial \theta} E[(\xi, \omega) | z] =0.\]

optimalIV approximates the optimal instruments by taking as initial estimate of $\theta$, computing $\frac{\partial (\xi, \omega)}{\partial \theta}$ for each observation in the data, and then regressing this on a polynomial function of $z$. Using the fitted values from this regression as $f(z)$ results in much more precise estimates of $\theta$.[z]

sim=optimalIV(mpec.β, max.(mpec.σ, 0.1), # calculating optimal IV with σ near 0 gives poor peformance
              mpec.γ, sim, degree=3);
nothing
julia> @time nfxp = estimateBLP(sim, method=:NFXP, verbose=false) 19.953054 seconds (43.33 M allocations: 5.011 GiB, 6.86% gc time, 0.07% compilation time)
(β = [-1.014806082826089, 2.2026055688140755], σ = [1.338493518142783, 0.010000000009683587], γ = [0.12333794024876811], ξ = [[0.24541729723239225, -0.06209555527330579, -0.1733402330510312, 0.1584858263838117, -0.7403294557982105, 0.4684691707096146], [-0.40858569039050097, -0.235800884676181, -0.1872943562042192, 0.6091864168274445, -0.2112791471110983, 0.27794613415203834], [0.13193890404727515, -0.3710951606126075, 0.11580626446897435, 0.3449503045193919, -0.08721505660378281, 0.4794503209915496], [-0.3715954600810667, 0.09261485590722662, -0.8247878177265723, -0.9214480356048054, -0.976311332959875, 0.3730480414200133], [-1.7378663494005109, 0.12182077807470248, -0.6247835198298504, -0.06565192950485921, -0.7261625434597478, -0.7410458280741634], [-0.07194040784087885, -0.3370504631273652, 0.49611962325535064, -0.4575238192739661, -0.39821348799654, -0.011443699848175004], [0.5919392679944878, 0.7256254567594222, 0.5316539944818774, 0.5746253011243025, 0.8317794636270266, 0.26426649267423], [-0.9781665963999535, -0.011116015578684246, -0.36378115088547025, 1.4059944643340538, -0.7769791543773454, -0.5523814801947907], [1.0156905985358118, -0.7381733926165318, 0.01053037040056967, -0.013089772262896093, -0.5079838944768315, 0.4408800350418629], [0.3763096663986061, 0.19610296041015102, -0.4621499948061365, -0.18259997030836583, -0.07983932911588609, -0.5975582894319581]  …  [0.5147808212165981, 0.4816761172723523, 0.14508760624307548, 0.800966271418844, 0.6358563750316266, 0.037897353383143706], [0.20103633384451958, -0.07205862526044471, 0.04554272713533093, 0.19404056544947296, 0.49100346303331166, 0.4038433553639639], [-0.0314983285123264, -0.2684838075347308, 0.348586714787339, 0.515376641513722, -0.4466369989441218, -0.7659342444061263], [-0.3988377728307061, 0.2551887240951771, 0.2626779504772919, -0.13866308770580374, -0.4332337659598204, -0.1344725351236211], [-0.5789834658895165, -0.5021563808702509, -0.05395295358796476, 0.04607901285220861, -0.9279601420309562, 0.4541118643680243], [0.6556400398258493, -0.6266470754347113, 0.46930018330244183, 0.2303222450164537, 0.002529268544453178, -0.1719252593771885], [0.3974836093164539, 0.11749303741280581, 0.11602303156511096, -0.16895087398312114, 0.1274681746164793, 0.08473495466064734], [-0.8559334190199083, -0.24683823054342557, 0.06854975017090337, -0.001036417708678794, -0.1482217756209332, 0.1053271455732441], [-0.3299391326827413, 0.29808289234481966, -0.013233102691484033, -0.7114636234266212, -0.9164912113170112, 0.8742359553930722], [-0.17138580285357996, -0.48607997713797646, 0.061399049856756704, 0.3375729835688692, -0.9687793783611158, -0.588603794302342]], ω = [[-0.7162952052369547, -0.5101710236471856, -0.7189248434386766, -0.7748187483798694, -0.5945831837887501, -0.46955737451338203], [-0.14094787584593443, 0.6029473795932193, -0.189570314641898, -0.8797969982395286, -0.09301049005886979, -0.46487436510665797], [0.405116083605928, 0.4332770566027053, -0.18831386446422355, -0.8145888408555939, 0.07878734996571182, -0.5752438268021457], [1.1639856055244637, 0.5350867548208311, -0.06112917791949, 0.552370360092581, -0.22060324359758574, -0.8463501702441978], [-0.5419603809493574, 0.17051146722984195, 0.6292436348982529, 0.45750644480089375, 0.3851957806168497, 0.13605741149466716], [0.22009040417974168, 0.5065823354681931, 0.44720185816115043, 0.325712634905581, -0.09529595557742768, -0.3058543491679795], [-0.38544791564761055, 0.128116370400987, 0.8846643820841721, -0.42035161591695713, -0.49244908441573343, -0.2250659068471654], [-0.5753329521863201, -0.08587978148779107, 0.14276949809623105, -0.5924923216454628, 0.43603520783790645, -0.3316025014612532], [0.5201922003462494, 0.37999243770972846, 0.06215874027506176, -0.39136410077777467, -0.5734772174301778, 0.41399249665666626], [-0.11230533138865392, -0.4757185207789143, -0.08829413020493976, 0.44749071617973923, 0.6264191830144014, 0.17861036057050164]  …  [0.6533442966974748, -0.19013559242620062, 0.8020709789516944, -0.3880572654613293, 0.009027945068524904, 0.01647552544533529], [-0.2715889732606759, 0.046112208569858065, 0.04451837564732735, -0.8276414693322331, 0.19085327072183203, -0.09030731324887574], [-0.4728778380853805, -0.28611137423055594, -0.01740400697079171, -0.6697529354359991, -0.12585290915560585, -0.08834401075008255], [-0.3489534770208378, -0.5456082128735474, -0.8041188151332052, -0.461875817270165, -0.21550969201079873, -0.03832748510836329], [0.2402620707430901, 0.11836208777567908, -1.0815403502289873, -0.32731355082020624, 0.0325790957119462, 0.03577698429281801], [-0.10146391062940023, 0.7773466006600878, -0.334059485637359, -0.3236418875179034, 0.2500457183134002, 0.36157564469534514], [-0.1878649055988115, 0.5711899038295858, -0.1983823967959854, -0.3265313766216949, -0.060521065710365175, 0.46686311022616844], [-1.128984323899577, 0.6742401883381421, 0.5157577278576678, -0.3343449891843988, 0.7426913828598399, -0.2458629989491154], [-0.16009745609647585, -0.10907783530503787, 0.5402414853479542, 1.1106624356270873, 0.8068970754426392, 0.06300391578752128], [1.3042517907584683, 0.27714785138603065, -0.19492729427100491, 0.4593356988267128, 0.14507414072003963, 0.19372160927015247]], opt =  * Status: failure (line search failed)

 * Candidate solution
    Final objective value:     6.840687e-04

 * Found with
    Algorithm:     L-BFGS

 * Convergence measures
    |x - x'|               = 2.22e-16 ≰ 0.0e+00
    |x - x'|/|x'|          = 8.76e-18 ≰ 0.0e+00
    |f(x) - f(x')|         = 2.71e-18 ≰ 0.0e+00
    |f(x) - f(x')|/|f(x')| = 3.96e-15 ≰ 0.0e+00
    |g(x)|                 = 6.94e-05 ≰ 1.0e-08

 * Work counters
    Seconds run:   19  (vs limit Inf)
    Iterations:    654
    f(x) calls:    1835
    ∇f(x) calls:   1835
)
julia> @time mpec = estimateBLP(sim, method=:MPEC,verbose=false);  7.667301 seconds (2.40 M allocations: 178.270 MiB, 0.94% gc time)
julia> @time gel = estimateBLP(sim, method=:GEL,verbose=false);  7.513637 seconds (2.51 M allocations: 186.566 MiB, 1.08% gc time)
vnfxp = varianceBLP(nfxp.β, nfxp.σ, nfxp.γ, sim)
vmpec = varianceBLP(mpec.β, mpec.σ, mpec.γ, sim)
v = varianceBLP(gel.β, gel.σ, gel.γ, sim)
vgel = varianceBLP(gel.β, gel.σ, gel.γ, sim, W=inv(v.varm))

f(v) = @sprintf("(%.2f)", norm(sqrt(Complex(v))))
tbl = vcat(["" "True" "NFXP" "MPEC" "GEL"],
           [(i->"β[$i]").(1:length(β)) β nfxp.β mpec.β gel.β],
           [(i->"σ[$i]").(1:length(σ)) σ nfxp.σ mpec.σ gel.σ],
           [(i->"γ[$i]").(1:length(γ)) γ nfxp.γ mpec.γ gel.γ])
vtbl = permutedims(hcat(tbl[1,:],
                        [hcat(tbl[i+1,:],
                              ["", "", f.([vnfxp.Σ[i,i], vmpec.Σ[i,i], vgel.Σ[i,i]])...])
                         for i in 1:size(vgel.Σ,1)]...
                        ))
pretty_table(vtbl, noheader=true, formatter=ft_printf("%6.3f",3:5))

We see that with the optimal instruments all three methods produce essentially the same results (in this case, theoretically, NFXP=MPEC, and both are asymptotically equivalent to GEL, so this is expected). Moreover, the estimates are now much more precise and quite close to the true parameter values.

Calculating Elasticities

Using share and the ForwardDiff.jl package, we can calculate elasticities of demand with respect to each characteristic.

using ForwardDiff

function elasticity(β, σ, γ, dat, ξ)
  T = length(dat)
  K,J = size(dat[1].x)
  etype = typeof(dat[1].x[:,1]'*β+ξ[1][1])
  e = Array{Array{etype,3},1}(undef, T)
  for t in 1:T
    @views xt = dat[t].x
    st(x) = share(x'*β + ξ[t], σ, x, dat[t].ν)
    ∂s = reshape(ForwardDiff.jacobian(st, xt), J, K, J)
    s = st(xt)
    e[t] = zeros(etype, J,K,J)
    for k in 1:K
      for j in 1:J
        e[t][:,k,j] .= ∂s[:,k,j]./s .* xt[k,j]
      end
    end
  end
  return e
end

function avg_price_elasticity(β,σ,γ, dat)
  ξ = Vector{Vector{eltype(β)}}(undef, T)
  for t in 1:T
    ξ[t] = delta(dat[t].s, dat[t].x, dat[t].ν, σ) - dat[t].x'*β
  end
  e=elasticity(β,σ,γ,dat,ξ)
  avge = zeros(eltype(e[1]),size(e[1]))
  # this assumes J is the same for all markets
  for t in 1:T
    avge .+= e[t]
  end
  avge /= T
  avge[:,1,:]
end

price_elasticity = avg_price_elasticity(mpec.β, mpec.σ, mpec.γ, sim)
6×6 Matrix{Float64}:
 -2.9088     0.22363    0.262108   0.272599   0.229125   0.26681
  0.219268  -2.83299    0.255089   0.287176   0.232711   0.254208
  0.230881   0.224091  -2.77334    0.275359   0.238741   0.263066
  0.218435   0.228432   0.26313   -2.54149    0.231777   0.257383
  0.220412   0.227279   0.269268   0.287435  -2.75028    0.251296
  0.227462   0.228777   0.268632   0.272097   0.229481  -2.67363

Standard errors of elasticities and other quantities calculated from estimates can be calculated using the delta method.

D = ForwardDiff.jacobian(θ->avg_price_elasticity(θ[1:K], θ[(K+1):(2K)], θ[(2K+1):end], sim), [mpec.β;mpec.σ;mpec.γ])
V = D*vmpec.Σ*D'
se = reshape(sqrt.(diag(V)), size(price_elasticity))

tbl = Array{Any, 2}(undef, 2J+1, J+1)
tbl[1,1]=""
for j in 1:J
  tbl[1,j+1] = "price $j"
  tbl[2j,1] = "share $j"
  tbl[2j+1,1] = ""
  for l in 1:J
    tbl[2j, l+1] = price_elasticity[j,l]
    tbl[2j+1, l+1] = @sprintf("(%.3f)",se[j,l])
  end
end
pretty_table(tbl, noheader=true, formatter=ft_printf("%6.3f", 2:(J+1)))

The table above shows the estimated average elasticity of each good's share with respect to each good's price. Standard errors are in parentheses.

  • zThe choice of z still matters when using optimalIV. For firm j, the values of x[2:end,l,t] (for both l=j and l != j) are all potential instruments for x[1,j,t]. However, makeivblp does not use all of these. It instead uses their sum and the sum of exp(-(x[2:end,j,t] - x[2:end,l,t])^2). Adjusting these choices might give better results.