A Cornerstone in Structural Estimation: BLP's method ¶
Alfred Galichon (NYU & Sciences Po) and Ugo Arena (Sciences Po) ¶
'math+econ+code' masterclass series ¶
With python code examples ¶
© 2018-2026 by Alfred Galichon. Past and present support from NSF grant DMS-1716489, ERC grant CoG-866274 are acknowledged, as well as inputs from contributors listed here.
If you reuse material from this masterclass, please cite as:
Alfred Galichon, 'math+econ+code' masterclass. https://github.com/math-econ-code/mec_notebooks
References¶
Lin, Ranie. "Berry-Levinsohn-and-Pakes-1995-Replication." GitHub. https://github.com/ranielin/Berry-Levinsohn-and-Pakes-1995-Replication.
Galichon, Alfred. Discrete Choice Models: Mathematical Methods, Data Science, and Econometrics. Princeton Univeristy Press.
Galichon, Alfred. "math-econ-code." GitHub. https://github.com/math-econ-code.
Berry, Steven, James Levinsohn, and Ariel Pakes. "Automobile Prices in Market Equilibrium." Econometrica 63, no. 4 (1995): 841–90. https://doi.org/10.2307/2171802.
Conlon, Christopher and Jeff Gortmaker (2020). Best practices for differentiated products demand estimation with PyBLP. In: The RAND Journal of Economics 51.4, pp. 1108–1161.
The aim of this notebook is to replicate BLP's results, from Berry, Levinsohn, and Pakes, Automobile Prices in Market Equilibrium (1995)
import pandas as pd, numpy as np
from scipy.optimize import minimize
import scipy, datetime
from sklearn.linear_model import LinearRegression
from scipy.linalg import block_diag
!pip install pyblp
import pyblp
!pip install mec
from mec.data import load_blp_car_data
from mec.blp import create_blp_instruments, organize_markets, collapse_markets,build_nus,build_dnudps
import pandas_datareader.data as web
import statsmodels.api as sm
from statsmodels.sandbox.regression.gmm import IV2SLS
pyblp.options.verbose = False
Collecting pyblp Downloading pyblp-1.2.0-py3-none-any.whl.metadata (8.5 kB) Requirement already satisfied: numpy>=1.17.0 in /usr/local/lib/python3.12/dist-packages (from pyblp) (2.0.2) Requirement already satisfied: patsy>=0.5.1 in /usr/local/lib/python3.12/dist-packages (from pyblp) (1.0.2) Collecting pyhdfe>=0.1.0 (from pyblp) Downloading pyhdfe-0.2.0-py3-none-any.whl.metadata (4.0 kB) Requirement already satisfied: scipy>=1.2.0 in /usr/local/lib/python3.12/dist-packages (from pyblp) (1.16.3) Requirement already satisfied: sympy>=1.1.0 in /usr/local/lib/python3.12/dist-packages (from pyblp) (1.14.0) Requirement already satisfied: mpmath<1.4,>=1.1.0 in /usr/local/lib/python3.12/dist-packages (from sympy>=1.1.0->pyblp) (1.3.0) Downloading pyblp-1.2.0-py3-none-any.whl (2.1 MB) ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 2.1/2.1 MB 12.6 MB/s eta 0:00:00 Downloading pyhdfe-0.2.0-py3-none-any.whl (19 kB) Installing collected packages: pyhdfe, pyblp Successfully installed pyblp-1.2.0 pyhdfe-0.2.0 Collecting mec Downloading mec-0.231-py3-none-any.whl.metadata (392 bytes) Collecting gurobipy (from mec) Downloading gurobipy-13.0.1-cp312-cp312-manylinux2014_x86_64.manylinux_2_17_x86_64.whl.metadata (16 kB) Collecting PyGithub (from mec) Downloading pygithub-2.9.1-py3-none-any.whl.metadata (3.9 kB) Requirement already satisfied: python-dotenv in /usr/local/lib/python3.12/dist-packages (from mec) (1.2.2) Collecting pynacl>=1.4.0 (from PyGithub->mec) Downloading pynacl-1.6.2-cp38-abi3-manylinux_2_34_x86_64.whl.metadata (10.0 kB) Requirement already satisfied: requests>=2.14.0 in /usr/local/lib/python3.12/dist-packages (from PyGithub->mec) (2.32.4) Requirement already satisfied: pyjwt>=2.4.0 in /usr/local/lib/python3.12/dist-packages (from pyjwt[crypto]>=2.4.0->PyGithub->mec) (2.12.1) Requirement already satisfied: typing-extensions>=4.5.0 in /usr/local/lib/python3.12/dist-packages (from PyGithub->mec) (4.15.0) Requirement already satisfied: urllib3>=1.26.0 in /usr/local/lib/python3.12/dist-packages (from PyGithub->mec) (2.5.0) Requirement already satisfied: cryptography>=3.4.0 in /usr/local/lib/python3.12/dist-packages (from pyjwt[crypto]>=2.4.0->PyGithub->mec) (43.0.3) Requirement already satisfied: cffi>=2.0.0 in /usr/local/lib/python3.12/dist-packages (from pynacl>=1.4.0->PyGithub->mec) (2.0.0) Requirement already satisfied: charset_normalizer<4,>=2 in /usr/local/lib/python3.12/dist-packages (from requests>=2.14.0->PyGithub->mec) (3.4.7) Requirement already satisfied: idna<4,>=2.5 in /usr/local/lib/python3.12/dist-packages (from requests>=2.14.0->PyGithub->mec) (3.13) Requirement already satisfied: certifi>=2017.4.17 in /usr/local/lib/python3.12/dist-packages (from requests>=2.14.0->PyGithub->mec) (2026.4.22) Requirement already satisfied: pycparser in /usr/local/lib/python3.12/dist-packages (from cffi>=2.0.0->pynacl>=1.4.0->PyGithub->mec) (3.0) Downloading mec-0.231-py3-none-any.whl (3.9 MB) ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 3.9/3.9 MB 12.9 MB/s eta 0:00:00 Downloading gurobipy-13.0.1-cp312-cp312-manylinux2014_x86_64.manylinux_2_17_x86_64.whl (14.8 MB) ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 14.8/14.8 MB 61.3 MB/s eta 0:00:00 Downloading pygithub-2.9.1-py3-none-any.whl (449 kB) ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 449.7/449.7 kB 29.3 MB/s eta 0:00:00 Downloading pynacl-1.6.2-cp38-abi3-manylinux_2_34_x86_64.whl (1.4 MB) ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 1.4/1.4 MB 50.4 MB/s eta 0:00:00 Installing collected packages: gurobipy, pynacl, PyGithub, mec Successfully installed PyGithub-2.9.1 gurobipy-13.0.1 mec-0.231 pynacl-1.6.2
Loading our Data¶
We import BLP's original dataset using the math+econ+code package (specifically the one from Voluntary Export Restraints on Automobiles: Evaluating a Trade Policy, BLP 1999), which is identical to the one used in their 1995 paper. We extract the number of US households (in thousands) from 1971 to 1991 using FRED. Because BLP used another source for this demographic data, there will be some discrepancies in the descriptive statistics.
prod, _ = load_blp_car_data()
print(prod.head())
name car_ids market_ids cy dr at ps air drv prices ... \ 0 AMGREM 129 1971 6 2 0 0 0 -9 4.935802 ... 1 AMHORN 130 1971 6 4 0 0 0 -9 5.516049 ... 2 AMJAVL 132 1971 6 2 0 0 0 -9 7.108642 ... 3 AMMATA 134 1971 6 4 0 0 0 -9 6.839506 ... 4 AMAMBS 136 1971 6 4 1 0 0 -9 8.928395 ... supply_instruments2 supply_instruments3 supply_instruments4 \ 0 0.0 1.705933 1.595656 1 0.0 1.680910 1.490295 2 0.0 1.801067 1.357703 3 0.0 1.818061 1.261347 4 0.0 1.933210 1.237365 supply_instruments5 supply_instruments6 supply_instruments7 \ 0 87.0 -61.959985 0.0 1 87.0 -61.959985 0.0 2 87.0 -61.959985 0.0 3 87.0 -61.959985 0.0 4 87.0 -61.959985 0.0 supply_instruments8 supply_instruments9 supply_instruments10 \ 0 46.060389 29.786989 0 1 46.060389 29.786989 0 2 46.060389 29.786989 0 3 46.060389 29.786989 0 4 46.060389 29.786989 0 supply_instruments11 0 1.888146 1 1.935989 2 1.716799 3 1.687871 4 1.504286 [5 rows x 52 columns]
start = datetime.datetime(1971, 1, 1)
end = datetime.datetime(1990, 1, 1)
households_df = web.DataReader('TTLHH', 'fred', start, end)
households_df['Year'] = households_df.index.year
households_1000s = pd.Series(households_df.TTLHH.values, index=households_df.Year).to_dict()
Descriptive Statistics (Replication of Table I & II)¶
Description of the BLP dataset
The original data on product characteristics were obtained from annual issues of the Automotive News Market Data Book. The included product characteristics are the number of cylinders, doors, weight, engine displacement, horsepower, length, width, wheelbase, MPG rating, and dummy variables for standard equipment such as air conditioning. The price used is the base model's list retail price in thousands of 1983 dollars. The sales variable corresponds to U.S. sales in thousands by nameplate, and the characteristics reflect the base model for a given nameplate.
The dataset contains all models marketed from 1971 to 1990, excluding models with insignificant market shares. The authors treat each model/year combination as an observation, yielding a sample size of 2,217. Two observations in adjacent years represent the same model if they meet two criteria: (i) they share the same nameplate, and (ii) their horsepower, width, length, and wheelbase do not change by more than 10%. This classification results in a total of 997 distinct models.
In the dataset, a car's price is endogenous because it is correlated with unobserved characteristics ($\xi$). For example, a car might have a great reputation for reliability or a beautiful design that isn't captured by the data's basic metrics. If we don't correct for this, our price sensitivity estimate ($\alpha$) will be severely biased. To fix this, the authors constructed instruments directly from the dataset using the physical characteristics of the cars. The dataset constructs three categories of instruments:
- Own Characteristics: Physical traits of the car itself.
- Same-Firm Characteristics : The sum of the characteristics of all other cars produced by the same automaker in that market
- Rival-Firm Characteristics (Market Competition) : The sum of the characteristics of all cars produced by competing automakers in that market.
Additional data was gathered and transformed; for instance, miles per dollar (MP$) was calculated by dividing MPG by the real price of unleaded gasoline (sourced from Business Statistics). The market size ($M$) was defined as the number of US households in thousands, originally sourced from the Statistical Abstract of the US (though we draw this from FRED). Finally, we provide some summary descriptive statistics for the variables used below.
Table I : Descriptive statistics¶
Table I provides a year-by-year summary of the U.S. automobile market's overarching trends from 1971 to 1990.
It displays the total number of distinct car models available each year alongside their average sales quantities and real prices.
Additionally, it tracks the evolution of sales-weighted vehicle characteristics, such as horsepower-to-weight ratio, size, and fuel efficiency.
Finally, it highlights the shifting competitive landscape by showing the market shares of Domestic, Japanese, and European manufacturers.
if 'region' in prod.columns:
region_str = prod['region'].astype(str).str.upper()
prod['Domestic'] = region_str.str.contains('US').astype(int)
prod['Japan'] = region_str.str.contains('JP').astype(int)
prod['European'] = region_str.str.contains('EU').astype(int)
# We reconstruct Quantity (Sales in 1000s)
prod['market_size'] = prod['market_ids'].map(households_1000s)
prod['Quantity'] = prod['shares'] * prod['market_size']
# We define columns for the sales-weighted mean
cols_to_weight = ['prices', 'Domestic', 'Japan', 'European', 'hpwt', 'space', 'air', 'mpg', 'mpd']
def sales_weighted_mean(df, value_col):
return np.average(df[value_col], weights=df['Quantity'])
# We build the table row by year
table_rows = []
years = sorted(prod['market_ids'].unique())
for year in years:
group = prod[prod['market_ids'] == year]
row = {
'Year': str(year),
'No. of Models': len(group),
'Quantity': group['Quantity'].mean()
}
for col in cols_to_weight:
row[col] = sales_weighted_mean(group, col)
table_rows.append(row)
row_all = {
'Year': 'All',
'No. of Models': len(prod),
'Quantity': prod['Quantity'].mean()
}
for col in cols_to_weight:
row_all[col] = sales_weighted_mean(prod, col)
table_rows.append(row_all)
We format and display the table, which we compare with the table in the original paper.
table_1_df = pd.DataFrame(table_rows)
rename_map = {
'prices': 'Price',
'hpwt': 'HP/Wt',
'space': 'Size',
'air': 'Air',
'mpg': 'MPG',
'mpd': 'MP$'
}
table_1_df.rename(columns=rename_map, inplace=True)
format_cols = [col for col in table_1_df.columns if col not in ['Year', 'No. of Models']]
formatters = {col: '{:.3f}'.format for col in format_cols}
pd.set_option('display.max_rows', None)
print("- Table I: Descriptive Statistics -")
print(table_1_df.to_string(formatters=formatters, index=False))
- Table I: Descriptive Statistics - Year No. of Models Quantity Price Domestic Japan European HP/Wt Size Air MPG MP$ 1971 92 84.418 7.868 0.866 0.057 0.077 0.490 1.496 0.000 1.662 1.849 1972 89 96.347 7.979 0.892 0.042 0.066 0.391 1.510 0.014 1.619 1.875 1973 86 90.650 7.535 0.932 0.040 0.028 0.364 1.529 0.022 1.589 1.818 1974 72 102.330 7.506 0.887 0.050 0.064 0.347 1.510 0.026 1.567 1.452 1975 93 82.742 7.821 0.853 0.083 0.064 0.337 1.479 0.054 1.584 1.503 1976 99 91.776 7.787 0.876 0.081 0.043 0.338 1.508 0.059 1.759 1.696 1977 95 95.300 7.651 0.837 0.112 0.051 0.340 1.467 0.032 1.947 1.835 1978 95 97.773 7.645 0.855 0.107 0.039 0.346 1.405 0.034 1.982 1.929 1979 102 79.212 7.599 0.803 0.158 0.038 0.348 1.343 0.047 2.061 1.657 1980 103 70.184 7.718 0.773 0.191 0.036 0.350 1.296 0.078 2.215 1.466 1981 116 61.170 8.349 0.741 0.213 0.046 0.349 1.286 0.094 2.363 1.559 1982 110 61.604 8.831 0.714 0.235 0.051 0.347 1.277 0.134 2.440 1.817 1983 115 66.786 8.821 0.734 0.215 0.051 0.351 1.276 0.126 2.601 2.087 1984 113 84.449 8.870 0.783 0.179 0.038 0.361 1.293 0.129 2.469 2.117 1985 136 76.669 8.938 0.761 0.191 0.048 0.372 1.265 0.140 2.261 2.024 1986 130 82.800 9.382 0.733 0.216 0.050 0.379 1.249 0.176 2.416 2.856 1987 143 66.446 9.965 0.702 0.245 0.052 0.395 1.246 0.229 2.327 2.789 1988 150 65.845 10.069 0.717 0.237 0.045 0.396 1.251 0.237 2.334 2.919 1989 147 62.565 10.321 0.690 0.261 0.049 0.406 1.259 0.289 2.310 2.806 1990 131 65.698 10.337 0.682 0.276 0.043 0.419 1.270 0.308 2.270 2.852 All 2217 77.314 8.609 0.789 0.162 0.049 0.372 1.356 0.116 2.100 2.088
Table II : The range of continuous demand characteristics¶
Table II illustrates the overall distribution and spread of the continuous variables used to model consumer demand.
It reports the minimum (0th), 25th percentile, median (50th), 75th percentile, and maximum (100th) values for each characteristic.
These variables include critical decision factors like price, total sales volume, physical vehicle size, and miles per dollar.
The table identifies the specific car models that represent the absolute minimum and maximum extremes.
# We define variables matching Table II
continuous_vars = ['prices', 'Quantity', 'hpwt', 'space', 'mpg', 'mpd']
# The paper reports the 0th (Min), 25th, 50th, 75th, and 100th (Max) percentiles, we will do the same
quantiles = [0.0, 0.25, 0.50, 0.75, 1.0]
# We calculate the quantiles for each variable
table_2_data = {}
for var in continuous_vars:
table_2_data[var] = prod[var].quantile(quantiles).values
# We create the DataFrame and transpose it so Variables are rows
table_2_df = pd.DataFrame(table_2_data, index=['0', '25', '50', '75', '100']).T
# We rename the index to match the exact Variable names in the paper
rename_map = {
'prices': 'Price',
'Quantity': 'Sales',
'hpwt': 'HP/Wt',
'space': 'Size',
'mpg': 'MPG',
'mpd': 'MP$'
}
table_2_df.rename(index=rename_map, inplace=True)
for col in table_2_df.columns:
table_2_df[col] = table_2_df[col].apply(lambda x: f"{float(x):.3f}")
# We dynamically append car names for the Min (0) and Max (100) percentiles
if 'name' in prod.columns:
for var, row_name in zip(continuous_vars, table_2_df.index):
# Find the index of the min and max values
min_idx = prod[var].idxmin()
max_idx = prod[var].idxmax()
min_name = prod.loc[min_idx, 'name']
min_year = str(prod.loc[min_idx, 'market_ids'])[-2:]
max_name = prod.loc[max_idx, 'name']
max_year = str(prod.loc[max_idx, 'market_ids'])[-2:]
table_2_df.at[row_name, '0'] += f" ('{min_year} {min_name})"
table_2_df.at[row_name, '100'] += f" ('{max_year} {max_name})"
print("- TABLE II -")
print("THE RANGE OF CONTINUOUS DEMAND CHARACTERISTICS\n")
print(table_2_df.to_string(justify='left'))
- TABLE II -
THE RANGE OF CONTINUOUS DEMAND CHARACTERISTICS
0 25 50 75 100
Price 3.393 ('90 YGGVPL) 6.714 8.729 13.074 68.597 ('89 PS911C)
Sales 0.048 ('73 TY1600) 15.161 46.533 107.839 631.606 ('72 BKRIVE)
HP/Wt 0.170 ('80 PLGRNF) 0.337 0.375 0.428 0.948 ('89 PS911T)
Size 0.756 ('73 HDCIVI) 1.131 1.270 1.453 1.888 ('72 IMLEBA)
MPG 0.913 ('74 CDELDR) 1.700 2.000 2.500 5.300 ('89 GEOMET)
MP$ 0.846 ('74 CDELDR) 1.557 2.010 2.483 6.437 ('89 GEOMET)
The Logit and the IV Logit¶
Table III : Results with Logit Demand and Marginal Cost Pricing¶
This initial set of results is based on the simplest logit utility specification. This model is estimated via OLS by regressing the difference in log-market shares, $\ln(s_j) - \ln(s_0)$, on product characteristics and price.
Although one could argue for a more extensive set of attributes, the authors adopt a concise list for the sake of comparability.A key feature of this specification is the treatment of unobserved product attributes ($\xi_j$). The model's error term explicitly accounts for characteristics that are observed by consumers and influence their purchasing decisions but remain unobserved by the econometrician.
# We Prepare Dependent Variables
# Logit Demand Dependent Variable: ln(s_j) - ln(s_0)
prod['sum_shares'] = prod.groupby('market_ids')['shares'].transform('sum')
prod['s_0'] = 1 - prod['sum_shares']
prod['y_demand'] = np.log(prod['shares']) - np.log(prod['s_0'])
# Cost Function Dependent Variable: ln(price)
prod['y_supply'] = np.log(prod['prices'])
# We prepare Independent Variables
prod['const'] = 1
# Trend in BLP is in our case the two digits of the year
prod['trend'] = prod['market_ids'] - 1900
# We apply log transformations for the cost function covariates
prod['ln_hpwt'] = np.log(prod['hpwt'])
prod['ln_space'] = np.log(prod['space'])
prod['ln_mpg'] = np.log(prod['mpg'])
# We define covariate columns exactly as they appear in Table III
dem_cols = ['const', 'hpwt', 'air', 'mpd', 'space', 'prices']
sup_cols = ['const', 'ln_hpwt', 'air', 'ln_mpg', 'ln_space', 'trend']
X_dem = prod[dem_cols].values
y_dem = prod['y_demand'].values
X_sup = prod[sup_cols].values
y_sup = prod['y_supply'].values
# We prepare Instruments for IV
dem_inst_cols = [c for c in prod.columns if c.startswith('demand_instruments')]
sup_inst_cols = [c for c in prod.columns if c.startswith('supply_instruments')]
if len(dem_inst_cols) > 0 and len(sup_inst_cols) > 0:
Z_dem_cols = ['const', 'hpwt', 'air', 'mpd', 'space'] + dem_inst_cols
Z_sup_cols = ['const', 'ln_hpwt', 'air', 'ln_mpg', 'ln_space', 'trend'] + sup_inst_cols
Z_dem = prod[Z_dem_cols].values
Z_sup = prod[Z_sup_cols].values
else:
def build_blp_instruments(df, X_cols):
firm_sums, market_sums = [], []
for m in df['market_ids'].unique():
market_df = df[df['market_ids'] == m]
for f in market_df['firm_ids'].unique():
firm_df = market_df[market_df['firm_ids'] == f]
firm_sums.append(firm_df[X_cols].sum() - firm_df[X_cols])
other_firm_df = market_df[market_df['firm_ids'] != f]
market_sums.append(pd.DataFrame([other_firm_df[X_cols].sum()] * len(firm_df), index=firm_df.index))
return np.hstack([df[X_cols].values, pd.concat(firm_sums).sort_index().values, pd.concat(market_sums).sort_index().values])
Z_dem = build_blp_instruments(prod, ['const', 'hpwt', 'air', 'mpd', 'space'])
Z_sup = build_blp_instruments(prod, ['const', 'ln_hpwt', 'air', 'ln_mpg', 'ln_space', 'trend'])
Next, estimate paramters using Statsmodels.
# Demand (OLS)
res_dem_ols = sm.OLS(y_dem, X_dem).fit()
beta_dem_ols, se_dem_ols = res_dem_ols.params, res_dem_ols.bse
r2_dem_ols = res_dem_ols.rsquared
# Demand (IV/2SLS)
res_dem_iv = IV2SLS(endog=y_dem, exog=X_dem, instrument=Z_dem).fit()
beta_dem_iv, se_dem_iv = res_dem_iv.params, res_dem_iv.bse
# Supply / Marginal Cost (OLS)
res_sup_ols = sm.OLS(y_sup, X_sup).fit()
beta_sup_ols, se_sup_ols = res_sup_ols.params, res_sup_ols.bse
r2_sup_ols = res_sup_ols.rsquared
# ==
# We format and display Table III
# ==
param_names = ['Constant', 'HP/Weight', 'Air', 'MP$', 'Size', 'Price',
'ln(HP/Weight)', 'ln(MPG)', 'ln(Size)', 'Trend']
table_3 = pd.DataFrame(index=param_names, columns=['Demand (OLS)', 'Demand (IV)', 'ln(MC) (OLS)'])
table_3.fillna('', inplace=True)
def fmt(b, s):
return f"{b:8.3f}\n({s:6.3f})"
for i, var in enumerate(['Constant', 'HP/Weight', 'Air', 'MP$', 'Size', 'Price']):
table_3.at[var, 'Demand (OLS)'] = fmt(beta_dem_ols[i], se_dem_ols[i])
table_3.at[var, 'Demand (IV)'] = fmt(beta_dem_iv[i], se_dem_iv[i])
for i, var in enumerate(['Constant', 'ln(HP/Weight)', 'Air', 'ln(MPG)', 'ln(Size)', 'Trend']):
table_3.at[var, 'ln(MC) (OLS)'] = fmt(beta_sup_ols[i], se_sup_ols[i])
print(f"R^2 for Demand (OLS): {r2_dem_ols:.4f}")
print(f"R^2 for ln(Price) (OLS): {r2_sup_ols:.4f}\n")
print("--- TABLE III ---")
print("RESULTS WITH LOGIT DEMAND AND MARGINAL COST PRICING")
print(table_3)
R^2 for Demand (OLS): 0.3871
R^2 for ln(Price) (OLS): 0.6564
--- TABLE III ---
RESULTS WITH LOGIT DEMAND AND MARGINAL COST PRICING
Demand (OLS) Demand (IV) ln(MC) (OLS)
Constant -10.072\n( 0.253) -9.921\n( 0.262) 1.882\n( 0.119)
HP/Weight -0.124\n( 0.277) 1.179\n( 0.403)
Air -0.034\n( 0.073) 0.468\n( 0.133) 0.680\n( 0.019)
MP$ 0.265\n( 0.043) 0.175\n( 0.049)
Size 2.342\n( 0.125) 2.293\n( 0.129)
Price -0.089\n( 0.004) -0.134\n( 0.011)
ln(HP/Weight) 0.520\n( 0.035)
ln(MPG) -0.471\n( 0.049)
ln(Size) 0.125\n( 0.063)
Trend 0.013\n( 0.002)
The Demand Side: OLS vs. IV¶
In the Demand (OLS) column, the model produces counterintuitive results. The coefficient on Price is relatively small (-0.089). Worse, the coefficients for desirable attributes like HP/Weight (-0.124) and Air Conditioning (-0.034) are negative, falsely implying that consumers dislike those attributes. It is caused by an endogeneity problem. OLS cannot account for unobserved quality (like brand prestige, ride smoothness, or aesthetic design), biasing price towards zero and distorting attributes.
The Demand (IV) column uses instruments to isolate the variation in price that is strictly due to market competition, uncorrelated from unobserved quality. Once corrected, the Price coefficient becomes much more negative (-0.134), showing that consumers are correctly identified as being more price-sensitive. Furthermore, the coefficients for HP/Weight (1.179) and Air (0.468) flip to their expected positive signs, confirming that consumers do indeed value characteristics.
The Supply Side (Marginal Cost)¶
The ln(MC) (OLS) column represents a basic cost function.
The positive coefficients for HP/Weight (0.520), Air (0.680), and Size (0.125) make perfect economic sense: it costs manufacturers more money to build larger, heavier, more powerful, and better-equipped vehicles.
If we do not properly instrument for price, the demand model will generate biased consumer preferences and unrealistic substitution patterns.
Results from the Full Model¶
Table IV : Estimated Parameters of the Demand and Pricing Equations¶
Table IV presents the core structural parameter estimates from the full random coefficients logit model, covering both consumer demand and firm pricing equations.
The demand side reports the mean consumer preferences ($\beta$) for vehicle attributes alongside the standard deviations ($\sigma$) that capture unobserved consumer heterogeneity.
It also highlights the crucial price coefficient ($\alpha$), which interacts with consumer income to reflect realistic, income-based price sensitivities.
Finally, the supply side details the marginal cost parameters ($\gamma$), showing how features like vehicle size, air conditioning, and horsepower directly impact manufacturing costs.
Following the recommendations of Conlon and Gortmaker (2020), we manually construct the agent data using 2,000 Halton draws via PyBLP’s built-in tools.
# Get the Halton nodes from PyBLP
base_integration = pyblp.build_integration(
integration=pyblp.Integration('halton', size=2000),
dimensions=6
)
agent_dict = {
'weights': base_integration['weights'].flatten()
}
for i in range(6):
agent_dict[f'nodes{i}'] = base_integration['nodes'][:, i].flatten()
base_df = pd.DataFrame(agent_dict)
# Historical log-income means (approximate CPS data for 1971-1990)
D_mean = {
1971: 2.01156, 1972: 2.06526, 1973: 2.07843, 1974: 2.05775,
1975: 2.02915, 1976: 2.05346, 1977: 2.06745, 1978: 2.09805,
1979: 2.10404, 1980: 2.07208, 1981: 2.06019, 1982: 2.06561,
1983: 2.07672, 1984: 2.10437, 1985: 2.12608, 1986: 2.16426,
1987: 2.18071, 1988: 2.18856, 1989: 2.21250, 1990: 2.18377
}
D_var = 2.9584
agent_frames = []
np.random.seed(123)
for t in prod['market_ids'].unique():
market_agents = base_df.copy()
market_agents['market_ids'] = t
log_y = np.random.normal(loc=D_mean[int(t)], scale=np.sqrt(D_var), size=len(market_agents))
market_agents['income'] = np.exp(log_y)
agent_frames.append(market_agents)
final_agent_data = pd.concat(agent_frames, ignore_index=True)
To calculate the predicted market shares in the model, we cannot simply look at the "average" consumer. Instead, the model requires us to simulate a diverse population of consumers for every market and calculate each simulated person's probability of buying each car.
nodes (Unobserved Heterogeneity): These represent a consumer's hidden, personal tastes for vehicle characteristics. We need one node column for every random coefficient in your model. The original authors simply used a random number generator to pull 200 random points from a normal distribution, but we're following Conlon and Gortmaker and use Halton Draws to avoid simulation errors.
income : We draw random samples from a log-normal distribution that exactly matches the CPS mean and variance for that specific year.
pyblp.options.weights_tol = 1e-3
pyblp.options.verbose = False
if 'trend' not in prod.columns:
prod['trend'] = prod['market_ids'] - 1900
# We define Structural Formulations matching BLP (1995)
X1_formulation = pyblp.Formulation('1 + hpwt + air + mpd + space')
X2_formulation = pyblp.Formulation('1 + prices + hpwt + air + mpd + space')
X3_formulation = pyblp.Formulation('1 + log(hpwt) + air + log(mpg) + log(space) + trend')
product_formulations = (X1_formulation, X2_formulation, X3_formulation)
agent_formulation = pyblp.Formulation('0 + I(1 / income)')
# We initialize the PyBLP Problem
problem = pyblp.Problem(
product_formulations=product_formulations,
product_data=prod,
agent_formulation=agent_formulation,
agent_data=final_agent_data,
costs_type='log'
)
# We set Starting Values (from original paper Table IV)
initial_sigma = np.diag([3.612, 0.0, 4.628, 1.818, 1.050, 2.056])
pi_start = np.zeros((6, 1))
pi_start[1, 0] = -43.501 # Price interacts with 1/income
sigma_lb = np.zeros_like(initial_sigma)
sigma_ub = np.full_like(initial_sigma, np.inf)
sigma_lb[1, 1] = 0.0
sigma_ub[1, 1] = 0.0
# We solve using 2-Step GMM and L-BFGS-B Optimizer
optimization = pyblp.Optimization('l-bfgs-b', {'gtol': 1e-3})
results = problem.solve(
sigma=initial_sigma,
pi=pi_start,
sigma_bounds=(sigma_lb, sigma_ub),
optimization=optimization,
method='2s',
initial_update = True,
W_type = 'clustered',
se_type = 'clustered'
)
We now formating and print Table IV.
# Extract linear parameters
betas = results.beta.flatten()
beta_se = results.beta_se.flatten()
gammas = results.gamma.flatten()
gamma_se = results.gamma_se.flatten()
# Extract non-linear standard deviations
sigmas = np.diag(results.sigma)
sigma_se = np.diag(results.sigma_se)
# Extract price demographic interaction
alpha = results.pi[1, 0]
alpha_se = results.pi_se[1, 0]
demand_vars = ['Constant', 'HP/Weight', 'Air', 'MP$', 'Size']
supply_vars = ['Constant', 'ln(HP/Weight)', 'Air', 'ln(MPG)', 'ln(Size)', 'Trend']
sigma_mapped = [sigmas[0], sigmas[2], sigmas[3], sigmas[4], sigmas[5]]
sigma_se_mapped = [sigma_se[0], sigma_se[2], sigma_se[3], sigma_se[4], sigma_se[5]]
def fmt(est, se):
return f"{est:8.3f}\n({se:6.3f})"
print("\n" + "="*2)
print("TABLE IV: ESTIMATED PARAMETERS OF THE DEMAND AND PRICING EQUATIONS")
print("="*2)
# --- DEMAND EQUATION ---
demand_df = pd.DataFrame(index=demand_vars, columns=['Means (Beta)', 'Std. Dev. (Sigma)'])
for i, var in enumerate(demand_vars):
demand_df.at[var, 'Means (Beta)'] = fmt(betas[i], beta_se[i])
demand_df.at[var, 'Std. Dev. (Sigma)'] = fmt(sigma_mapped[i], sigma_se_mapped[i])
print("\n--- DEMAND EQUATION ---")
print(demand_df.to_string())
# --- TERM ON PRICE ---
print("\n--- TERM ON PRICE ---")
price_df = pd.DataFrame(index=['Price (Alpha)'], columns=['Estimate (Pi)'])
price_df.at['Price (Alpha)', 'Estimate (Pi)'] = fmt(alpha, alpha_se)
print(price_df.to_string(header=False))
# --- PRICING (SUPPLY) EQUATION ---
supply_df = pd.DataFrame(index=supply_vars, columns=['Estimate (Gamma)'])
for i, var in enumerate(supply_vars):
supply_df.at[var, 'Estimate (Gamma)'] = fmt(gammas[i], gamma_se[i])
print("\n--- PRICING (SUPPLY) EQUATION ---")
print(supply_df.to_string())
print("="*2)
==
TABLE IV: ESTIMATED PARAMETERS OF THE DEMAND AND PRICING EQUATIONS
==
--- DEMAND EQUATION ---
Means (Beta) Std. Dev. (Sigma)
Constant -5.577\n( 1.479) 3.750\n( 1.734)
HP/Weight 7.615\n( 1.857) 1.468\n( 6.998)
Air -4.073\n( 1.817) 7.270\n( 1.641)
MP$ 0.058\n( 0.515) 0.183\n( 1.542)
Size 2.945\n( 1.625) 4.595\n( 1.932)
--- TERM ON PRICE ---
Price (Alpha) -62.852\n(11.689)
--- PRICING (SUPPLY) EQUATION ---
Estimate (Gamma)
Constant 1.261\n( 0.138)
ln(HP/Weight) 0.744\n( 0.063)
Air 0.602\n( 0.035)
ln(MPG) -0.259\n( 0.063)
ln(Size) 0.128\n( 0.136)
Trend 0.017\n( 0.002)
==
This table is the structural "engine" of the paper. It displays the raw mathematical coefficients estimated by the GMM algorithm:
Demand Equation: Shows the average consumer's baseline preference for a feature ($\beta$), and how much individual consumers vary in that preference ($\sigma$).
Term on Price: Shows the critical price sensitivity coefficient ($\alpha$). Because it is interacted with consumer income ($1/y$), it ensures that wealthier consumers are less sensitive to high car prices than lower-income consumers.
Pricing (Supply) Equation: Shows the parameters ($\gamma$) for how much each physical characteristic (like horsepower or size) adds to the underlying manufacturing marginal cost of the vehicle.
Table V¶
Table V translates the parameter estimates from Table IV into practical, interpretable economic metrics by displaying estimated demand elasticities.
It focuses specifically on a representative sample of individual car models from the 1990 market.bFor each featured vehicle, the table lists its actual physical attributes and base price in the top row.
Directly beneath those values, the own-elasticities reveal exactly how sensitive that specific car's market share is to a hypothetical change in its own price or characteristics.
# We compute Own-Elasticities using PyBLP for Table V attributes
prod['elas_hpwt'] = results.extract_diagonals(results.compute_elasticities(name='hpwt'))
prod['elas_air'] = results.extract_diagonals(results.compute_elasticities(name='air'))
prod['elas_mpd'] = results.extract_diagonals(results.compute_elasticities(name='mpd'))
prod['elas_space'] = results.extract_diagonals(results.compute_elasticities(name='space'))
prod['elas_prices'] = results.extract_diagonals(results.compute_elasticities(name='prices'))
prod['elas_prices'] = prod['elas_prices'].abs()
# We filter strictly for the year 1990
df_1990 = prod[prod['market_ids'] == 1990].copy()
target_cars = ['323', 'SENT', 'ESCO', 'CAVA', 'ACCO', 'TAUR', 'MAXI', 'SEVI', 'LS40', '735']
selected_rows = []
for car in target_cars:
match = df_1990[df_1990['name'].str.contains(car, case=False, na=False)]
if not match.empty:
selected_rows.append(match.iloc[0])
sample_cars = pd.DataFrame(selected_rows)
# We format and Print the Replica of Table V
print("="*2)
print("TABLE V: A SAMPLE FROM 1990 OF ESTIMATED DEMAND ELASTICITIES")
print("="*2)
for _, row in sample_cars.iterrows():
print(f"{row['name']:<10} {row['hpwt']:>10.3f} {row['air']:>10.3f} {row['mpd']:>10.3f} {row['space']:>10.3f} {row['prices']:>10.3f}")
print(f"{'':<10} ({row['elas_hpwt']:>8.3f}) ({row['elas_air']:>8.3f}) ({row['elas_mpd']:>8.3f}) ({row['elas_space']:>8.3f}) ({row['elas_prices']:>8.3f})\n")
print("="*2)
==
TABLE V: A SAMPLE FROM 1990 OF ESTIMATED DEMAND ELASTICITIES
==
MZ323 0.366 0.000 3.645 1.075 5.049
( 2.805) ( -0.000) ( 0.376) ( 7.219) ( 6.932)
NISENT 0.391 0.000 3.645 1.092 5.661
( 2.996) ( -0.000) ( 0.362) ( 6.731) ( 6.252)
FDESCO 0.401 0.000 4.022 1.116 5.663
( 2.939) ( -0.000) ( 0.418) ( 6.849) ( 6.047)
CVCAVA 0.385 0.000 3.142 1.179 5.797
( 2.792) ( -0.000) ( 0.257) ( 7.904) ( 6.172)
HDACCO 0.457 0.000 3.016 1.255 9.292
( 3.269) ( -0.000) ( 0.165) ( 7.278) ( 4.256)
FDTAUR 0.304 0.000 2.262 1.334 9.671
( 2.129) ( -0.000) ( 0.069) ( 8.796) ( 4.257)
NIMAXI 0.518 1.000 2.513 1.300 13.695
( 4.003) ( 5.749) ( 0.257) ( 3.703) ( 6.109)
CDSEVI 0.517 1.000 2.011 1.374 24.353
( 3.847) ( 5.060) ( 0.197) ( 3.618) ( 4.460)
LXLS40 0.665 1.000 2.262 1.410 27.544
( 5.093) ( 4.907) ( 0.226) ( 4.097) ( 4.149)
BW735i 0.542 1.000 1.885 1.403 37.490
( 4.200) ( 5.321) ( 0.183) ( 4.226) ( 3.744)
==
It focuses on a representative sample of specific cars from the year 1990, ranging from cheap compacts (Mazda 323) to luxury sedans (BMW 735i).
For each car, it lists the actual physical attributes and the price on the top row.
On the bottom row, it provides the Own-Elasticities. This tells you exactly how sensitive a specific car's market share is to a 1% change in its own price or characteristics. For example, it proves mathematically that cheap, compact cars have highly elastic (price-sensitive) demand, while luxury cars have inelastic (price-insensitive) demand.
Table VIII¶
The primary purpose of Table VIII is to prove that the model respects a fundamental law of microeconomics: profit margins are driven by price sensitivity. Since, they assume a Bertrand-Nash oligopoly, the markups should mirror the elasticities of the previous table.
# We compute the implied markups for every car in the dataset
# This relies on the Bertrand-Nash supply-side inversion we estimated
markups = results.compute_markups().flatten()
# We calculate Marginal Costs (MC = Price - Markup)
costs = prod['prices'] - markups
df_econ = prod[['market_ids', 'name', 'prices']].copy()
df_econ['markups'] = markups
df_econ['costs'] = costs
# We isolate the year 1990
df_1990 = df_econ[df_econ['market_ids'] == 1990].reset_index(drop=True)
# We extract the data for these specific cars
table_vi_data = []
for car in target_cars:
match = df_1990[df_1990['name'].str.contains(car, case=False, na=False)]
row = match.iloc[0]
table_vi_data.append({
'Car Model': car,
'Price ($)': row['prices'],
'Marginal Cost ($)': row['costs'],
'Markup ($)': row['markups']
})
table_vi = pd.DataFrame(table_vi_data).set_index('Car Model')
table_vi_formatted = (table_vi * 1000).astype(int)
print("="*65)
print("TABLE VIII: ESTIMATED PRICE-MARGINAL COST MARKUPS (1990)")
print("="*65)
print(table_vi_formatted.to_string())
print("="*65)
=================================================================
TABLE VIII: ESTIMATED PRICE-MARGINAL COST MARKUPS (1990)
=================================================================
Price ($) Marginal Cost ($) Markup ($)
Car Model
323 5048 4903 145
SENT 5661 5498 163
ESCO 5663 5468 195
CAVA 5797 5564 232
ACCO 9292 9050 241
TAUR 9671 9410 260
MAXI 13694 13527 167
SEVI 24353 24055 297
LS40 27543 27299 244
735 37490 37215 275
=================================================================
