Coding: Logit Demand

Intro

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

Model

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

Setup

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)

Code Demand

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.4120077746005573, 0.26650568009936, 0.24027826270165709], 0.08120828259842561)

Code Supply

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.20289186252172428
##  0.1596422305025479
##  0.1270874470740512

Code Best Reply

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.8019637881982011

What are the implied profits now?

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

## Profits old: [0.2029, 0.1596, 0.1271]

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

## Profits new: [0.3095, 0.2073, 0.1651]

Indeed firm 1 has increased its profits.

Code Equilibrium

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;

Code Equilibrium

Let’s test it

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

## Equilibrium prices: [1.9764, 1.9602, 1.8366]

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

## Equilibrium profits: [0.484, 0.3612, 0.3077]

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

DGP

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;

Equilibrium

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;

Simulate Dataset

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;

Simulate Dataset (2)

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ξ);

The Data

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

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

Estimation

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.

IV 1: Cost Shifters

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.

IV 2: BLP Instruments

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.

Appendix

References