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

• Compute equilibrium outcomes with Logit demand
• Simulate a dataset
• Estimate Logit demand
• Compare different instruments
• Include supply

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\) has dimension \(K\)
  • i.e. goods have \(K\) characteristics

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

J = 3;                            # 3 firms == products
K = 2;                            # 2 product characteristics
c = rand(Uniform(0, 1), J);       # Random uniform marginal costs
ξ = rand(Normal(0, 1), J+1);      # Random normal individual shocks
X = rand(Exponential(1), J, K);   # Random exponential product characteristics
β = [.5, 2, -1];                  # Preferences (last one is for prices, i.e. alpha)

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

We can try with an example.

p = 2 .* c;
demand(p, X, β, ξ)

## ([0.5491221990381039, 0.24359787592994134, 0.1810840379517347], 0.026195887080220043)

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

We can try with an example.

profits(p, c, X, β, ξ)

## 3-element Array{Float64,1}:
##  0.16600481498889777
##  0.0889567306013599
##  0.12280474674253772

We first code the best reply of firm \(j\)

function profits_j(pj::Number, j::Int, p::Vector, c::Vector, X::Matrix, β::Vector, ξ::Vector)::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;

Let’s test it.

j = 1;
obj_fun(pj) = - profits_j(pj[1], j, copy(p), c, X, β, ξ);
pj = optimize(x -> obj_fun(x), [1.0], LBFGS()).minimizer[1]

## 1.706820256123432

What are the implied profits now?

print("Profits old: ",  round.(profits(p, c, X, β, ξ), digits=4))

## Profits old: [0.166, 0.089, 0.1228]

p_new = copy(p);
p_new[j] = pj;
print("Profits new: ",  round.(profits(p_new, c, X, β, ξ), digits=4))

## Profits new: [0.4045, 0.1405, 0.1939]

Indeed firm 1 has increased its profits.

We can now compute equilibrium prices

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

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

        # Compute best reply for each firm
        p1 = 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()).minimizer[1];
        end

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

Let’s test it

# Compute equilibrium prices
p_eq = equilibrium(c, X, β, ξ);
print("Equilibrium prices: ",  round.(p_eq, digits=4))

## Equilibrium prices: [1.9722, 1.7173, 2.0356]

# And profits
pi_eq = profits(p_eq, c, X, β, ξ);
print("Equilibrium profits: ",  round.(pi_eq, digits=4))

## Equilibrium profits: [0.6699, 0.3521, 0.3574]

As expected the prices of the first 2 firms are lower and their profits are higher.

Let’s generate our Data Generating Process (DGP).

• \(\boldsymbol x \sim exp(V_{x})\)
• \(\xi \sim N(0, V_{\xi})\)
• \(w \sim N(0, 1)\)
• \(\omega \sim N(0, 1)\)

function draw_data(J::Int, K::Int, rangeJ::Vector, 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)            # Product-level utility shocks
    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;

We first compute the equilibrium in one market.

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

    # Initialize variables
    K = size(β, 1) - 1
    X_, ξ_, w_, c_, j_ = draw_data(J, K, rangeJ, 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 can now write the code to simulate the whole dataset.

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

We generate the dataset by simulating many markets that differ by

• number of firms (and their identity)
• their marginal costs
• their product characteristics

# Set parameters
J = 10;                 # Number of firms
K = 2;                  # Product caracteristics
T = 500;                # Markets
β = [.5, 2, -1];        # Preferences
rangeJ = [2, 6];        # Min and max firms per market
varX = 1;               # Variance of X
varξ = 2;               # Variance of xi

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

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

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

First we need to compute the dependent variable

df$y = log(df$q) - log(df$q0)

Now we can estimate the logit model. The true values are \(alpha=1\).

ols <- lm(y ~ x1 + x2 + p, data=df)
kable(tidy(ols), digits=4)

The estimate of \(\alpha = 1\) is biased (positive and significant) since \(p\) is endogenous. We need instruments.

First set of instruments: cost shifters.

fm_costiv <- ivreg(y ~ x1 + x2 + p | x1 + x2 + w, data=df)
kable(tidy(fm_costiv), digits=4)

Now the estimate of \(\alpha\) is negative and significant.

Second set of instruments: product characteristics of other firms in the same market.

fm_blpiv <- ivreg(y ~ x1 + x2 + p | x1 + x2 + z1 + z2, data=df)
kable(tidy(fm_blpiv), digits=4)

Also the BLP instruments deliver an estimate of \(\alpha\) is negative and significant.