In this session, I am going to cover demand estimation.

• Compute equilibrium outcomes with RCL demand
• Simulate market-level data
  • Extremely similar to the logit demand simulation
• Build the BLP estimator from Berry, Levinsohn, and Pakes (1995)

In this first part, we are going to assume that consumer \(i \in \lbrace1,...,I\rbrace\) utility from good \(j \in \lbrace1,...,J\rbrace\) in market \(t \in \lbrace1,...,T\rbrace\) takes the form

\[ u_{ijt} = \boldsymbol x_{jt} \boldsymbol \beta_{it} - \alpha p_{jt} + \xi_{jt} + \epsilon_{ijt} \]

where

• \(\xi_{jt}\) is type-1 extreme value distributed
• \(\boldsymbol \beta_{it}\): has dimension \(K\) \[\beta_{it}^k = \beta_0^k + \sigma_k \zeta_{it}^k\]
  • \(\beta_0^k\): fixed taste for characteristic \(k\) (the usual \(\beta\))
  • \(\zeta_{it}^k\): random taste, i.i.d. across consumers and markets \(t\)

We have \(J\) firms and each product has \(K\) characteristics.

i = 100;                # Number of consumers
J = 10;                 # Number of firms
K = 2;                  # Product characteristics
T = 100;                # Number of markets
β = [.5, 2, -1];        # Preferences
varζ = 5;               # Variance of the random taste
rangeJ = [2, 6];        # Min and max firms per market
varX = 1;               # Variance of X
varξ = 2;               # Variance of xi

Demand is the main difference w.r.t. the logit model. Now we have individual shocks \(\zeta\) we have to integrate over.

function demand(p::Vector, X::Matrix, β::Vector, ξ::Matrix, ζ::Matrix)::Tuple{Vector, Number}
    """Compute demand"""
    δ = [X p] * (β .+ ζ)                    # Mean value
    δ0 = zeros(1, size(ζ, 2))               # Mean value of the outside option
    u = [δ; δ0] + ξ                         # Utility
    e = exp.(u)                             # Take exponential
    q = mean(e ./ sum(e, dims=1), dims=2)   # Compute demand
    return q[1:end-1], q[end]
end;

Computing profits is instead exactly the same as before. We just have to save the shocks \(\zeta\) to be sure demand is stable.

function profits(p::Vector, c::Vector, X::Matrix, β::Vector, ξ::Matrix, ζ::Matrix)::Vector
    """Compute profits"""
    q, _ = demand(p, X, β, ξ, ζ)            # Compute demand
    pr = (p - c) .* q                       # Compute profits
    return pr
end;

function profits_j(pj::Number, j::Int, p::Vector, c::Vector, X::Matrix, β::Vector, ξ::Matrix, ζ::Matrix)::Number
    """Compute profits of firm j"""
    p[j] = pj                               # Insert price of firm j
    pr = profits(p, c, X, β, ξ, ζ)          # Compute profits
    return pr[j]
end;

We can now compute the equilibrium for a specific market, as before.

function equilibrium(c::Vector, X::Matrix, β::Vector, ξ::Matrix, ζ::Matrix)::Vector
    """Compute equilibrium prices and profits"""
    p = 2 .* c;
    dist = 1;
    iter = 0;

    # Iterate until convergence
    while (dist > 1e-8) && (iter<1000)

        # Compute best reply for each firm
        p_old = copy(p);
        for j=1:length(p)
            obj_fun(pj) = - profits_j(pj[1], j, p, c, X, β, ξ, ζ);
            optimize(x -> obj_fun(x), [1.0], LBFGS());
        end

        # Update distance
        dist = max(abs.(p - p_old)...);
        iter += 1;
    end
    return p
end;

We are now ready to simulate the data, i.e. equilibrium outcomes across different markets. We first draw all the variables.

function draw_data(I::Int, J::Int, K::Int, rangeJ::Vector, varζ::Number, varX::Number, varξ::Number)::Tuple
    """Draw data for one market"""
    J_ = rand(rangeJ[1]:rangeJ[2])              # Number of firms (products)
    X_ = rand(Exponential(varX), J_, K)         # Product characteristics
    ξ_ = rand(Normal(0, varξ), J_+1, I)         # Product-level utility shocks
    # Consumer-product-level preference shocks
    ζ_ = [rand(Normal(0,1), 1, I) * varζ; zeros(K,I)]
    w_ = rand(Uniform(0, 1), J_)                # Cost shifters
    ω_ = rand(Uniform(0, 1), J_)                # Cost shocks
    c_ = w_ + ω_                                # Cost
    j_ = sort(sample(1:J, J_, replace=false))   # Subset of firms
    return X_, ξ_, ζ_, w_, c_, j_
end;

Then we simulate the data for one market.

function compute_mkt_eq(I::Int, J::Int, β::Vector, rangeJ::Vector, varζ::Number, varX::Number, varξ::Number)::DataFrame
    """Compute equilibrium one market"""

    # Initialize variables
    K = size(β, 1) - 1
    X_, ξ_, ζ_, w_, c_, j_ = draw_data(I, J, K, rangeJ, varζ, varX, varξ)

    # Compute equilibrium
    p_ = equilibrium(c_, X_, β, ξ_, ζ_)    # Equilibrium prices
    q_, q0 = demand(p_, X_, β, ξ_, ζ_)     # Demand with shocks
    pr_ = (p_ - c_) .* q_                       # Profits

    # Save to data
    q0_ = ones(length(j_)) .* q0
    df = DataFrame(j=j_, w=w_, p=p_, q=q_, q0=q0_, pr=pr_)
    for k=1:K
      df[!,"x$k"] = X_[:,k]
      df[!,"z$k"] = sum(X_[:,k]) .- X_[:,k]
    end
    return df
end;

We repeat for \(T\) markets.

function simulate_data(I::Int, J::Int, β::Vector, T::Int, rangeJ::Vector, varζ::Number, varX::Number, varξ::Number)
    """Simulate full dataset"""
    df = compute_mkt_eq(I, J, β, rangeJ, varζ, varX, varξ)
    df[!, "t"] = ones(nrow(df)) * 1
    for t=2:T
        df_temp = compute_mkt_eq(I, J, β, rangeJ, varζ, varX, varξ)
        df_temp[!, "t"] = ones(nrow(df_temp)) * t
        append!(df, df_temp)
    end
    CSV.write("../data/blp.csv", df)
    return df
end;

Now let’s run the code

# Simulate
df = simulate_data(i, J, β, T, rangeJ, varζ, varX, varξ);

What does the data look like? Let’s switch to R!

# Read data
df = fread("../data/blp.csv")
kable(df[1:6,], digits=4)

The BLP estimation procedure

First, we need to compute the shares implied by aspecific vector of \(\delta\)s

function implied_shares(Xt_::Matrix, ζt_::Matrix, δt_::Vector, δ0::Matrix)::Vector
    """Compute shares implied by deltas and shocks"""
    u = [δt_ .+ (Xt_ * ζt_); δ0]                  # Utility
    e = exp.(u)                                 # Take exponential
    q = mean(e ./ sum(e, dims=1), dims=2)       # Compute demand
    return q[1:end-1]
end;

We can now compute the inner loop and invert the demand function: from shares \(q\) to \(\delta\)s

function inner_loop(qt_::Vector, Xt_::Matrix, ζt_::Matrix)::Vector
    """Solve the inner loop: compute delta, given the shares"""
    δt_ = ones(size(qt_))
    δ0 = zeros(1, size(ζt_, 2))
    dist = 1

    # Iterate until convergence
    while (dist > 1e-8)
        q = implied_shares(Xt_, ζt_, δt_, δ0)
        δt2_ = δt_ + log.(qt_) - log.(q)
        dist = max(abs.(δt2_ - δt_)...)
        δt_ = δt2_
    end
    return δt_
end;

We can now repeat the inversion for every market and get the vector of mean utilities \(\delta\)s from the observed market shares \(q\).

function compute_delta(q_::Vector, X_::Matrix, ζ_::Matrix, T::Vector)::Vector
    """Compute residuals"""
    δ_ = zeros(size(T))

    # Loop over each market
    for t in unique(T)
        qt_ = q_[T.==t]                             # Quantity in market t
        Xt_ = X_[T.==t,:]                           # Characteristics in mkt t
        δ_[T.==t] = inner_loop(qt_, Xt_, ζ_)        # Solve inner loop
    end
    return δ_
end;

Now that we have \(\delta\), it is pretty straightforward to compute \(\xi\). We just need to perform a linear regression (with instruments) of mean utilities \(\delta\) on prices \(p\) and product characteristics \(X\) and compute the residuals \(\xi\).

function compute_xi(X_::Matrix, IV_::Matrix, δ_::Vector)::Tuple
    """Compute residual, given delta (IV)"""
    β_ = inv(IV_' * X_) * (IV_' * δ_)           # Compute coefficients (IV)
    ξ_ = δ_ - X_ * β_                           # Compute errors
    return ξ_, β_
end;

We now have all the ingredients to set up the GMM objective function.

function GMM(varζ_::Number)::Tuple
    """Compute GMM objective function"""
    δ_ = compute_delta(q_, X_, ζ_ * varζ_, T)   # Compute deltas
    ξ_, β_ = compute_xi(X_, IV_, δ_)            # Compute residuals
    gmm = ξ_' * Z_ * Z_' * ξ_ / length(ξ_)^2    # Compute ortogonality condition
    return gmm, β_
end;

First, we need to set up our objects

# Retrieve data
T = Int.(df.t)
X_ = [df.x1 df.x2 df.p]
q_ = df.q
q0_ = df.q0
IV_ = [df.x1 df.x2 df.w]
Z_ = [df.x1 df.x2 df.z1 df.z2]

What would a logit regression estimate?

# Compute logit estimate
y = log.(df.q) - log.(df.q0);
β_logit = inv(IV_' * X_) * (IV_' * y);
print("Estimated logit coefficients: $β_logit")

## Estimated logit coefficients: [1.6707629684340313, 1.090127651610374, -0.8676967920899061]

We can now run the BLP machinery

# Draw shocks (less)
ζ_ = [rand(Normal(0,1), 1, i); zeros(K, i)];

# Minimize GMM objective function
varζ_ = optimize(x -> GMM(x[1])[1], [2.0], LBFGS()).minimizer[1];
β_blp = GMM(varζ_)[2];
print("Estimated BLP coefficients: $β_blp")

## Estimated BLP coefficients: [0.5404717706555511, 1.1612479679061647, -0.5876407838141476]