Update: I have added a post on importance sampling here.

Quick note: This code was written in Julia 0.6 and so will not be compatible with some aspects of my previous posts.

I recently replicated Petrin and Seo’s 2016 paper on using optimization as a source of identification in the random coefficients framework. Their method is more computationally complex than BLP, making efficiency a primary concern. Because it is built off of BLP’s basic methodology (it still uses the contraction mapping and calculates marginal costs analogously), I had to go back and streamline my baseline random coefficients estimator. This page only introduces the new code and what changed from my previous version. For a detailed discussion of the mechanics of BLP, see my previous post and the resources included there.

###### Data Set-up

I have to include the data-cleaning portion of my code as the previous version is not compatible with Julia 0.6. However, it is nearly identical to what I have previously posted, so you only need to worry about this to the extent that the import and cleaning are no longer working.

using DataFrames using DataFramesMeta using Base.LinAlg.BLAS using CSV # Read Data using BenchmarkTools using NLopt using Distributions # For Halton Draws using StatsFuns using Primes cars = CSV.read("[filepath]/Automobile Data.csv", header = true, allowmissing = :none); cars[:ln_hpwt] = log.(cars[:hpwt]) cars[:ln_space] = log.(cars[:space]) cars[:ln_mpg] = log.(cars[:mpg]) cars[:ln_mpd] = log.(cars[:mpd]) cars[:ln_price] = log.(cars[:price]) cars[:trend] = cars[:market] - 1 + 70 cars[:cons] = 1; regSet = @linq cars |> @by(:model_year, s_0 = log.(1 - sum(:share))) regSet = join(cars, regSet, on = :model_year); regSet = @linq regSet |> @transform(s_i = log.(:share)) |> @transform(dif = :s_i - :s_0); regSet[:dif_2] = log.(regSet[:share]) - log.(regSet[:share_out]) ; regSet[:ln_price] = log.(regSet[:price]); regSet = sort!(regSet, [:market, :firmid]) markets = convert(Matrix, regSet[:, [ :market] ]); marks = unique(markets) firms = convert(Matrix, regSet[:, [ :firmid] ]); modlist = unique(regSet[:newmodv]) modLoc = regSet[:newmodv]; X = convert(Array{Float64,2}, regSet[:, [ :cons, :hpwt, :air, :mpd, :space] ]) # Demand Variables W = convert(Array{Float64,2}, regSet[:, [ :cons, :ln_hpwt, :air, :ln_mpg, :ln_space, :trend] ]) # Supply Variables price = convert(Array{Float64,2}, regSet[:, [ :price] ]) # Price delta_0 = convert(Array, round.(regSet[:,:dif_2], 20)); βs = [-7.061, 2.883, 1.521, -0.122, 3.46] γs = [0.952, 0.477, 0.619, -0.415, -0.046, 0.019] α = 43.501 σs = [3.612, 4.628, 1.818, 1.050, 2.056];

Here’s where the changes actually start to matter. My previous function for computing marginal costs was very inefficient because it performed a lot of unnecessary calculations. In it, I calculated the derivative of shares with respect to price between all models in the market. However, it is only necessary to calculate these derivatives for models owned by the same company (e.g. Ford only cares about the substitution between cars it owns, and not about the substitution to other brands). The below code determines the start and end rows for each firm, so that I can efficiently index the data and perform the relevant calculations for each firm individually.

J_f = Array{Int64}(0,2) count = 0 for ms = 1:20 midx = find(markets .== ms) tfirms = unique(firms[midx]) for f in tfirms J_f = vcat(J_f, [count + findfirst(firms[midx], f) count + findlast(firms[midx], f)]) end count += length(midx) end;

Similarly, I find the first and last index of each market in order to more efficiently index the data. However, there isn’t that much of a speed difference between using the find function once vs. thousands of times.

# Beginning and end of the market ids idxs = [findfirst(markets, i) for i = 1:20] idxe = [findlast(markets, i) for i = 1:20] midx = hcat(idxs, idxe);

###### Model Type

As always, I define a custom type to hold all the variables I will use throughout estimation. This way nothing will be defined in global scope (in Julia, using variables in a function that are defined globally can seriously slow down your program).

type ModelData price::Array{Float64, 1} X::Array{Float64, 2} alpha::Float64 beta::Array{Float64, 1} gamma::Array{Float64, 1} sigma::Array{Float64, 1} delta::Array{Float64, 1} xi::Array{Float64, 1} marketidx::Array{Int64, 1} # Market IDs in a vector midx::Array{Int64, 2} # Beginning and end index for each market J_f::Array{Int64, 2} # Beginning and end index for each firm's products y_it::Array{Float64, 2} # Random income draws v_ik::Array{Float64, 2} # Random coefficient draws unobs_weight::Array{Float64, 1} p_ijt::Array{Float64, 2} # Array for individual choice probability s_jt::Array{Float64, 1} # Array for implied market share act_s::Array{Float64, 1} # Array for actual market shares mc::Array{Float64, 1} # Array for marginal costs optinst::Array{Float64,2} # Array of Optimal Instruments W::Array{Float64, 2} function ModelData( price , X , alpha , beta , gamma, sigma, delta, xi, marketidx, midx, J_f, y_it, v_ik, unobs_weight, p_ijt, s_jt, act_s, mc, optinst) return new( price , X , alpha , beta , gamma, sigma, delta, xi, marketidx, midx, J_f, y_it, v_ik, unobs_weight, p_ijt, s_jt, act_s, mc, optinst, zeros(X) ) end end

###### Generating Instruments

The instruments are generated in the same manner as before. I include two options for calculating the instruments. First, when normal == 1, they are calculated consistently with BLP 1995 (based on the mistake they made in their code), and can be used to recreate their Table III. Otherwise, the instruments are calculated following their description in Section 5.1 (for more information, see the discussion in their working paper as well).

function gen_inst( inX, normal = 1 ) totMarket = similar(inX) totFirm = similar(inX) for m in marks sub = inX[find(markets .== m),:] firminfo = firms[find(markets .== m),:] #modelinfo = modLoc[find(markets .== m),:] sameFirm = convert(Array{Float64,2}, firminfo .== firminfo') #sameFirm = sameFirm - diagm(diag(sameFirm)) #sameProduct = ones(sameFirm) - diagm(ones(size(sub,1), 1)[:]) z_1 = similar(sub) for i = 1:size(sub, 2) if normal == 1 z_1[:,i] = sum((sub[:,i] .* sameFirm), 1)' # Prime on sameFirm is key else z_1[:,i] = sum((sub[:,i] .* (sameFirm .- eye(sameFirm)))', 1)' # Prime on sameFirm is key end end totFirm[find(markets .== m),:] = z_1 # Within Market sub = inX[find(markets .== m),:] z_1 = similar(sub) sameFirm = firminfo .== firminfo' for i = 1:size(sub, 2) if normal == 1 z_1[:,i] = sum(sub[:,i] .* (sameFirm .+ .!sameFirm ),1) else z_1[:,i] = sum((sub[:,i] .* (.!sameFirm )),1) end end totMarket[find(markets .== m),:] = z_1 end return [totFirm, totMarket] end Z_x_1, Z_x_2 = gen_inst(X,0) Z_w_1, Z_w_2 = gen_inst(W,0) Z_x = hcat(X, Z_x_1, Z_x_2, convert(Array{Float64,1}) Z_w = hcat(W, Z_w_1, Z_w_2, convert(Array{Float64,1}, regSet[:mpd])) Z_ = [Z_x zeros(Z_w); zeros(Z_x) Z_w];

###### Random Draws

Following a recommendation in Skraina and Judd, I use different draws for each market. BLP used the same set of draws for each market, which is not ideal because the simulation error will be correlated across markets. By using a different set of draws, the simulation error will tend to cancel out, resulting in more precise estimates. Because we use the same number of draws in each market, this we get more precise estimates for little additional computational cost. I scale the number of draws and the weights by the number of markets (in this case 20) to ensure that the set of draws can be partitioned by market.

srand(41318) ns = 1000*20 v_ik1 = randn(6, ns )' m_t = repeat(incomeMeans, inner = [ns, 1]) #y_it = exp(m_t + sigma_v * repeat(v_ik[:,end], outer = [length(incomeMeans),1])); y_it1 = exp.(incomeMeans .+ sigma_v * v_ik1[:,end]'); unobs_weight = ones(ns)'/(ns)*20;

###### Contraction Mapping and Marginal Costs

Here is the bulk of the code. I have combined the calculation of the contraction mapping and the marginal costs into a single function. While I have written extensively about the bulk of this code already (again, here) there is one significant change. Raynaerts, Varadhan, and Nash (2012) note that the contraction mapping can also be formulated as a root-finding problem. With this insight, they explore three methods of solving for the : Newton’s Method, SQUAREM, and DF-SANE. Having implemented each of these methods, I found that a combination of the contraction mapping and Newton’s method worked best for this problem. Newton’s method has quadratic convergence but is unstable from starting points far away from the solution. I therefore use the contraction mapping to find a region of attraction and then switch to Newton’s method to quickly converge to the solution (for those of you familiar with John Rust’s bus engine replacement model, this is almost identical to what he used in his nested fixed point estimator).

While I don’t use the spectral methods SQUAREM or DF-SANE here, I have found them to be superior in a few contexts. First, with large datasets calculating the Jacobian for Newton’s method can be very expensive. SQUAREM and DF-SANE approximate the Jacobian with a diagonal matrix and so are much less computationally burdensome. You are then trading off between computational complexity and rates of convergence (Newton’s method converges quadratically while the other two converge superlinearly), often in favor of the spectral methods. Second, when using sparse grids for numerical integration, you potentially have to deal with negative market shares. This is becomes an issue for the contraction mapping, which we need to combine with Newton’s method for stability. The spectral methods are far more stable than Newton’s method, and therefore are more robust when using sparse grids. I have not found a situation in which the contraction mapping by itself was superior to any of these methods.

function gen_mm2(m::ModelData) α = m.alpha v_ik = m.v_ik J_f = m.J_f midx = m.midx for i = 1:length(m.p_ijt) m.p_ijt[i] = 0.0 end m.delta .= exp.(m.delta) step = Int(size(m.v_ik,1)/size(m.midx,1)) # iterate over markets for mi = 1:size(m.midx,1) muidx1 = step*(mi - 1) + 1; muidx2 = step*mi eps = 1.0 count = 1 r1 = midx[mi,1]; r2 = midx[mi,2]; # precalculate random coeffs mui = zeros(r2-r1+1, size(m.p_ijt,2)) # Allocate Array randcoeff = (m.sigma' .* m.v_ik[muidx1:muidx2,:]) yp = m.y_it[r1:r2,muidx1:muidx2].*m.price[r1:r2] # calculate random coefficient gemm!('N', 'T', 1.0, m.X[r1:r2,:], randcoeff, 0.0, mui) mui .-= α*yp mui .= exp.(mui); marketmat = ones(r2-r1+1, r2-r1+1) # contraction mapping/fixed point while (eps > 1e-14)*(200 > count) numer = mui.*m.delta[r1:r2]; s0 = 1.0./(1.0 + sum(numer,1)) m.p_ijt[r1:r2,:] .= numer .* s0; m.s_jt[r1:r2] .= m.p_ijt[r1:r2,:] * m.unobs_weight[muidx1:muidx2] # Finding the Fixed Point if eps > 0.1 # Nevo Contraction m.delta[r1:r2] .= m.delta[r1:r2] .* m.act_s[r1:r2]./m.s_jt[r1:r2] eps = maximum(abs.((m.s_jt[r1:r2]./m.act_s[r1:r2]) - 1)) else # Newton Root-finding m.delta[r1:r2] .= log.(m.delta[r1:r2]) tmp2 = m.p_ijt[r1:r2,:] .* m.unobs_weight[muidx1:muidx2]' Jf = Base.LinAlg.I - (tmp2*m.p_ijt[r1:r2,:]')./m.s_jt[r1:r2,:] diff = log.(m.s_jt[r1:r2]) .- log.(m.act_s[r1:r2]) stp = -Jf\diff m.delta[r1:r2,:] .= m.delta[r1:r2,:] .+ stp eps = maximum(abs.(stp)) m.delta[r1:r2] .= exp.(m.delta[r1:r2]) end count += 1 end end m.delta .= log.(m.delta) # calculate mc for j = 1:size(J_f, 1) idx1 = J_f[j,1]; idx2 = J_f[j,2]; muidx = searchsortedlast(m.midx[:,1], idx1 ) muidx_1 = (muidx-1)*step + 1 muidx_2 = (muidx)*step rp = m.p_ijt[idx1:idx2,:] wα = m.unobs_weight[muidx_1:muidx_2] .* m.y_it[idx1,muidx_1:muidx_2] * α nfirms = idx2 - idx1 + 1; tmp2 = rp .* wα' tmpvec = diagm(rp * wα); gemm!('N', 'T', -1.0, m.p_ijt[idx1:idx2,:], tmp2 , 1.0, tmpvec) b = tmpvec\m.s_jt[idx1:idx2] m.mc[idx1:idx2] .= m.price[idx1:idx2] .- b end for i = 1:length(m.mc) m.mc[i] = ifelse(m.mc[i] < 0, 0.001, m.mc[i]) end return nothing end

###### Estimation

The rest of the code is pretty much the same as before. The only real change I made was to use an initial approximation of the optimal GMM weighting matrix rather than the identity matrix in the first step of my estimation routine. This appears to better approximate BLP’s method.

d_share = convert(Array{Float64,1}, regSet[:share]) m = ModelData( price[:], X, α, βs, γs, σs, deepcopy(delta_0), deepcopy(delta_0), markets[:], midx, J_f, incy, v_ik1, unobs_weight1[:], zeros(size(X,1), size(v_ik1,1)/20), convert(Array{Float64,1}, deepcopy(d_share)), convert(Array{Float64,1}, deepcopy(d_share)), deepcopy(delta_0), zeros(size(X,1), 19)); @time gen_mm2(m) mean(m.delta)

0.256679 seconds (10.21 k allocations: 152.647 MiB, 14.09% gc time) -0.5655685632777482

function GMM(param, m, wght) m.alpha = abs(param[1]) m.sigma = abs.(param[2:end]) gen_mm2(m) Z_ = m.optinst y_ = vcat(m.delta, log.(m.mc)) X_ = [m.X zeros(m.W); zeros(m.X) m.W] zx = Z_'*X_ βols = (zx'wght*zx)\(zx'wght*Z_'y_) ϵ = y_ - X_*βols g = Z_'*ϵ/length(ϵ) m.beta = βols[1:size(m.X,2)] m.gamma = βols[size(m.X,2)+1:end] res = g'wght*g*length(g) println(m.alpha, m.sigma, res) return res end

m.W = W m.optinst = Z_; tp = [α, σs...]; wght = inv(Z_'Z_); y_ = vcat(m.delta, log.(m.mc)) X_ = [m.X zeros(m.W); zeros(m.X) m.W] zx = Z_'*X_ βols = (zx'wght*zx)\(zx'wght*Z_'y_) ϵ = y_ - X_*βols; mm = Z_ .* ϵ g = mean(mm, 1) wght = inv(cholfact(Hermitian(mm'mm/size(mm,1) .- g.*g'))); opt = Opt(:LN_COBYLA, 6) initial_step!(opt, abs.(tp)*0.1 ) ftol_rel!(opt,1e-4) xtol_rel!(opt,1e-4) maxeval!(opt, 250) min_objective!(opt, (x,y) -> GMM(x, m, wght)) @time minf_,minx_,ret_ = NLopt.optimize(opt, tp )

323.130396 seconds (3.91 M allocations: 130.944 GiB, 6.67% gc time) (4.330589470676522, [66.6355, 0.578001, 12.0295, 3.8127, 0.0768498, 5.59248], :FTOL_REACHED)

The code in my previous post took twice as long to run for the same number of draws, so this is a modest improvement. We can just repeat the above for the second stage to get your final estimate.

m.alpha = minx_[1] m.sigma = abs.(minx_[2:end]) gen_mm2(m); y_ = vcat(m.delta, log.(m.mc)) X_ = [m.X zeros(m.W); zeros(m.X) m.W] zx = Z_'*X_ βols = (zx'wght*zx)\(zx'wght*Z_'y_) ϵ = y_ - X_*βols; mm = Z_ .* ϵ g = mean(mm, 1) wght2 = inv(cholfact(Hermitian(mm'mm/size(mm,1) .- g.*g'))); opt = Opt(:LN_COBYLA, 6) initial_step!(opt, abs.(tp)*0.01 ) ftol_rel!(opt,1e-4) maxeval!(opt, 250) min_objective!(opt, (x,y) -> GMM(x, m, wght2)) @time minf_,minx_,ret_ = NLopt.optimize(opt, tp )

314.578066 seconds (3.66 M allocations: 121.229 GiB, 6.55% gc time) (4.786798337602201, [60.0746, 4.58298, 6.35377, 2.11136, -1.02675e-6, 4.18812], :FTOL_REACHED)

Note that the parameter estimates are pretty off from what BLP actually estimate. This is because 1,000 random draws does a poor job of approximating the market share integrals. I have previously discussed some methods to reduce simulation error, such as Halton draws and quadrature, and will provide some notes on importance sampling in a future post.

###### References

- Berry, Steven, James Levinsohn, and Ariel Pakes. “Automobile prices in market equilibrium.” Econometrica: Journal of the Econometric Society (1995): 841-890.
- Reynaerts, Jo, R. Varadha, and John C. Nash. “Enhencing the Convergence Properties of the BLP (1995) Contraction Mapping.” (2012).
- Skrainka, Benjamin S., and Kenneth L. Judd. “High performance quadrature rules: How numerical integration affects a popular model of product differentiation.” (2011).