knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )
BLPestimatoR
provides an efficient estimation algorithm to perform the demand estimation described in @BLP1995. The routine uses analytic gradients and offers a large number of optimization routines and implemented integration methods as discussed in @Brunner2017.
This extended documentation demonstrates the steps of a typical demand estimation with the package:
BLP_data
(includes the specification of a model and providing integration draws for observed or unobserved heterogeneity)estimate_BLP
summary
get_elasticities
For this purpose the well-known training datasets for the cereal market [@Nevo2001] and the car market [@BLP1995] are included in the package. Loading the package is therefore the very first step of the demand estimation:
library(BLPestimatoR)
Since version 0.1.6 the model is provided in R’s formula syntax and consists of five parts. The variable
to be explained is given by observed market shares. Explanatory variables are grouped into four (possibly
overlapping) categories separated by |
:
The first part of this documentation starts with the cereal data example from @Nevo2001. Nevo's model can be translated into the following formula syntax:
nevos_model <- as.formula("share ~ price + productdummy | 0+ productdummy | price + sugar + mushy | 0+ IV1 + IV2 + IV3 + IV4 + IV5 + IV6 + IV7 + IV8 + IV9 + IV10 + IV11 + IV12 + IV13 + IV14 + IV15 + IV16 + IV17 + IV18 + IV19 + IV20")
The model is directly related to consumer $i$'s indirect utility from purchasing cereal $j$ in market $t$:
$$u_{ijt}=\sum_{m=1}^M x^{(m)}{jt} \beta{i,m}+\xi_{jt}+\epsilon_{ijt} \;\; \text{with}$$ $$\beta_{i,m}= \bar{\beta}m + \sum{r=1}^R \gamma_{m,r} d_{i,r} + \sigma_m \nu_{i,m}$$ and
price
, sugar
, mushy
and an intercept)income
, incomesq
, age
, child
)$$\theta_2 = \begin{pmatrix} \sigma_1 & \gamma_{1,1} & \cdots & \gamma_{1,R} \ \sigma_2 & \gamma_{2,1} & \cdots & \gamma_{2,R} \ \vdots & \vdots & \ddots & \vdots \ \sigma_M & \gamma_{M,1} & \cdots & \gamma_{M,R} \end{pmatrix}$$
Product related variables are collected in the dataframe productData
with the following requirements:
market_identifier
)A variable that uniquely identifies a product in a market (product_identifier
) is optional, but enhances clarity (interpreting elasticities, for example, is much easier). market_identifier
and product_identifier
together uniquely identify an observation, which is used by the function update_BLP_data
to update any variable in the data (in this case product_identifier
is mandatory).
In the cereal example, this gives the following dataframe:
head(productData_cereal)
The arguments related to the numerical integration problem are of particular importance when providing own integration draws and weights, which is most relevant for observed heterogeneity (for unobserved heterogeneity, the straightforward approach is the use of automatic integration).
In the cereal data, both, observed and unobserved heterogeneity, is used for the random coefficients. Starting with observed heterogeneity, user provided draws are collected in a list. Each list entry must be named according to the name of a demographic. Each entry contains the following variables:
market_identifier
that matches each line to a market (same variable name as in productData
)In the cereal example, observed heterogeneity is provided as follows (list names correspond to the demographics):
demographicData_cereal$income[1:4, 1:5] demographicData_cereal$incomesq[1:4, 1:5] demographicData_cereal$age[1:4, 1:5] demographicData_cereal$child[1:4, 1:5]
If demographic input (demographicData
) is missing, the estimation routine considers only coefficients for unobserved heterogeneity. This can be done by already implemented integration methods via integration_method
as shown in the estimation section. In Nevo's cereal example however, a specific set of 20 draws is
given. For this situation, draws are also provided as a list (list names correspond to the formula's random coefficients and each list entry has a variable market_identifier
):
originalDraws_cereal$constant[1:4, 1:5] # renaming constants: names(originalDraws_cereal)[1] <- "(Intercept)" originalDraws_cereal$price[1:4, 1:5] originalDraws_cereal$sugar[1:4, 1:5] originalDraws_cereal$mushy[1:4, 1:5]
As demonstrated above, list entries for draws of constants must be named (Intercept)
. Other names of list entries must match the random coefficients specified in the formula.
Calling BLP_data
structures and prepares the data for estimation and creates the data object:
productData_cereal$startingGuessesDelta <- c(log(w_guesses_cereal)) # include orig. draws in the product data cereal_data <- BLP_data( model = nevos_model, market_identifier = "cdid", par_delta = "startingGuessesDelta", product_identifier = "product_id", productData = productData_cereal, demographic_draws = demographicData_cereal, blp_inner_tol = 1e-6, blp_inner_maxit = 5000, integration_draws = originalDraws_cereal, integration_weights = rep(1 / 20, 20) )
The arguments in greater detail:
model
provides the utility model as explained above
market_identifier
gives the name of the variable in productData
that matches
each observation to a market
product_identifier
gives the name of the variable in productData
that matches
each observation to a product (must be unique in a market)
productData
is given as a dataframe and demographicData
as a list as described above
par_delta
gives the name of the variable in productData
for mean utilities
blp_inner_tol
, blp_inner_maxit
: arguments related to be BLP algorithm include the convergence threshold and the maximum number of iterations in the contraction mapping
if integration draws are provided manually, integration_draws
and
integration_weights
need to be specified
for automatic integration the user specifies integration_method
, for example
integration_method= "MLHS"
, and the accuracy of the integration method by
integration_accuracy
(for stochastic integration methods this equals the number
of draws)
If you decide to update your data later, you can use the function update_BLP_data
.
The provided set of starting guesses par_theta2
is matched with formula input and demographic data:
par_theta2
must match with the random coefficients specified in the formula (note: constants must be named (Intercept)
) par_theta2
must match with list entry names of demographicData
and a column for unobserved heterogeneity (must be named `unobs_sd)NA
s in par_theta2
indicate the exclusion from estimation, i.e. the coefficient is assumed to be zero.These requirements are demonstrated with a set of exemplary starting guesses:
# before: theta_guesses_cereal theta_guesses_cereal[theta_guesses_cereal == 0] <- NA colnames(theta_guesses_cereal) <- c("unobs_sd", "income", "incomesq", "age", "child") rownames(theta_guesses_cereal) <- c("(Intercept)", "price", "sugar", "mushy") # correctly named: theta_guesses_cereal
The following code performs the demand estimation:
cereal_est <- estimateBLP( blp_data = cereal_data, par_theta2 = theta_guesses_cereal, solver_method = "BFGS", solver_maxit = 1000, solver_reltol = 1e-6, standardError = "heteroskedastic", extremumCheck = FALSE, printLevel = 1 ) summary(cereal_est)
The arguments in greater detail:
par_theta2
gives initial values for non-linear parameters to be optimized
over. Correct naming of columns and rows is important to allow correct matching.
solver_method
, solver_maxit
, solver_reltol
: solver related arguments that specify the R internal optimization (optim
function). Additional arguments can be passed to optim via ...
standardError
can be specified as homoskedastic
, heteroskedastic
or
cluster
. The latter requires the variable group_structure
in productData
giving the related
cluster.
if extremumCheck
is TRUE
, numerical derivatives at the solver optimum are used
to check, if a local minimum was found
printLevel
controls for the amount of information that is provided during the
estimation
Many of these arguments have default values. In the following setting you see a minimum of necessary arguments with an automatic generation of integration draws and just unobserved heterogeneity. The summary output informs you about the most important default values.
cereal_data2 <- BLP_data( model = nevos_model, market_identifier = "cdid", product_identifier = "product_id", productData = productData_cereal, integration_method = "MLHS", integration_accuracy = 20, integration_seed = 213 ) cereal_est2 <- estimateBLP(blp_data = cereal_data2, printLevel = 1) summary(cereal_est2)
Standard errors can be computed with three options that control for the unobserved characteristic $\xi$, which consists of $N$ elements. $\Omega$ denotes the variance covariance matrix of $\xi$.
option homoskedastic
requires the standard deviation $\sigma_i$ for each $\xi_i \;\forall
i\in 1,\cdots,N$ to be identical:
$$\Omega = \begin{pmatrix} \sigma & 0 &
\dots & 0\ 0 & \sigma & & 0\ \vdots & & \ddots & 0\
0 & 0 & 0 & \sigma \ \end{pmatrix}$$
option heteroskedastic
allows for individual standard deviations $\sigma_i$ for each
$\xi_i$ :
$$\Omega = \begin{pmatrix} \sigma_1 & 0 &
\dots & 0\
0 & \sigma_2 & & 0\
\vdots & & \ddots & 0\
0 & 0 & 0 & \sigma_N \ \end{pmatrix}$$
option cluster
allows for cluster individual variance covariance matrices in each of $M$ cluster groups. For this option the argument group_structure
needs to be specified in the function BLP_data
to determine the cluster group. This gives the block-diagonal form with $\Sigma_m$ as the variance covariance matrix for all $\xi_i$ in cluster $m$:
$$\Omega = \begin{pmatrix} \Sigma_1 & 0 & \dots & 0\ 0 & \Sigma_2 & & 0\ \vdots & & \ddots & 0\ 0 & 0 & 0 & \Sigma_M \ \end{pmatrix}$$
The following code demonstrates the calculation of elasticities for the estimation object cereal_est
.
# extract parameters from output theta1_price <- cereal_est$theta_lin["price", ] theta2 <- matrix(NA, nrow = 4, ncol = 5) colnames(theta2) <- c("unobs_sd", "income", "incomesq", "age", "child") rownames(theta2) <- c("(Intercept)", "price", "sugar", "mushy") for (i in 1:13) { theta2[cereal_est$indices[i, 1], cereal_est$indices[i, 2]] <- cereal_est$theta_rc[i] } delta_data <- data.frame( "product_id" = cereal_data$parameters$product_id, "cdid" = cereal_data$parameters$market_id_char_in, "startingGuessesDelta" = cereal_est$delta ) # always use update_BLP_data() to update data object to maintain consistent data cereal_data <- update_BLP_data( data_update = delta_data, blp_data = cereal_data ) shareObj <- getShareInfo( blp_data = cereal_data, par_theta2 = theta2, printLevel = 1 ) get_elasticities( blp_data = cereal_data, share_info = shareObj, theta_lin = theta1_price, variable = "price", products = c("cereal_1", "cereal_4"), market = "market_2" )
The value of the elasticity matrix in row $j$ and column $i$ for a variable $x$, gives the effect of a change in product $i$'s characteristic $x$ on the share of product $j$.
Further analysis like incorporating a supply side or performing a merger simulation often requires access to building blocks of the BLP algorithm. The following wrappers insure correct data inputs and access the internal functions of the algorithm.
In the following, you find an example of the contraction mapping and an evaluation of the GMM function at the starting guess:
delta_eval <- getDelta_wrap( blp_data = cereal_data, par_theta2 = theta_guesses_cereal, printLevel = 4 ) productData_cereal$startingGuessesDelta[1:6] delta_eval$delta[1:6] delta_eval$counter gmm <- gmm_obj_wrap( blp_data = cereal_data, par_theta2 = theta_guesses_cereal, printLevel = 2 ) gmm$local_min
Printed distances in the contraction mapping are maximum absolute distances between the current vector of mean utilities and the previous one.
For any $\theta_2$, you can compute predicted shares:
shareObj <- getShareInfo( blp_data = cereal_data, par_theta2 = theta_guesses_cereal, printLevel = 4 ) shareObj$shares[1:6]
The object contains a list of outputs that are useful for further economic analysis. For example, the list element sij
contains share probabilities for every individual and needs to be given to calculate elasticities.
The gradient contains two important building blocks as explained in the appendix of @Nevo2001:
$\frac{\partial s_{ijt}}{\partial \theta_2}$ , i.e. the derivative of individual $i$'s share of product $j$ in market $t$ with respect to non-linear parameters
$\frac{\partial s_{ijt}}{\partial \delta}$ , i.e. the derivative of individual $i$'s share of product $j$ in market $t$ with respect to mean utilities
Both are used to compute the jacobian and are easy to obtain with the package as the following example demonstrates:
# market 2: derivatives1 <- dstdtheta_wrap( blp_data = cereal_data, par_theta2 = theta_guesses_cereal, market = "market_2" ) derivatives2 <- dstddelta_wrap( blp_data = cereal_data, par_theta2 = theta_guesses_cereal, market = "market_2" ) jac_mkt2 <- -solve(derivatives2) %*% derivatives1 jac_mkt2[1:5, 1:4] # all markets jacobian_nevo <- getJacobian_wrap( blp_data = cereal_data, par_theta2 = theta_guesses_cereal, printLevel = 2 ) jacobian_nevo[25:29, 1:4] # compare to jac_mkt2
Analyzing a hypothetical merger is demonstrated by the car data of @BLP1995. In this case, the preparation of product data comprises the computation of instruments as a function of product characteristics of competitors' products (for details, check @BLP1995). This example is based on data and documentation of @KM2014.
# add owner matix to productData own_pre <- dummies_cars colnames(own_pre) <- paste0("company", 1:26) productData_cars <- cbind(productData_cars, own_pre) # construct instruments nobs <- nrow(productData_cars) X <- data.frame( productData_cars$const, productData_cars$hpwt, productData_cars$air, productData_cars$mpg, productData_cars$space ) sum_other <- matrix(NA, nobs, ncol(X)) sum_rival <- matrix(NA, nobs, ncol(X)) sum_total <- matrix(NA, nobs, ncol(X)) for (i in 1:nobs) { other_ind <- productData_cars$firmid == productData_cars$firmid[i] & productData_cars$cdid == productData_cars$cdid[i] & productData_cars$id != productData_cars$id[i] rival_ind <- productData_cars$firmid != productData_cars$firmid[i] & productData_cars$cdid == productData_cars$cdid[i] total_ind <- productData_cars$cdid == productData_cars$cdid[i] sum_other[i, ] <- colSums(X[other_ind == 1, ]) sum_rival[i, ] <- colSums(X[rival_ind == 1, ]) sum_total[i, ] <- colSums(X[total_ind == 1, ]) } colnames(sum_other) <- paste0("IV", 1:5) colnames(sum_rival) <- paste0("IV", 6:10) productData_cars <- cbind(productData_cars, sum_other, sum_rival) head(productData_cars) # To show similarities between implementations of other authors, # the variable "const" is used, although constants are considered by default. blps_model <- as.formula("share ~ 0 + const + price + hpwt + air + mpg + space | 0 + const + hpwt + air + mpg + space | 0 + price + const + hpwt + air + mpg | 0 + IV1 + IV2 + IV3 + IV4 + IV5 + IV6 + IV7 + IV8 + IV9 + IV10") car_data <- BLP_data( model = blps_model, market_identifier = "cdid", product_identifier = "id", additional_variables = paste0("company", 1:26), # check reordering works productData = productData_cars, blp_inner_tol = 1e-9, blp_inner_maxit = 5000, integration_method = "MLHS", integration_accuracy = 50, integration_seed = 48 )
In the next step, starting guesses for random coefficients are generated from a standard normal distribution. The estimation of the model works like before.
set.seed(121) theta_guesses <- matrix(rnorm(5)) rownames(theta_guesses) <- c("price", "const", "hpwt", "air", "mpg") colnames(theta_guesses) <- "unobs_sd" car_est <- estimateBLP( blp_data = car_data, par_theta2 = theta_guesses, solver_method = "BFGS", solver_maxit = 1000, solver_reltol = 1e-6, extremumCheck = FALSE, printLevel = 0 ) summary(car_est)
Next, all parameters that are required by the subsequent merger analysis are extracted. Note that all extracted data is based on the estimation object car_est
or the data object car_data
to maintain data consistency (for example, the order of data in product_data_cars
might differ from car_data
). Moreover, mean utilities are updated in car_data
by the values in the estimation object car_est
.
## Pre-Merger data own_pre <- as.matrix(car_data$data$additional_data[, paste0("company", 1:26)]) delta_pre <- car_est$delta theta1_price <- car_est$theta_lin["price", ] theta2_price <- car_est$theta_rc["unobs_sd*price"] theta2_all <- matrix(car_est$theta_rc) rownames(theta2_all) <- c("price", "const", "hpwt", "air", "mpg") colnames(theta2_all) <- "unobs_sd" ## update mean utility in data ( always use update_BLP_data() to update data object to maintain consistent data ) delta_data <- data.frame( "id" = car_data$parameters$product_id, "cdid" = car_data$parameters$market_id, "delta" = delta_pre ) car_data_updated <- update_BLP_data( data_update = delta_data, blp_data = car_data )
In the next step, an estimate for marginal costs $mc$ before the merger is computed. The following is based on the FOC of a Bertrand equilibrium with prices $p$ before the merger:
$$ p^{pre} - \widehat{mc} = \Omega^{pre}(p^{pre})^{-1} \hat{s}(p^{pre}) $$ $\Omega^{pre}(p^{pre})^{-1}$ is defined marketwise as the inverse of
$$ \Omega^{pre}(p^{pre}) = \pmatrix{ -\frac{\partial s_{1}}{\partial p_{1}} (p^{pre}) \cdot D_{1,1} & -\frac{\partial s_{2}}{\partial p_{1}} (p^{pre}) \cdot D_{1,2} & \cdots & -\frac{\partial s_{j}}{\partial p_{1}} (p^{pre}) \cdot D_{1,j} & \cdots & -\frac{\partial s_{J}}{\partial p_{1}} (p^{pre}) \cdot D_{1,J}\ -\frac{\partial s_{1}}{\partial p_{2}} (p^{pre}) \cdot D_{2,1} & -\frac{\partial s_{2}}{\partial p_{2}} (p^{pre}) \cdot D_{2,2} & \cdots & -\frac{\partial s_{j}}{\partial p_{2}} (p^{pre}) \cdot D_{2,j} & \cdots & -\frac{\partial s_{J}}{\partial p_{2}} (p^{pre}) \cdot D_{2,J}\ \vdots & \vdots & \ddots & \vdots & \ddots & \vdots\ -\frac{\partial s_{1}}{\partial p_{k}} (p^{pre}) \cdot D_{k,1} & -\frac{\partial s_{2}}{\partial p_{k}} (p^{pre}) \cdot D_{k,2} & \cdots & -\frac{\partial s_{j}}{\partial p_{k}} (p^{pre}) \cdot D_{k,j} & \cdots & -\frac{\partial s_{J}}{\partial p_{k}} (p^{pre}) \cdot D_{k,J}\ \vdots & \vdots & \ddots & \vdots & \ddots & \vdots\ -\frac{\partial s_{1}}{\partial p_{J}} (p^{pre}) \cdot D_{J,1} & -\frac{\partial s_{2}}{\partial p_{J}} (p^{pre}) \cdot D_{J,2} & \cdots & -\frac{\partial s_{j}}{\partial p_{J}} (p^{pre}) \cdot D_{J,j}& \cdots & -\frac{\partial s_{J}}{\partial p_{J}} (p^{pre}) \cdot D_{J,J}\ } $$
with $$D_{k,j} = \begin{cases} 1 & \text{if products k and j are produced by the same firm} \ 0 & \text{otherwise} \ \end{cases}$$
Partial derivatives $\frac{\partial s_j}{\partial p_k}$ can be calculated based on the elasticity $\eta_{jk} = \frac{\partial s_j }{\partial p_k }\frac{ p_k}{ s_j}$, so $$ \frac{\partial s_j}{\partial p_k} = \eta_{jk} \cdot \frac{ s_j}{ p_k} $$
In the following code chunk, these objects in a market i
are labeled as follows:
own_prod_pre_i
($D_{k,j}$) elasticities_i
($\eta_{jk}$)derivatives_i
($\eta_{jk} \cdot \frac{ s_j}{ p_k}$)-solve(t(derivatives_i) * own_prod_pre_i)
($\Omega^{pre}(p^{pre})^{-1}$)shareObj$shares
($\hat{s}(p^{pre})$).## calculate sij shareObj <- getShareInfo( blp_data = car_data_updated, par_theta2 = theta2_all, printLevel = 0 ) ## computation of marginal costs market_id <- car_data$parameters$market_id nmkt <- length(unique(market_id)) markups <- numeric(length(market_id)) sh <- shareObj$shares prices_pre <- car_data$data$X_rand[, "price"] for (i in 1:nmkt) { mkt_ind <- market_id == i share_i <- sh[ mkt_ind ] price_pre_i <- prices_pre[ mkt_ind ] scalar_i <- matrix(1 / share_i) %*% matrix(price_pre_i, nrow = 1) elasticities_i <- get_elasticities( blp_data = car_data_updated, share_info = shareObj, theta_lin = theta1_price, variable = "price", market = i, printLevel = 0 ) derivatives_i <- elasticities_i / scalar_i # partial derivatives of shares wrt price own_pre_i <- own_pre[ mkt_ind, ] own_prod_pre_i <- own_pre_i %*% t(own_pre_i) # if element (i,j) equals 1, that means that prod i and j are produced by same firm markups[mkt_ind] <- c(-solve(t(derivatives_i) * own_prod_pre_i) %*% share_i) } marg_cost <- prices_pre - markups
The ownership matrix is adjusted to implement a hypothetical merger between Chrysler and GM:
# Merger between company 16 and 19 (i.e. GM and Chrysler) prices_post <- numeric(2217) own_post <- cbind( own_pre[, 1:15], own_pre[, 16] + own_pre[, 19], own_pre[, 17:18], own_pre[, 20:26] )
To analyze the effect on prices the FOC of the new equilibrium must be solved: $$ p^{post} - \widehat{mc} = \Omega^{post}(p^{post})^{-1} \hat{s}(p^{post}) $$
The solution of this set of non-linear equations is obtained by the function foc_bertrand_mkt
and the package nleqslv
:
foc_bertrand_mkt <- function(par, own_prod, blp_data, mkt, marg_cost, theta_lin, theta_rc) { # argument par: candidate for post merger prices # arguments own_prod, blp_data, mkt, marg_cost, theta_lin, theta_rc: see previous code blocks # post merger updates: update the BLP_data object for market i tmp <- data.frame( "id" = blp_data$parameters$product_id, "cdid" = blp_data$parameters$market_id, "delta" = blp_data$data$delta, "price" = blp_data$data$X_rand[, "price"] ) market_ind <- blp_data$parameters$market_id == mkt delta_old <- blp_data$data$delta prices_pre <- blp_data$data$X_rand[, "price"] tmp$price[ market_ind ] <- par tmp$delta[ market_ind ] <- delta_old[market_ind] - prices_pre[market_ind] * theta_lin + par * theta_lin new_blp_data <- update_BLP_data( blp_data = blp_data, data_update = tmp ) ShareObj <- getShareInfo( blp_data = new_blp_data, par_theta2 = theta_rc, printLevel = 0 ) implied_shares <- as.matrix(ShareObj$shares[market_ind]) elasticities_post_mkt <- get_elasticities( blp_data = new_blp_data, share_info = ShareObj, theta_lin = theta_lin, variable = "price", market = mkt, printLevel = 0 ) scalar_mkt <- matrix(1 / implied_shares) %*% matrix(par, nrow = 1) derivatives_mkt <- elasticities_post_mkt / scalar_mkt markups_post <- c(-solve(t(derivatives_mkt) * own_prod) %*% implied_shares) differences <- par - marg_cost[market_ind] - markups_post return(differences) }
Finally, the function is used to compute the new equilibrium:
library(nleqslv) # to solve non linear first order conditions for (i in 1:nmkt) { mkt_ind <- market_id == i own_post_i <- own_post[ mkt_ind, ] own_prod_post_i <- own_post_i %*% t(own_post_i) price_pre_i <- prices_pre[ mkt_ind ] solution <- nleqslv( x = price_pre_i, foc_bertrand_mkt, # startingguesses: price_pre_i own_prod = own_prod_post_i, blp_data = car_data_updated, mkt = i, marg_cost = marg_cost, theta_lin = theta1_price, theta_rc = theta2_all ) prices_post[ market_id == i ] <- solution$x }
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.