This package is in beta. The API may change in future versions. Please use the GitHub issue tracker to report bugs or to request features.
The pyblp package is a Python 3 implementation of routines for estimating the demand for differentiated products with BLPtype random coefficients logit models. This package was created by Jeff Gortmaker in collaboration with Chris Conlon.
Development of the package has been guided by the work of many researchers and practitioners. For a full list of references, including the original work of Berry, Levinsohn, and Pakes (1995), refer to the references section of the documentation.
The pyblp package has been tested on Python versions 3.6 and 3.7. The SciPy instructions for installing related packages is a good guide for how to install a scientific Python environment. A good choice is the Anaconda Distribution, since, along with many other packages that are useful for scientific computing, it comes packaged with pyblpβs only required dependencies: NumPy, SciPy, SymPy, and Patsy.
However, pyblp may not work with old versions of its dependencies. You can update pyblpβs dependencies in Anaconda with:
condaupdatenumpyscipysympypatsy
You can install the current release of pyblp with pip:
pipinstallpyblp
You can upgrade to a newer release with the upgrade flag:
pipinstallupgradepyblp
If you lack permissions, you can install pyblp in your user directory with the user flag:
pipinstalluserpyblp
Alternatively, you can download a wheel or source archive from PyPI. You can find the latest development code on GitHub and the latest development documentation here.
The following sections provide a very brief overview of the BLP model and how it is estimated. This goal is to concisely introduce the notation and terminology used throughout the rest of the documentation. For a more indepth overview, refer to Conlon and Gortmaker (2019).
There are \(t = 1, 2, \dotsc, T\) markets, each with \(j = 1, 2, \dotsc, J_t\) products produced by \(f = 1, 2, \dotsc, F_t\) firms, for a total of \(N\) products across all markets. There are \(i = 1, 2, \dotsc, I_t\) agents who choose among the \(J_t\) products and an outside good \(j = 0\).
Observed demandside product characteristics are contained in the \(N \times K_1\) matrix of linear characteristics, \(X_1\), and the \(N \times K_2\) matrix of nonlinear characteristics, \(X_2\), which is typically a subset of \(X_1\). Unobserved demandside product characteristics, \(\xi\), are a \(N \times 1\) vector.
In market \(t\), observed agent characteristics are a \(I_t \times D\) matrix called demographics, \(d\). Unobserved agent characteristics are a \(I_t \times K_2\) matrix, \(\nu\).
The indirect utility of agent \(i\) from purchasing product \(j\) in market \(t\) is
The \(K_1 \times 1\) vector of demandside linear paramterers, \(\beta\), is partitioned into two components: \(\alpha\) is a \(K_1^p \times 1\) vector of parameters on the \(N \times K_1^p\) submatrix of endogenous characteristics, \(X_1^p\), and \(\beta^x\) is a \(K_1^x \times 1\) vector of parameters on the \(N \times K_1^x\) submatrix of exogenous characteristics, \(X_1^x\). Usually, \(X_1^p = p\), prices, so \(\alpha\) is simply a scalar.
The agentspecific portion of utility in a single market is
The model incorporates both observable (demographic) and unobservable taste heterogeneity though random coefficients. For the unobserved heterogeneity, we let \(\nu\) denote independent draws from the standard normal distribution. These are scaled by a \(K_2 \times K_2\) upper triangular matrix \(\Sigma\), which denotes the Cholesky root of the covariance matrix for unobserved taste heterogeneity. The \(K_2 \times D\) matrix \(\Pi\) measures how agent tastes vary with demographics.
Random idiosyncratic preferences, \(\epsilon_{jti}\), are assumed to be Type I Extreme Value, so that conditional on the heterogeneous coefficients, marketshares follow the well known logit form. Aggregate marketshares are obtained by integrating over the distribution of individual heterogeneity. They are approximated with Monte Carlo integration or quadrature rules defined by the \(I_t \times K_2\) matrix of integration nodes, \(\nu\), and a \(I_t \times 1\) vector of integration weights, \(w\):
Observed supplyside product characteristics are contained in the \(N \times K_3\) matrix of cost characteristics, \(X_3\). Prices cannot be cost characteristics, but nonprice product characteristics often overlap with the demandside characteristics in \(X_1\) and \(X_2\). Unobserved supplyside product characteristics, \(\omega\), are a \(N \times 1\) vector.
Firm \(f\) chooses prices in market \(t\) to maximize the profits of its products \(\mathscr{J}_{ft} \subset \{1, 2, \ldots, J_t\}\):
Here, \(O\) denotes the marketlevel ownership matrix, where \(O_{jk}\) is typically \(1\) if the same firm produces products \(j\) and \(k\), and \(0\) otherwise.
To include a supply side, we must specify a functional form for marginal costs:
A demand side is always estimated but including a supply side is optional. With only a demand side, there are three sets of parameters to be estimated: \(\beta\) (which may include \(\alpha\)), \(\Sigma\) and \(\Pi\). With a supply side, there is also \(\gamma\). The linear parameters, \(\beta\) and \(\gamma\), are typically concentrated out of the problem. The exception is \(\alpha\), which cannot be concentrated out when there is a supply side because it is needed to compute demand derivatives and hence marginal costs. Linear parameters that are not concentrated out along with unknown nonlinear parameters in \(\Sigma\) and \(\Pi\) are collectively denoted \(\theta\), a \(P \times 1\) vector.
in which \(W\) is a \((M_D + M_S) \times (M_D + M_S)\) weighting matrix and \(\bar{g}\) is a \((M_D + M_S) \times 1\) vector of sample moment conditions:
in which \(Z_D\) and \(Z_S\) are \(N \times M_D\) and \(N \times M_S\) matrices of demand and supplyside instruments containing excluded instruments along with \(X_1^x\) and \(X_3\), respectively.
The vector \(\bar{g}\) is the sample analogue of the GMM moment conditions \(E[g_{jt}] = 0\) where
Given a \(\hat{\theta}\), the first step to computing the objective \(q(\hat{\theta})\) is to compute \(\delta(\hat{\theta})\) in each market with the following standard contraction:
where \(s\) are the marketβs observed shares and \(s(\delta, \hat{\theta})\) are calculated marketshares. Iteration terminates when the norm of the change in \(\delta(\hat{\theta})\) is less than a small number.
With a supply side, marginal costs are then computed according to (7):
With only a demand side, \(\alpha\) can be concentrated out, so \(X = X_1\), \(Z = Z_D\), and \(Y = \delta(\hat{\theta})\) recover the full \(\hat{\beta}\) in (15).
Finally, the unobserved product characteristics (structural errors),
are interacted with the instruments to form the sample moment conditions \(\bar{g}(\hat{\theta})\) in (11), which give the GMM objective \(q(\hat{\theta})\) in (10).
where \(k_1, k_2, \dotsc, k_{E_D}\) and \(\ell_1, \ell_2, \dotsc, \ell_{E_S}\) index unobserved characteristics that are fixed across \(E_D\) and \(E_S\) dimensions. For example, with \(E_D = 1\) dimension of product fixed effects, \(\xi_{jt} = \xi_j + \Delta\xi_{jt}\).
Small numbers of fixed effects can be estimated with dummy variables in \(X_1\), \(X_3\), \(Z_D\), and \(Z_S\). However, this approach does not scale with high dimensional fixed effects because it requires constructing and inverting an infeasibly large matrix in (15).
Instead, fixed effects are typically absorbed into \(X\), \(Z\), and \(Y(\hat{\theta})\) in (15). With one fixed effect, these matrices are simply demeaned within each level of the fixed effect. Both \(X\) and \(Z\) can be demeaned just once, but \(Y(\hat{\theta})\) must be demeaned for each new \(\hat{\theta}\).
This procedure is equivalent to replacing each column of the matrices with residuals from a regression of the column on the fixed effect. The FrishWaughLovell (FWL) theorem of Frisch and Waugh (1933) and Lovell (1963) guarantees that using these residualized matrices gives the same results as including fixed effects as dummy variables. When \(E_D > 1\) or \(E_S > 1\), the matrices are residualized with more involved algorithms.
Once fixed effects have been absorbed, estimation is as described above with the structural errors \(\Delta\xi\) and \(\Delta\omega\).
Incorporating parameters that measure within nesting group correlation gives the random coefficients nested logit (RCNL) model of Brenkers and Verboven (2006) and Grigolon and Verboven (2014). There are \(h = 1, 2, \dotsc, H\) nesting groups and each product \(j\) is assigned to a group \(h(j)\). The set \(\mathscr{J}_{ht} \subset \{1, 2, \ldots, J_t\}\) denotes the products in group \(h\) and market \(t\).
In the RCNL model, idiosyncratic preferences are partitioned into
where \(\bar{\epsilon}_{jti}\) is Type I Extreme Value and \(\bar{\epsilon}_{h(j)ti}\) is distributed such that \(\epsilon_{jti}\) is still Type I Extreme Value.
The nesting parameters, \(\rho\), can either be a \(H \times 1\) vector or a scalar so that for all groups \(\rho_h = \rho\). Letting \(\rho \to 0\) gives the standard BLP model and \(\rho \to 1\) gives division by zero errors. With \(\rho_h \in (0, 1)\), the expression for choice probabilities in (5) becomes more complicated:
In both models, a supply side can still be estimated jointly with demand. Estimation is as described above with a representative agent in each market: \(I_t = 1\) and \(w_1 = 1\).
Counterfactual evaluation, synthetic data simulation, and optimal instrument generation often involve solving for prices implied by the Bertrand first order conditions in (7). Solving this system with Newtonβs method is slow and iterating over \(p \leftarrow c + \eta(p)\) may not converge because it is not a contraction.
which, unlike (7), is a contraction. Iteration terminates when the norm of firmsβ first order conditions, \(\Lambda(p)(p  c  \zeta(p))\), is less than a small number.
This section uses a series of Jupyter Notebooks to explain how pyblp can be used to solve example problems, compute postestimation outputs, and simulate problems. Other examples can be found in the API Documentation.
In this tutorial, weβll use data from Nevo (2000) to solve the paperβs fake cereal problem. Locations of CSV files that contain the data are in the data module.
We will compare two simple models, the plain (IIA) logit model and the nested logit (GEV) model using the fake cereal dataset of Nevo (2000).
where \(\epsilon_{jti}\) is distributed IID with the Type I Extreme Value (Gumbel) distribution. It is common to normalize the mean utility of the outside good to zero so that \(U_{0ti} = \epsilon_{0ti}\). This gives us aggregate marketshares
A Logit Problem can be created by simply excluding the formulation for the nonlinear parameters, \(X_2\), along with any agent information. In other words, it requires only specifying the linear component of demand.
Weβll set up and solve a simple version of the fake data cereal problem from Nevo (2000). Since we wonβt include any nonlinear characteristics or parameters, we donβt have to worry about configuring an optimization routine.
There are some important reserved variable names:
market_ids are the unique market identifiers which we subscript with \(t\).
shares specifies the marketshares which need to be between zero and one, and within a market ID, \(\sum_{j} s_{jt} \leq 1\).
prices are prices \(p_{jt}\). These have some special properties and are always treated as endogenous.
demand_instruments0, demand_instruments1, and so on are numbered demand instruments. These represent only the excluded instruments. The exogenous regressors in \(X_1\) will be automatically added to the set of instruments.
We begin with two steps:
Load the product data which at a minimum consists of market_ids, shares, prices, and at least a single column of demand instruments, demand_instruments0.
Define a Formulation for the \(X_1\) (linear) demand model.
This and all other formulas are similar to R and patsy formulas.
It includes a constant by default. To exclude the constant, specify either a 0 or a 1.
To efficiently include fixed effects, use the absorb option and specify which categorical variables you would like to absorb.
Some model reduction may happen automatically. The constant will be excluded if you include fixed effects and some precautions are taken against collinearity. However, you will have to make sure that differentlynamed variables are not collinear.
The product_data argument of Problem should be a structured arraylike object with fields that store data. Product data can be a structured NumPy array, a pandas DataFrame, or other similar objects.
The product data contains market_ids, product_ids, firm_ids, shares, prices, a number of other IDs and product characteristics, and some precomputed excluded demand_instruments0, demand_instruments1, and so on. The product_ids will be incorporated as fixed effects.
For more information about the instruments and the example data as a whole, refer to the data module.
We can combine the Formulation and product_data to construct a Problem. We pass the Formulation first and the product_data second. We can also display the properties of the problem by typing its name.
The Problem.solve method always returns a ProblemResults class, which can be used to compute postestimation outputs. See the post estimation tutorial for more information.
[5]:
logit_results=problem.solve()logit_results
[5]:
Problem Results Summary:
=========================================
Computation GMM Objective Objective
Time Step Evaluations Value
   
00:00:01 2 2 +4.2E+05
=========================================
Beta Estimates (Robust SEs in Parentheses):
==========
prices

3.0E+01
(+1.0E+00)
==========
Now, we require that \(\epsilon_{jti} = \bar{\epsilon}_{h(j)ti} + (1  \rho) \bar{\epsilon}_{jti}\) is distributed IID with the Type I Extreme Value (Gumbel) distribution. As \(\rho \rightarrow 1\), all consumers stay within their group. As \(\rho \rightarrow 0\), this collapses to the IIA logit. Note that if we wanted, we could allow \(\rho\) to differ between groups with the notation \(\rho_{h(j)}\).
This gives us aggregate marketshares as the product of two logits, the within group logit and the across group logit:
where \(s_{jh(j)t} = s_{jt} / s_{h(j)t}\) and \(s_{h(j)t}\) is the share of group \(h\) in market \(t\). See Berry (1994) or Cardell (1997) for more information.
Again, the left hand side is data, though the \(\ln s_{jh(j)t}\) is clearly endogenous which means we must instrument for it. Rather than include \(\ln s_{jh(j)t}\) along with the linear components of utility, \(X_1\), whenever nesting_ids are included in product_data, \(\rho\) is treated as a nonlinear \(X_2\) parameter. This means that the linear component is given instead by
By including nesting_ids (another reserved name) as a field in product_data, we tell the package to estimate a nested logit model, and we donβt need to change any of the formulas. We show how to construct the category groupings in two different ways:
We put all products in a single nest (only the outside good in the other nest).
We put products into two nests (either mushy or nonmushy).
We also construct an additional instrument based on the number of products per nest. Typically this is useful as a source of exogenous variation in the within group share \(\ln s_{jh(j)t}\). However, in this example because the number of products per nest does not vary across markets, if we include product fixed effects, this instrument is irrelevant.
Weβll define a function that constructs the additional instrument and solves the nested logit problem. Weβll exclude product ID fixed effects, which are collinear with mushy, and weβll choose \(\rho = 0.7\) as the initial value at which the optimization routine will start.
Problem Results Summary:
==================================================================================
Computation GMM Optimization Objective Objective Gradient Hessian
Time Step Iterations Evaluations Value Infinity Norm Eigenvalue
      
00:00:11 2 3 8 +4.6E+05 +3.1E06 +2.4E+07
==================================================================================
Rho Estimates (Robust SEs in Parentheses):
==========
All Groups

+9.8E01
(+1.4E02)
==========
Beta Estimates (Robust SEs in Parentheses):
==========
prices

1.2E+00
(+4.0E01)
==========
When we inspect the Problem, the only changes from the plain logit model is the additional instrument that contributes to \(M_D\) and the inclusion of \(H\), the number of nesting categories.
[8]:
nl_results1.problem
[8]:
Dimensions:
===============================
T N F K1 MD H
     
94 2256 5 1 21 1
===============================
Formulations:
=====================================
Column Indices: 0
 
X1: Linear Characteristics prices
=====================================
Next, weβll solve the problem when there are two nests for mushy and nonmushy.
The package is designed to prevent the user from treating the within group share, \(\log s_{jh(j)t}\), as an exogenous variable. For example, if we were to compute a group_share variable and use the algebraic functionality of Formulation by including the expression log(shares/group_share) in our formula for \(X_1\), the package would raise an error because the package knows that shares should not be included in this
formulation.
To demonstrate why this is a bad idea, weβll override this feature by calculating \(\log s_{jh(j)t}\) and including it as an additional variable in \(X_1\). To do so, weβll first redefine our function for setting up and solving the nested logit problem.
In this tutorial, weβll use data from Nevo (2000) to solve the paperβs fake cereal problem. Locations of CSV files that contain the data are in the data module.
The random coefficients model extends the plain logit model by allowing for correlated tastes for different product characteristics. In this model (indirect) utility is given by
The main addition is that \(\beta_i = (\alpha_i, \beta_i^x)\) have individual specific subscripts \(i\).
Conditional on \(\beta_i\), the individual marketshares follow the same logit form as before. But now we must integrate over heterogeneous individuals to get the aggregate marketshares:
In general, this integral needs to be calculated numerically. This also means that we canβt directly linearize the model. It is common to reparametrize the model to separate the aspects of mean utility that all individuals agree on, \(\delta_{jt} = \alpha p_{jt} + x_{jt} \beta^x + \xi_{jt}\), from the individual specific heterogeneity, \(\mu_{jti}(\theta)\). This gives us
Given a guess of \(\theta\) we can solve the system of nonlinear equations for the vector \(\delta\) which equates observed and predicted marketshares \(s_{jt} = s_{jt}(\delta, \theta)\). Now we can perform a linear IV GMM regression of the form
Application of Random Coefficients Logit with the Fake Cereal DataΒΆ
To include random coefficients we need to add a Formulation for the nonlinear characteristics \(X_2\).
Just like in the logit case we have the same reserved field names in product_data:
market_ids are the unique market identifiers which we subscript \(t\).
shares specifies the marketshares which need to be between zero and one, and within a market ID, \(\sum_{j} s_{jt} < 1\).
prices are prices \(p_{jt}\). These have some special properties and are always treated as endogenous.
demand_instruments0, demand_instruments1, and so on are numbered demand instruments. These represent only the excluded instruments. The exogenous regressors in \(X_1\) (of which \(X_2\) is typically a subset) will be automatically added to the set of instruments.
We proceed with the following steps:
Load the productdata which at a minimum consists of market_ids, shares, prices, and at least a single column of demand instruments, demand_instruments0.
Define a Formulation for the \(X_1\) (linear) demand model.
This and all other formulas are similar to R and patsy formulas.
It includes a constant by default. To exclude the constant, specify either a 0 or a 1.
To efficiently include fixed effects, use the absorb option and specify which categorical variables you would like to absorb.
Some model reduction may happen automatically. The constant will be excluded if you include fixed effects and some precautions are taken against collinearity. However, you will have to make sure that differentlynamed variables are not collinear.
Define a Formulation for the \(X_2\) (nonlinear) demand model.
Include only the variables over which we want random coefficients.
Do not absorb or include fixed effects.
It will include a random coefficient on the constant (to capture inside good vs. outside good preference) unless you specify not to with a 0 or a 1.
Define an Integration configuration to solve the market share integral from several available options:
It is required to specify an initial guess of the nonlinear parameters. This serves two primary purposes: speeding up estimation and indicating to the solver through initial values of zero which parameters are restricted to be always zero.
Problem.solve takes an initial guess \(\Sigma_0\) of \(\Sigma\). It guarantees that \(\hat{\Sigma}\) (the estimated parameters) will have the same sparsity structure as \(\Sigma_0\). So any zero element of \(\Sigma\) is restricted to be zero in the solution \(\hat{\Sigma}\). For example, a popular restriction is that \(\Sigma\) is diagonal, this can be achieved by passing a diagonal matrix as
\(\Sigma_0\).
As is common with multivariate normal distributions, \(\Sigma\) is not estimated directly. Rather, its matrix square (Cholesky) root \(L\) is estimated where \(LL' = \Sigma\).
The product_data argument of Problem should be a structured arraylike object with fields that store data. Product data can be a structured NumPy array, a pandas DataFrame, or other similar objects.
The product data contains market_ids, product_ids, firm_ids, shares, prices, a number of other firm IDs and product characteristics, and some precomputed excluded demand_instruments0, demand_instruments1, and so on. The product_ids will be incorporated as fixed effects.
For more information about the instruments and the example data as a whole, refer to the data module.
Setting Up and Solving the Problem Without DemographicsΒΆ
Formulations, product data, and an integration configuration are collectively used to initialize a Problem. Once initialized, Problem.solve runs the estimation routine. The arguments to Problem.solve configure how estimation is performed. For example, optimization and iteration arguments configure the optimization and
iteration routines that are used by the outer and inner loops of estimation.
Weβll specify Formulation configurations for \(X_1\), the linear characteristics, and \(X_2\), the nonlinear characteristics.
The formulation for \(X_1\) consists of prices and fixed effects constructed from product_ids, which we will absorb using absorb argument of Formulation.
If we were interested in reporting estimates for each fixed effect, we could replace the formulation for \(X_1\) with Formulation('prices+C(product_ids)').
Because sugar, mushy, and the constant are collinear with product_ids, we can include them in \(X_2\) but not in \(X_1\).
We also need to specify an Integration configuration. We consider two alternatives:
Monte Carlo draws: we simulate 50 individuals from a random normal distribution.
Product rules: we construct nodes and weights according to a product rule that exactly integrates polynomials of degree \(5 \times 2  1 = 9\) or less.
Dimensions:
=============================================
T N F I K1 K2 MD ED
       
94 2256 5 58750 1 4 20 1
=============================================
Formulations:
===========================================================
Column Indices: 0 1 2 3
    
X1: Linear Characteristics prices
X2: Nonlinear Characteristics 1 prices sugar mushy
===========================================================
As an illustration of how to configure the optimization routine, weβll use a simpler, nondefault Optimization configuration that doesnβt support parameter bounds.
[8]:
bfgs=pyblp.Optimization('bfgs')bfgs
[8]:
Configured to optimize using the BFGS algorithm implemented in SciPy with analytic gradients and options {}.
We estimate three versions of the model:
An unrestricted covariance matrix for random tastes using Monte Carlo integration.
An unrestricted covariance matrix for random tastes using the product rule.
A restricted diagonal matrix for random tastes using Monte Carlo Integration.
Notice that the only thing that changes when we estimate the restricted covariance is our initial guess of \(\Sigma_0\). The lower diagonal in this initial guess is ignored because we are optimizing over the upper triangular (Cholesky) root of \(\Sigma\).
We see that all three models give similar estimates of the price coefficient \(\hat{\alpha} \approx 30\). Note a few of the estimated terms on the diagonal of \(\Sigma\) are negative. Since the diagonal consists of standard deviations, negative values are unrealistic. When using another optimization routine that supports bounds (like the default LBFGSB routine), these diagonal elements are by default bounded from below by zero.
To add demographic data we need to make two changes:
We need to load agent_data, which for this cereal problem contains precomputed Monte Carlo draws and demographics.
We need to add an agent_formulation to the model.
The agentdata has several reserved column names.
market_ids are the index that link the agentdata to the market_ids in productdata.
weights are the weights \(w_{ti}\) attached to each agent. In each market, these should sum to one so that \(\sum_i w_{ti} = 1\). It is often the case the \(w_{ti} = 1 / I_t\) where \(I_t\) is the number of agents in market \(t\), so that each agent gets equal weight. Other possibilities include quadrature nodes and weights.
nodes0, nodes1, and so on are the nodes at which the unobserved agent tastes \(\mu_{jti}\) are evaluated. The nodes should be labeled from \(0, \ldots, K_2  1\) where \(K_2\) is the number of random coefficients.
Other fields are the realizations of the demographics \(d\) themselves.
The agentformulation tells us which columns of demographic information to interact with \(X_2\).
[13]:
agent_formulation=pyblp.Formulation('0 + income + income_squared + age + child')agent_formulation
[13]:
income + income_squared + age + child
This tells us to include demographic realizations for income, income_squared, age, and the presence of children, child, but to ignore other possible demographics when interacting demographics with \(X_2\). We should also generally exclude the constant from the demographic formula.
Now we configure the Problem to include the agent_formulation and agent_data, which follow the product_formulations and product_data.
When we display the class, it lists the demographic interactions in the table of formulations and reports \(D = 4\), the dimension of the demographic draws.
Dimensions:
=================================================
T N F I K1 K2 D MD ED
        
94 2256 5 1880 1 4 4 20 1
=================================================
Formulations:
===================================================================
Column Indices: 0 1 2 3
    
X1: Linear Characteristics prices
X2: Nonlinear Characteristics 1 prices sugar mushy
d: Demographics income income_squared age child
===================================================================
The initialized problem can be solved with Problem.solve. Weβll use the same starting values as Nevo (2000). By passing a diagonal matrix as starting values for \(\Sigma\), weβre choosing to ignore covariance terms. Similarly, zeros in the starting values for \(\Pi\) mean that those parameters will be fixed at zero.
To replicate common estimates, weβll use the nondefault BFGS optimization routine, and weβll configure Problem.solve to use 1step GMM instead of 2step GMM.
Results are similar to those in the original paper with a GMM objective value of \(q(\hat{\theta}) = 4.56\) and a price coefficient of \(\hat{\alpha} = 62.7\).
Because the interactions between price, income, and income_squared are potentially collinear, we might worry that \(\hat{\Pi}_{21} = 588\) and \(\hat{\Pi}_{22} = 30.2\) are pushing the price coefficient in opposite directions. Both are large in magnitude but statistically insignficant. One way of dealing with this is to restrict \(\Pi_{22} = 0\).
There are two ways we can do this:
Change the initial \(\Pi_0\) values to make this term zero.
Change the agent formula to drop income_squared.
First, weβll change the initial \(\Pi_0\) values.
Now weβll drop both income_squared and the corresponding column in \(\Pi_0\).
[17]:
restricted_formulation=pyblp.Formulation('0 + income + age + child')deleted_pi=np.delete(initial_pi,1,axis=1)restricted_problem=pyblp.Problem(product_formulations,product_data,restricted_formulation,agent_data)restricted_problem.solve(initial_sigma,deleted_pi,optimization=bfgs,method='1s')
The parameter estimates and standard errors are identical for both approaches. Based on the number of fixed point iterations, there is some evidence that the solver took a slightly different path for each problem, but both restricted problems arrived at identical answers. The \(\hat{\Pi}_{12}\) interaction term is still insignificant.
The product_data contains market IDs, product IDs, firm IDs, shares, prices, a number of product characteristics, and some precomputed excluded instruments. The product IDs are called clustering_ids because they will be used to compute clustered standard errors. For more information about the instruments and the example data as a whole, refer to the data module.
The agent_data contains market IDs, integration weights \(w_{ti}\), integration nodes \(\nu_{ti}\), and demographics \(d_{ti}\). Here we use \(I_{t} = 200\) equally weighted Monte Carlo draws per market.
In nonexample problems, it is usually a better idea to use many more draws, or a more sophisticated Integration configuration such as sparse grid quadrature.
Unlike the fake cereal problem, we wonβt absorb any fixed effects in the automobile problem, so the linear part of demand \(X_1\) has more components. We also need to specify a formula for the random coefficients \(X_2\), including a random coefficient on the constant, which captures correlation among all inside goods.
The primary new addition to the model relative to the fake cereal problem is that we add a supply side formula for product characteristics that contribute to marginal costs, \(X_3\). The patsystyle formulas support functions of regressors such as the log function used below.
We stack the three product formulations in order: \(X_1\), \(X_2\), and \(X_3\).
[4]:
product_formulations=(pyblp.Formulation('1 + hpwt + air + mpd + space'),pyblp.Formulation('1 + prices + hpwt + air + mpd + space'),pyblp.Formulation('1 + log(hpwt) + air + log(mpg) + log(space) + trend'))product_formulations
The original specification for the automobile problem includes the term \(\log(y_i  p_j)\), in which \(y\) is income and \(p\) are prices. Instead of including this term, which gives rise to a host of numerical problems, weβll follow Berry, Levinsohn, and Pakes (1999) and use its firstorder linear approximation, \(p_j / y_i\).
The agent formulation for demographics, \(d\), includes a column of \(1 / y_i\) values, which weβll interact with \(p_j\). To do this, we will treat draws of \(y_i\) as demographic variables.
Choosing \(\Sigma\) and \(\Pi\) starting values, \(\Sigma_0\) and \(\Pi_0\).
Potentially choosing bounds for \(\Sigma\) and \(\Pi\).
Choosing a form for marginal costs, \(c_{jt}\): either a linear or loglinear functional form.
The decisions we will use are:
Use published estimates as our starting values in \(\Sigma_0\).
Interact the inverse of income, \(1 / y_i\), only with prices, and use the published estimate on \(\log(y_i  p_j)\) as our starting value for \(\alpha\) in \(\Pi_0\).
Bound \(\Sigma_0\) to be positive since it is a diagonal matrix where the diagonal consists of standard deviations.
Constrain the \(p_j / y_i\) coefficient to be negative. Specifically, weβll use a bound thatβs slightly smaller than zero because when the parameter is exactly zero, there are matrix inversion problems with estimating marginal costs.
When using a routine that supports bounds, Problem chooses some default bounds. These bounds are often very wide, so itβs usually a good idea to set your own more restrictive bounds.
Note that there are only 5 nonzeros on the diagonal of \(\Sigma\), which means that we only need 5 columns of integration nodes to integrate over these 5 dimensions of unobserved heterogeneity. Indeed, agent_data contains exactly 5 columns of nodes. If we were to ignore the \(\log(y_i  p_j)\) term (by not configuring \(\Pi\)) and include a term on prices in \(\Sigma\) instead, we would have needed 6 columns of integration nodes in our agent_data.
A linear marginal cost specification is the default setting, so weβll need to use the costs_type argument of Problem.solve to employ the loglinear specification used by Berry, Levinsohn, and Pakes (1995). A downside of this specification is that nonpositive estimated marginal costs can create problems for the optimization routine when computing \(\log c(\hat{\theta})\).
Weβll use the costs_bounds argument to bound marginal costs from below by a small number.
Finally, as in the original paper, weβll use W_type and se_type to cluster by product IDs, which were specified as clustering_ids in product_data.
Problem Results Summary:
==================================================================================================================================
Smallest Largest Clipped
Computation GMM Optimization Objective Fixed Point Contraction Objective Gradient Hessian Hessian Marginal
Time Step Iterations Evaluations Iterations Evaluations Value Infinity Norm Eigenvalue Eigenvalue Costs
          
00:31:03 2 39 49 72322 217903 +2.1E+05 +2.8E+03 +1.8E+02 +1.8E+04 0
==================================================================================================================================
Nonlinear Coefficient Estimates (Robust SEs Adjusted for 999 Clusters in Parentheses):
=====================================================================================================
Sigma: 1 prices hpwt air mpd space  Pi: 1/income
         
1 +1.2E+00 +0.0E+00 +0.0E+00 +0.0E+00 +0.0E+00 +0.0E+00  1 +0.0E+00
(+2.1E+00) 

prices +0.0E+00 +0.0E+00 +0.0E+00 +0.0E+00 +0.0E+00  prices 1.2E+01
 (+4.8E+00)

hpwt +0.0E+00 +0.0E+00 +0.0E+00 +0.0E+00  hpwt +0.0E+00
(+1.0E+01) 

air +0.0E+00 +0.0E+00 +0.0E+00  air +0.0E+00
(+2.6E+00) 

mpd +2.6E+00 +0.0E+00  mpd +0.0E+00
(+1.4E+00) 

space +1.1E+00  space +0.0E+00
(+1.9E+00) 
=====================================================================================================
Beta Estimates (Robust SEs Adjusted for 999 Clusters in Parentheses):
==========================================================
1 hpwt air mpd space
    
7.6E+00 +4.8E+00 +2.6E+00 2.2E+00 +1.8E+00
(+1.2E+00) (+2.7E+00) (+1.4E+00) (+1.5E+00) (+1.3E+00)
==========================================================
Gamma Estimates (Robust SEs Adjusted for 999 Clusters in Parentheses):
======================================================================
1 log(hpwt) air log(mpg) log(space) trend
     
+2.3E+00 +6.4E01 +8.8E01 3.9E01 +2.8E01 +1.4E02
(+1.4E01) (+1.2E01) (+7.1E02) (+1.2E01) (+1.2E01) (+3.5E03)
======================================================================
There are some discrepancies between our results and the original paper. The instruments we constructed to are meant to mimic the original instruments, which we were unable to reconstruct perfectly. We also use different agent data and the firstorder linear approximation to the \(\ln(y_i  p_j)\) term.
As in the fake cereal tutorial, weβll first solve the fake cereal problem from Nevo (2000). We load the fake data and estimate the model as in the previous tutorial. We output the setup of the model to confirm we have correctly configured the Problem
[2]:
product_data=pd.read_csv(pyblp.data.NEVO_PRODUCTS_LOCATION)agent_data=pd.read_csv(pyblp.data.NEVO_AGENTS_LOCATION)product_formulations=(pyblp.Formulation('0 + prices',absorb='C(product_ids)'),pyblp.Formulation('1 + prices + sugar + mushy'))agent_formulation=pyblp.Formulation('0 + income + income_squared + age + child')problem=pyblp.Problem(product_formulations,product_data,agent_formulation,agent_data)problem
[2]:
Dimensions:
=================================================
T N F I K1 K2 D MD ED
        
94 2256 5 1880 1 4 4 20 1
=================================================
Formulations:
===================================================================
Column Indices: 0 1 2 3
    
X1: Linear Characteristics prices
X2: Nonlinear Characteristics 1 prices sugar mushy
d: Demographics income income_squared age child
===================================================================
Weβll solve the problem in the same way as before. The Problem.solve method returns a ProblemResults class, which displays basic estimation results. The results that are displayed are simply formatted information extracted from various class attributes such as ProblemResults.sigma and
ProblemResults.sigma_se.
Postestimation outputs are computed for each market and stacked. Weβll use matplotlib functions to display the matrices associated with a single market.
The diagonal of the first image consists of own elasticities and the diagonal of the second image consists of diversion ratios to the outside good. As one might expect, own price elasticities are large and negative while crossprice elasticities are positive but much smaller.
An alternative to summarizing full elasticity matrices is to use ProblemResults.compute_aggregate_elasticities to estimate aggregate elasticities of demand, \(E\), in each market, which reflect the change in total sales under a proportional sales tax of some factor.
Since demand for an entire product category is generally less elastic than the average elasticity of individual products, mean own elasticities are generally larger in magnitude than aggregate elasticities.
[9]:
plt.hist([means.flatten(),aggregates.flatten()],color=['red','blue'],bins=50);plt.legend(['Mean Own Elasticities','Aggregate Elasticities']);
To compute marginal costs, \(c\), the product_data passed to Problem must have had a firm_ids field. Since we included firm IDs when configuring the problem, we can use ProblemResults.compute_costs.
Other methods that compute supplyside outputs often compute marginal costs themselves. For example, ProblemResults.compute_markups will compute marginal costs when estimating markups, \(\mathscr{M}\), but computation can be sped up if we just use our precomputed values.
Before computing postmerger outputs, weβll supplement our premerger markups with some other outputs. Weβll compute HerfindahlHirschman Indices, \(\text{HHI}\), with ProblemResults.compute_hhi; populationnormalized gross expected profits, \(\pi\), with ProblemResults.compute_profits; and
populationnormalized consumer surpluses, \(\text{CS}\), with ProblemResults.compute_consumer_surpluses.
We can use ProblemResults.compute_approximate_prices or ProblemResults.compute_prices to estimate postmerger prices. The first method, which is discussed, for example, in Nevo (1997), assumes that shares and their price derivatives are unaffected by
the merger. The second method does not make these assumptions and iterates over the \(\zeta\)markup equation from Morrow and Skerlos (2011) to solve the full system of \(J_t\) equations and \(J_t\) unknowns in each market \(t\). Weβll use the latter, since it is fast enough for this example problem.
If the problem was configured with more than two columns of firm IDs, we could estimate postmerger prices for the other mergers with the firms_index argument, which is by default 1.
Postestimation outputs can be informative, but they donβt mean much without a sense sampletosample variability. One way to estimate confidence intervals for postestimation outputs is with a standard bootstrap procedure:
Construct a large number of bootstrap samples by sampling with replacement from the original product data.
Initialize and solve a Problem for each bootstrap sample.
Compute the desired postestimation output for each bootstrapped ProblemResults and from the resulting empirical distribution, construct boostrap confidence intervals.
Although appealing because of its simplicity, the computational resources required for this procedure are often prohibatively expensive. Furthermore, human oversight of the optimization routine is often required to determine whether the routine ran into any problems and if it successfully converged. Human oversight of estimation for each bootstrapped problem is usually not feasible.
A more reasonable alternative is a parametric bootstrap procedure:
Construct a large number of draws from the estimated joint distribution of parameters.
Compute the implied mean utility, \(\delta\), and shares, \(s\), for each draw. If a supply side was estimated, also computed the implied marginal costs, \(c\), and prices, \(p\).
Compute the desired postestimation output under each of these parametric bootstrap samples. Again, from the resulting empirical distribution, construct boostrap confidence intervals.
Compared to the standard bootstrap procedure, the parametric bootstrap requires far fewer computational resources, and is simple enough to not require human oversight of each bootstrap iteration. The primary complication to this procedure is that when supply is estimated, equilibrium prices and shares need to be computed for each parametric bootstrap sample by iterating over the \(\zeta\)markup equation from Morrow and Skerlos (2011).
Although nontrivial, this fixed point iteration problem is much less demanding than the full optimization routine required to solve the BLP problem from the start.
An empirical distribution of results computed according to this parametric bootstrap procedure can be created with the ProblemResults.bootstrap method, which returns a BootstrappedResults class that can be used just like ProblemResults to compute various postestimation outputs.
The difference is that BootstrappedResults methods return arrays with an extra first dimension, along which bootstrapped results are stacked.
Weβll construct 90% parametric bootstrap confidence intervals for estimated mean own elasticities in each market of the fake cereal problem. Usually, bootstrapped confidence intervals should be based on thousands of draws, but weβll only use a few for the sake of speed in this example.
Bootstrapped Problem Results Summary:
======================
Computation Bootstrap
Time Draws
 
00:00:24 100
======================
[21]:
bounds=np.percentile(bootstrapped_results.extract_diagonal_means(bootstrapped_results.compute_elasticities()),q=[10,90],axis=0)table=pd.DataFrame(index=problem.unique_market_ids,data={'Lower Bound':bounds[0].flatten(),'Mean Own Elasticity':aggregates.flatten(),'Upper Bound':bounds[1].flatten()})table.round(2).head()
Given a consistent estimate of \(\theta\), we may want to compute the optimal instruments of Chamberlain (1987) and use them to resolve the problem. Optimal instruments have been shown, for example, by Reynaert and Verboven (2014), to reduce bias, improve efficiency, and enhance stability of BLP estimates.
The ProblemResults.compute_optimal_instruments method computes the expected Jacobians that comprise the optimal instruments by integrating over the density of \(\xi\) (and \(\omega\) if a supply side was estimated). By default, the method approximates this integral by averaging over the Jacobian realizations computed under draws from the asymptotic normal distribution
of the error terms. Since this process is computationally expensive and often doesnβt make much of a difference, weβll use method='approximate' in this example to simply evaluate the Jacobians at the expected value of \(\xi\), zero.
Optimal Instrument Results Summary:
=======================
Computation Error Term
Time Draws
 
00:00:01 1
=======================
We can use the OptimalInstrumentResults.to_problem method to recreate the fake cereal problem with the estimated optimal excluded instruments.
Before configuring and solving a problem with real data, it may be a good idea to perform Monte Carlo analysis on simulated data to verify that it is possible to accurately estimate model parameters. For example, before configuring and solving the example problems in the prior tutorials, it may have been a good idea to simulate data according to the assumed models of supply and demand. During such Monte Carlo anaysis, the data would only be used to determine sample sizes and perhaps to choose
reasonable true parameters.
Simulations are configured with the Simulation class, which requires many of the same inputs as Problem. The two main differences are:
Variables in formulations that cannot be loaded from product_data or agent_data will be drawn from independent uniform distributions.
True parameters and the distribution of unobserved product characteristics are specified.
First, weβll use build_id_data to build market and firm IDs for a model in which there are \(T = 50\) markets, and in each market \(t\), a total of \(J_t = 20\) products produced by \(F = 10\) firms.
[2]:
id_data=pyblp.build_id_data(T=50,J=20,F=10)
Next, weβll create an Integration configuration to build agent data according to a GaussHermite product rule that exactly integrates polynomials of degree \(2 \times 9  1 = 17\) or less.
Configured to construct nodes and weights according to the level9 GaussHermite product rule.
Weβll then pass these data to Simulation. Weβll use Formulation configurations to create an \(X_1\) that consists of a constant, prices, and an exogenous characteristic; an \(X_2\) that consists only of the same exogenous characteristic; and an \(X_3\) that consists of the common exogenous characteristic and a costshifter.
The excluded instruments in Simulation.product_data include basic instruments computed with build_blp_instruments that are functions of all exogenous numerical variables in the problem. In this example, excluded demandside instruments are the costshifter z and traditional BLP instruments constructed from x. Excluded supplyside instruments
are traditional BLP instruments constructed from x and z.
The Simulation can be further configured with other arguments that determine how unobserved product characteristics are simulated and how marginal costs are specified.
Since at this stage prices and shares are all zeros, we still need to solve the simulation with Simulation.solve. This method computes synthetic prices and shares. Just like ProblemResults.compute_prices, to do so it iterates over the \(\zeta\)markup equation from Morrow and Skerlos
(2011).
Simulation Results Summary:
=====================================
Computation Fixed Point Contraction
Time Iterations Evaluations
  
00:00:01 697 697
=====================================
Now, we can try to recover the true parameters by creating and solving a Problem. By default, the convenience method SimulationResults.to_problem uses the same formulations and unobserved agent data as the simulation, so estimation is relatively easy. However, weβll choose starting values that are half the true parameters so that the optimization routine
has to do some work. Note that since weβre jointly estimating the supply side, we need to provide an initial value for the linear coefficient on prices because this parameter cannot be concentrated out of the problem (unlike linear coefficients on exogenous characteristics).
[6]:
problem=simulation_results.to_problem()problem
[6]:
Dimensions:
=================================================
T N F I K1 K2 K3 MD MS
        
50 1000 10 450 3 1 2 5 6
=================================================
Formulations:
===================================================
Column Indices: 0 1 2
   
X1: Linear Characteristics 1 prices x
X2: Nonlinear Characteristics x
X3: Cost Characteristics x z
===================================================
In addition to checking that the configuration for a model based on actual data makes sense, the Simulation class can also be a helpful tool for better understanding under what general conditions BLP models can be accurately estimated. Simulations are also used extensively in pyblpβs test suite.
The majority of the package consists of classes, which compartmentalize different aspects of the BLP model. There are some convenience functions as well.
Various components of the package require configurations for how to approximate integrals, solve fixed point problems, and solve optimimzation problems. Such configurations are specified with the following classes.
Approximate equilibrium prices after firm or cost changes, \(p^*\), under the assumption that shares and their price derivatives are unaffected by such changes.
A parametric bootstrap can be used, for example, to compute standard errors for the above postestimation outputs. The following method returns a results class with all of the above methods, which returns a distribution of postestimation outputs corresponding to different bootstrapped samples.
Product and agent data that are passed or constructed by Problem and Simulation are structured internally into classes with field names that more closely resemble BLP notation. Although these structured data classes are not directly constructable, they can be accessed with Problem and Simulation class attributes. It can be helpful to compare these structured data classes with the data or configurations used to create them.
In addition to classes and functions, there are also two modules that can be used to configure global package options and locate example data that comes with the package.
Failed to invert a Jacobian of shares with respect to \(\xi\) when computing the Jacobian of \(\xi\) (equivalently, of \(\delta\)) with respect to \(\theta\).
This page contains a full list of references cited in the documentation, including the original work of Berry, Levinsohn, and Pakes (1995). If you use pyblp in your research, we ask that you also cite the below Conlon and Gortmaker (2019), which describes the advances implemented in the package.
Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the βSoftwareβ), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED βAS ISβ, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
Testing is done with the tox automation tool, which runs a pytestbacked test suite in the tests/ directory. This FAQ contains some useful information about how to use tox on Windows.
In addition to the installation requirements for the package itself, running tests and building documentation requires additional packages specified by the tests and docs extras in setup.py, along with any other explicitly specified deps in tox.ini.
The full suite of tests also requires installation of the following software:
Artleys Knitro version 10.3 or newer: testing optimization routines.
MATLAB: comparing sparse grids with those created by the function nwspgr created by Florian Heiss and Viktor Winschel, which must be included in a directory on the MATLAB path.
If software is not installed, its associated tests will be skipped. Additionally, some tests that require support for extended precision will be skipped if on the platform running the tests, numpy.longdouble has the same precision as numpy.float64. This tends to be the case on Windows.
Defined in tox.ini are environments that test the package under different python versions, check types, enforce style guidelines, verify the integrity of the documentation, and release the package. The following command can be run in the toplevel pyblp directory to run all testing environments:
tox
You can choose to run only one environment, such as the one that builds the documentation, with the e flag:
Fixtures, which are defined in tests.conftest, configure the testing environment and simulate problems according to a range of specifications.
Most BLPspecific tests in tests.test_blp verify properties about results obtained by solving the simulated problems under various parameterizations. Examples include:
Reasonable formulations of problems should give rise to estimated parameters that are close to their true values.
Cosmetic changes such as the number of processes should not change estimates.
Postestimation outputs should satisfy certain properties.
Optimization routines should behave as expected.
Derivatives computed with finite differences should approach analytic derivatives.
Tests of generic utilities in tests.test_formulation, tests.test_integration, tests.test_iteration, and tests.test_optimization verify that matrix formulation, integral approximation, fixed point iteration, and nonlinear optimization all work as expected. Example include:
Nonlinear formulas give rise to expected matrices and derivatives.
GaussHermite integrals are better approximated with quadrature based on GaussHermite rules than with Monte Carlo integration.
To solve a fixed point iteration problem for which it was developed, SQUAREM requires fewer fixed point evaluations than does simple iteration.
All optimization routines manage to solve a wellknown optimization problem under different parameterizations.
"
],
"text/plain": [
" market_ids weights nodes0 nodes1 nodes2 nodes3 nodes4 \\\n",
"0 1971 0.005 0.548814 0.457760 0.564690 0.395537 0.392173 \n",
"1 1971 0.005 0.715189 0.376918 0.839746 0.844017 0.041157 \n",
"2 1971 0.005 0.602763 0.702335 0.376884 0.150442 0.923301 \n",
"3 1971 0.005 0.544883 0.207324 0.499676 0.306309 0.406235 \n",
"4 1971 0.005 0.423655 0.074280 0.081302 0.094570 0.944282 \n",
"\n",
" income \n",
"0 9.728478 \n",
"1 7.908957 \n",
"2 11.079404 \n",
"3 17.641671 \n",
"4 12.423995 "
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"agent_data = pd.read_csv(pyblp.data.BLP_AGENTS_LOCATION)\n",
"agent_data.head()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Setting up the Problem\n",
"\n",
"Unlike the fake cereal problem, we won't absorb any fixed effects in the automobile problem, so the linear part of demand $X_1$ has more components. We also need to specify a formula for the random coefficients $X_2$, including a random coefficient on the constant, which captures correlation among all inside goods.\n",
"\n",
"The primary new addition to the model relative to the fake cereal problem is that we add a supply side formula for product characteristics that contribute to marginal costs, $X_3$. The [patsy](https://patsy.readthedocs.io/en/stable/)style formulas support functions of regressors such as the `log` function used below.\n",
"\n",
"We stack the three product formulations in order: $X_1$, $X_2$, and $X_3$."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(1 + hpwt + air + mpd + space,\n",
" 1 + prices + hpwt + air + mpd + space,\n",
" 1 + log(hpwt) + air + log(mpg) + log(space) + trend)"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"product_formulations = (\n",
" pyblp.Formulation('1 + hpwt + air + mpd + space'),\n",
" pyblp.Formulation('1 + prices + hpwt + air + mpd + space'),\n",
" pyblp.Formulation('1 + log(hpwt) + air + log(mpg) + log(space) + trend')\n",
")\n",
"product_formulations"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The original specification for the automobile problem includes the term $\\log(y_i  p_j)$, in which $y$ is income and $p$ are prices. Instead of including this term, which gives rise to a host of numerical problems, we'll follow [Berry, Levinsohn, and Pakes (1999)](https://pyblp.readthedocs.io/en/stable/references.html#berrylevinsohnandpakes1999) and use its firstorder linear approximation, $p_j / y_i$. \n",
"\n",
"The agent formulation for demographics, $d$, includes a column of $1 / y_i$ values, which we'll interact with $p_j$. To do this, we will treat draws of $y_i$ as demographic variables."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"I(1 / income)"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"agent_formulation = pyblp.Formulation('0 + I(1 / income)')\n",
"agent_formulation"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As in the cereal example, the [`Problem`](https://pyblp.readthedocs.io/en/stable/_api/pyblp.Problem.html#pyblp.Problem) can be constructed by combining the `product_formulations`, `product_data`, `agent_formulation`, and `agent_data`. "
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Dimensions:\n",
"=======================================================\n",
" T N F I K1 K2 K3 D MD MS \n",
"         \n",
"20 2217 26 4000 5 6 6 1 11 12 \n",
"=======================================================\n",
"\n",
"Formulations:\n",
"======================================================================================\n",
" Column Indices: 0 1 2 3 4 5 \n",
"      \n",
" X1: Linear Characteristics 1 hpwt air mpd space \n",
"X2: Nonlinear Characteristics 1 prices hpwt air mpd space\n",
" X3: Cost Characteristics 1 log(hpwt) air log(mpg) log(space) trend\n",
" d: Demographics 1/income \n",
"======================================================================================"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"problem = pyblp.Problem(product_formulations, product_data, agent_formulation, agent_data)\n",
"problem"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The problem outputs a table of dimensions:\n",
"\n",
" $T$ denotes the number of markets.\n",
" $N$ is the length of the dataset (the number of products across all markets).\n",
" $F$ denotes the number of firms.\n",
" $I = \\sum_t I_t$ is the total number of agents across all markets (200 draws per market times 20 markets).\n",
" $K_1$ is the number of linear demand characteristics.\n",
" $K_2$ is the number of nonlinear demand characteristics.\n",
" $K_3$ is the number of linear supply characteristics.\n",
" $D$ is the number of demographic variables.\n",
" $M_D$ is the number of demand instruments, including exogenous regressors.\n",
" $M_S$ is the number of supply instruments, including exogenous regressors.\n",
"\n",
"The formulations table describes all four formulas for linear characteristics, nonlinear characteristics, cost characteristics, and demographics."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Solving the Problem"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The only remaining decisions are:\n",
"\n",
" Choosing $\\Sigma$ and $\\Pi$ starting values, $\\Sigma_0$ and $\\Pi_0$.\n",
" Potentially choosing bounds for $\\Sigma$ and $\\Pi$.\n",
" Choosing a form for marginal costs, $c_{jt}$: either a linear or loglinear functional form.\n",
"\n",
"The decisions we will use are:\n",
"\n",
" Use published estimates as our starting values in $\\Sigma_0$.\n",
" Interact the inverse of income, $1 / y_i$, only with prices, and use the published estimate on $\\log(y_i  p_j)$ as our starting value for $\\alpha$ in $\\Pi_0$.\n",
" Bound $\\Sigma_0$ to be positive since it is a diagonal matrix where the diagonal consists of standard deviations.\n",
" Constrain the $p_j / y_i$ coefficient to be negative. Specifically, we'll use a bound that's slightly smaller than zero because when the parameter is exactly zero, there are matrix inversion problems with estimating marginal costs.\n",
"\n",
"When using a routine that supports bounds, [`Problem`](https://pyblp.readthedocs.io/en/stable/_api/pyblp.Problem.html#pyblp.Problem) chooses some default bounds. These bounds are often very wide, so it's usually a good idea to set your own more restrictive bounds."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"initial_sigma = np.diag([3.612, 0, 4.628, 1.818, 1.050, 2.056])\n",
"initial_pi = np.c_[[0, 43.501, 0, 0, 0, 0]]\n",
"sigma_bounds = (\n",
" np.zeros_like(initial_sigma),\n",
" np.diag([100, 0, 100, 100, 100, 100])\n",
")\n",
"pi_bounds = (\n",
" np.c_[[0, 100, 0, 0, 0, 0]],\n",
" np.c_[[0, 0.1, 0, 0, 0, 0]]\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Note that there are only 5 nonzeros on the diagonal of $\\Sigma$, which means that we only need 5 columns of integration nodes to integrate over these 5 dimensions of unobserved heterogeneity. Indeed, `agent_data` contains exactly 5 columns of nodes. If we were to ignore the $\\log(y_i  p_j)$ term (by not configuring $\\Pi$) and include a term on prices in $\\Sigma$ instead, we would have needed 6 columns of integration nodes in our `agent_data`.\n",
"\n",
"A linear marginal cost specification is the default setting, so we'll need to use the `costs_type` argument of [`Problem.solve`](https://pyblp.readthedocs.io/en/stable/_api/pyblp.Problem.solve.html#pyblp.Problem.solve) to employ the loglinear specification used by [Berry, Levinsohn, and Pakes (1995)](https://pyblp.readthedocs.io/en/stable/references.html#berrylevinsohnandpakes1995). A downside of this specification is that nonpositive estimated marginal costs can create problems for the optimization routine when computing $\\log c(\\hat{\\theta})$. We'll use the `costs_bounds` argument to bound marginal costs from below by a small number. \n",
"\n",
"Finally, as in the original paper, we'll use `W_type` and `se_type` to cluster by product IDs, which were specified as `clustering_ids` in `product_data`."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Problem Results Summary:\n",
"==================================================================================================================================\n",
" Smallest Largest Clipped \n",
"Computation GMM Optimization Objective Fixed Point Contraction Objective Gradient Hessian Hessian Marginal\n",
" Time Step Iterations Evaluations Iterations Evaluations Value Infinity Norm Eigenvalue Eigenvalue Costs \n",
"          \n",
" 00:31:03 2 39 49 72322 217903 +2.1E+05 +2.8E+03 +1.8E+02 +1.8E+04 0 \n",
"==================================================================================================================================\n",
"\n",
"Nonlinear Coefficient Estimates (Robust SEs Adjusted for 999 Clusters in Parentheses):\n",
"=====================================================================================================\n",
"Sigma: 1 prices hpwt air mpd space  Pi: 1/income \n",
"         \n",
" 1 +1.2E+00 +0.0E+00 +0.0E+00 +0.0E+00 +0.0E+00 +0.0E+00  1 +0.0E+00 \n",
" (+2.1E+00)  \n",
"  \n",
"prices +0.0E+00 +0.0E+00 +0.0E+00 +0.0E+00 +0.0E+00  prices 1.2E+01 \n",
"  (+4.8E+00)\n",
"  \n",
" hpwt +0.0E+00 +0.0E+00 +0.0E+00 +0.0E+00  hpwt +0.0E+00 \n",
" (+1.0E+01)  \n",
"  \n",
" air +0.0E+00 +0.0E+00 +0.0E+00  air +0.0E+00 \n",
" (+2.6E+00)  \n",
"  \n",
" mpd +2.6E+00 +0.0E+00  mpd +0.0E+00 \n",
" (+1.4E+00)  \n",
"  \n",
"space +1.1E+00  space +0.0E+00 \n",
" (+1.9E+00)  \n",
"=====================================================================================================\n",
"\n",
"Beta Estimates (Robust SEs Adjusted for 999 Clusters in Parentheses):\n",
"==========================================================\n",
" 1 hpwt air mpd space \n",
"    \n",
" 7.6E+00 +4.8E+00 +2.6E+00 2.2E+00 +1.8E+00 \n",
"(+1.2E+00) (+2.7E+00) (+1.4E+00) (+1.5E+00) (+1.3E+00)\n",
"==========================================================\n",
"\n",
"Gamma Estimates (Robust SEs Adjusted for 999 Clusters in Parentheses):\n",
"======================================================================\n",
" 1 log(hpwt) air log(mpg) log(space) trend \n",
"     \n",
" +2.3E+00 +6.4E01 +8.8E01 3.9E01 +2.8E01 +1.4E02 \n",
"(+1.4E01) (+1.2E01) (+7.1E02) (+1.2E01) (+1.2E01) (+3.5E03)\n",
"======================================================================"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"results = problem.solve(\n",
" initial_sigma,\n",
" initial_pi,\n",
" sigma_bounds=sigma_bounds,\n",
" pi_bounds=pi_bounds,\n",
" costs_type='log',\n",
" costs_bounds=(0.001, None),\n",
" W_type='clustered',\n",
" se_type='clustered'\n",
")\n",
"results"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"There are some discrepancies between our results and the original paper. The instruments we constructed to are meant to mimic the original instruments, which we were unable to reconstruct perfectly. We also use different agent data and the firstorder linear approximation to the $\\ln(y_i  p_j)$ term."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/xpython",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.6"
}
},
"nbformat": 4,
"nbformat_minor": 2
}PK+N΅€[wϊ
ϊ
Jpyblpstable/_downloads/92fc72540c43fa666776d73d9225b023/formulation.ipynb{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Formulation Example"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"'0.7.0'"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import pyblp\n",
"\n",
"pyblp.__version__"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this example, we'll design a matrix without an intercept, but with both prices and another numeric size variable."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"prices + size"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"formulation = pyblp.Formulation('0 + prices + size')\n",
"formulation"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next, we'll design a second matrix with an intercept, with first and seconddegree size terms, with categorical product IDs and years, and with the interaction of the last two. The first formulation will include include the fixed effects as indicator variables, and the second will absorb them."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1 + size + I(size ** 2) + C(product) + C(year) + C(product):C(year)"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"formulation1 = pyblp.Formulation('size + I(size ** 2) + C(product) * C(year)')\n",
"formulation1"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"size + I(size ** 2) + Absorb[C(product)] + Absorb[C(year)] + Absorb[C(product):C(year)]"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"formulation2 = pyblp.Formulation('size + I(size ** 2)', absorb='C(product) * C(year)')\n",
"formulation2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally, we'll design a third matrix with an intercept and with a yearly trend interacted with the natural logarithm of income and categorical education. Absorption of continuous variables is not supported, so we need to use dummy variables."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1 + year:log(income) + year:C(education)"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"formulation = pyblp.Formulation('year:(log(income) + C(education))')\n",
"formulation"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/xpython",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.6"
}
},
"nbformat": 4,
"nbformat_minor": 2
}PK+NkΙ**Tpyblpstable/_downloads/eb4d6b8d4399a95225e8c60fca4edd43/build_blp_instruments.ipynb{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Building BLP Instruments Example"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"'0.7.0'"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import pyblp\n",
"import pandas as pd\n",
"\n",
"pyblp.__version__"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this example, we'll load the automobile product data from [Berry, Levinsohn, and Pakes (1995)](https://pyblp.readthedocs.io/en/stable/references.html#berrylevinsohnandpakes1995) and build some very simple excluded demandside instruments for the problem. These instruments are different from the prebuilt ones included in the automobile product data file."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1 + hpwt + air + mpd + space"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"formulation = pyblp.Formulation('1 + hpwt + air + mpd + space')\n",
"formulation"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plt.colorbar(plt.matshow(diversions[single_market]));"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The diagonal of the first image consists of own elasticities and the diagonal of the second image consists of diversion ratios to the outside good. As one might expect, own price elasticities are large and negative while crossprice elasticities are positive but much smaller.\n",
"\n",
"Elasticities and diversion ratios can be computed with respect to variables other than `prices` with the `name` argument of [`ProblemResults.compute_elasticities`](https://pyblp.readthedocs.io/en/stable/_api/pyblp.ProblemResults.compute_elasticities.html#pyblp.ProblemResults.compute_elasticities) and [`ProblemResults.compute_diversion_ratios`](https://pyblp.readthedocs.io/en/stable/_api/pyblp.ProblemResults.compute_diversion_ratios.html#pyblp.ProblemResults.compute_diversion_ratios). Additionally, [`ProblemResults.compute_long_run_diversion_ratios`](https://pyblp.readthedocs.io/en/stable/_api/pyblp.ProblemResults.compute_long_run_diversion_ratios.html#pyblp.ProblemResults.compute_long_run_diversion_ratios) can be used to used to understand substitution when products are eliminated from the choice set.\n",
"\n",
"The convenience methods [`ProblemResults.extract_diagonals`](https://pyblp.readthedocs.io/en/stable/_api/pyblp.ProblemResults.extract_diagonals.html#pyblp.ProblemResults.extract_diagonals) and [`ProblemResults.extract_diagonal_means`](https://pyblp.readthedocs.io/en/stable/_api/pyblp.ProblemResults.extract_diagonal_means.html#pyblp.ProblemResults.extract_diagonal_means) can be used to extract information about own elasticities of demand from elasticity matrices."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"means = results.extract_diagonal_means(elasticities)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"An alternative to summarizing full elasticity matrices is to use [`ProblemResults.compute_aggregate_elasticities`](https://pyblp.readthedocs.io/en/stable/_api/pyblp.ProblemResults.compute_aggregate_elasticities.html#pyblp.ProblemResults.compute_aggregate_elasticities) to estimate aggregate elasticities of demand, $E$, in each market, which reflect the change in total sales under a proportional sales tax of some factor."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"aggregates = results.compute_aggregate_elasticities(factor=0.1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Since demand for an entire product category is generally less elastic than the average elasticity of individual products, mean own elasticities are generally larger in magnitude than aggregate elasticities."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAHnZJREFUeJzt3Xt0FPX9//HnG0wFEQUlKhgt2C+IXEKAgHITqBUUEBXwHAGrLWL0WKvSL1TwBm21p176hUq9UbWg1XjBevSnfH+CAlXU1gbLHQWlsUT4QYySEi4C4f37YyfbBHbJZrPZhPH1OGdOZmY/O/PeIfvKMJfPmLsjIiJHv0b1XYCIiKSGAl1EJCQU6CIiIaFAFxEJCQW6iEhIKNBFREJCgS4iEhIKdBGRkFCgi4iExDHpXFmrVq28bdu26VyliMhRb/ny5V+6e2Z17dIa6G3btqWgoCCdqxQROeqZ2eeJtEv4kIuZNTazf5jZ68F0OzP7m5ltNLMXzOw7yRYrIiK1V5Nj6LcA6ytN3wfMdPf2wNfAtaksTEREaiahQDezLGA48EQwbcD3gflBk3nAZXVRoIiIJCbRY+izgJ8DzYPpk4Ed7n4gmC4CTo/1RjPLA/IAzjzzzMNe379/P0VFRezdu7cGZYtENGnShKysLDIyMuq7FJF6V22gm9kIYLu7LzezQRWzYzSN2bG6u88B5gDk5uYe1qaoqIjmzZvTtm1bIjv+Iolxd0pKSigqKqJdu3b1XY5IvUtkD70fMNLMhgFNgBOI7LG3MLNjgr30LGBLMgXs3btXYS5JMTNOPvlkiouL67sUkQah2mPo7j7N3bPcvS1wJbDY3ccDS4AxQbNrgFeTLUJhLsnS747If9TmTtHbgJ+Z2adEjqk/mZqSREQkGTUKdHdf6u4jgvFN7t7b3f/L3a9w929SUpFZaoeEVmn88Ic/jE4fOHCAzMxMRowYkZKPFI+7c88999C+fXs6dOjA4MGDWbt2bUrXsXTpUk488URycnKiw1tvvQXA8ccfn9QyZ82axe7du6PTw4YNY8eOHXHbP/bYYzz99NMAzJ07ly1b/nN0buLEiaxbty6pOkSkqrTeKdpQNWvWjDVr1rBnzx6aNm3KokWLOP30mBftpNTDDz/M+++/z8qVKznuuONYuHAhI0eOZO3atTRp0iRl6xkwYACvv/56ypY3a9YsrrrqKo477jgAFixYcMT2N9xwQ3R87ty5dOnShTZt2gDwxBNPpKwukbpWsY/oMS8BqX/qnCtw8cUX88YbbwCQn5/P2LFjo6/t2rWLCRMm0KtXL7p3786rr0ZOFxQWFjJgwAB69OhBjx49eP/994HIXvGgQYMYM2YMHTt2ZPz48XiM34D77ruP2bNnR4NxyJAh9O3bl2effZYXX3yRn/3sZwD87ne/46yzzgLgs88+o3///kCkK4Xp06fTo0cPunbtyscff5zUZy8rK+OCCy6ILqfi8+3atYvhw4fTrVs3unTpwgsvvMBDDz3Eli1bGDx4MIMHD47W8eWXXwLw9NNPk52dTbdu3aL/65kxYwYPPvgg8+fPp6CggPHjx5OTk8OePXsYNGhQtDuIhQsX0qdPH3r06MEVV1xBWVkZAFOnTqVTp05kZ2czefLkpD6jyLeCu6dt6Nmzpx9q3bp1VWdE/vilbkhAs2bNfOXKlT569Gjfs2ePd+vWzZcsWeLDhw93d/dp06b5M8884+7uX3/9tbdv397Lysp8165dvmfPHnd337Bhg1d8viVLlvgJJ5zgmzdv9vLycj/vvPP83XffrbLO0tJSb9my5WG1zJo1yydNmuRbt2713Nxcd3cfPXq05+bmelFRkc+dO9enTp3q7u7f/e53/aGHHnJ394cfftivvfbaw5ZXUUu3bt2iw6effhr93O7u+/fv99LSUnd3Ly4u9u9973t+8OBBnz9/vk+cODG6rB07dkTXW1xcHJ1fMb1mzRrv0KFD9LWSkhJ3d58+fbo/8MAD7u4+cOBA//vf/x59b8V0cXGxDxgwwMvKytzd/Te/+Y3/4he/8JKSEu/QoYMfPHgwuv0PddjvkEgdqUGspHi9FHgCGatDLoHs7GwKCwvJz89n2LBhVV5buHAhr732Gg8++CAQudTyX//6F23atOGmm25ixYoVNG7cmA0bNkTf07t3b7KysgDIycmhsLAwumd9JO6OmXHaaadRVlbGzp072bx5M+PGjeOdd97h3XffZdSoUdH2FeM9e/bkz3/+c8xlVnfIxd25/fbbeeedd2jUqBFffPEF27Zto2vXrkyePJnbbruNESNGMGDAgCPWvnjxYsaMGUOrVq0AOOmkk6r9vBX++te/sm7dOvr16wfAvn376NOnDyeccAJNmjRh4sSJDB8+vM7Pa4gczRTolYwcOZLJkyezdOlSSkpKovPdnZdffpmzzz67SvsZM2Zw6qmnsnLlSg4ePFjluPexxx4bHW/cuDEHDhyo8t4TTjiBZs2asWnTpujhFICPPvqIgQMHAtCnTx/++Mc/cvbZZzNgwACeeuopPvjgA377298etp5Y60jUs88+S3FxMcuXLycjI4O2bduyd+9eOnTowPLly1mwYAHTpk1jyJAh3H333XGXU/HHKBnuzoUXXkh+fv5hr3344Ye8/fbbPP/88/z+979n8eLFSa1DJOx0DL2SCRMmcPfdd9O1a9cq84cOHcrs2bOjx8H/8Y9/AFBaWkrr1q1p1KgRzzzzDOXl5TVa35QpU7j55pvZs2cPAG+99RbLli1j3LhxAJx//vk8+OCDnH/++XTv3p0lS5Zw7LHHcuKJJ9b2o1ZRWlrKKaecQkZGBkuWLOHzzyM9dW7ZsoXjjjuOq666ismTJ/PRRx8B0Lx5c3bu3HnYci644AJefPHF6B/Dr7766rA28d573nnn8d577/Hpp58CsHv3bjZs2EBZWRmlpaUMGzaMWbNmsWLFipR9bpGwaXh76PV4+jgrK4tbbrnlsPl33XUXt956K9nZ2bg7bdu25fXXX+fGG29k9OjRvPTSSwwePJhmzZrVaH0//elP+frrr+natSuNGzfmtNNO49VXX6Vp06ZA5FDJ5s2bOf/882ncuDFnnHEGHTt2rPHnevfdd8nJyYlO33nnnYwZMyY6PX78eC655BJyc3PJycmJrmP16tVMmTKFRo0akZGRwaOPPgpAXl4eF198Ma1bt2bJkiXR5XTu3Jk77riDgQMH0rhxY7p3787cuXOr1PKjH/2IG264gaZNm/LBBx9E52dmZjJ37lzGjh3LN99EroC95557aN68OZdeeil79+7F3Zk5c2aNP7/It4V5GgM0NzfXD33Axfr16znnnHPSVoOEj36HJF3q67JFM1vu7rnVtdMhFxGRkFCgi4iEhAJdRCQkFOgiIiGhQBcRCQkFuohISDS4QK+H3nMBeOWVVzCzpDu4SodDu55NxIwZMzj99NOrdJ+7Y8cOli5dmtRt9Dt27OCRRx6JTm/ZsqXKNe2xVO4i99e//nWV1/r27VvjGkQktgYX6PUlPz+f/v378/zzz6dsmcneih9PMoEOMGnSJFasWBEdWrRokXQNhwZ6mzZtmD9//hHf88QTT9CpUyfg8ECv6KFSRGpPgU6k+9j33nuPJ598skqgHzx4kBtvvJHOnTszYsQIhg0bFg2vBQsW0LFjR/r378/NN98c3dudMWMGeXl5DBkyhKuvvpry8nKmTJlCr169yM7O5vHHH6922b/85S/p1asXXbp0IS8vD3eP2fXs8uXLGThwID179mTo0KFs3bo1qc//4Ycf0rdvX7p3707fvn355JNPAFi7di29e/cmJyeH7OxsNm7cyNSpU/nss8/IyclhypQpFBYW0qVLFwDKy8uZPHkyXbt2JTs7m9mzZwNEu8idOnUqe/bsIScnh/HjxwNVH7LxwAMPRLfT9OnTgdhd+IpIHNV1x0jkwdAfAiuBtcAvgvlzgX8CK4Ihp7plJdJ9bj30nuvPPPOMT5gwwd3d+/Tp48uXL3d395deeskvvvhiLy8v961bt3qLFi38pZde8j179nhWVpZv2rTJ3d2vvPLKaFe706dP9x49evju3bvd3f3xxx/3X/3qV+7uvnfvXu/Zs6dv2rQp7rLd/9PtrLv7VVdd5a+99pq7V+16dt++fd6nTx/fvn27u7s///zz/uMf//iwzzZ9+nRv06ZNtOvcQYMGubtX6R64tLTU9+/f7+7uixYt8lGjRrm7+0033eR/+tOf3N39m2++8d27d/s///lP79y5c3T5lacfeeQRHzVqVHRZFZ+jct0VXfZWqJh+8803/brrrvODBw96eXm5Dx8+3P/yl7/E7cK3MnWfK+kShu5zvwG+7+5lZpYBLDOz/w1em+LuR/7/9lEgPz+fW2+9FYArr7yS/Px8evTowbJly7jiiito1KgRp512WvSBDh9//DFnnXUW7dq1A2Ds2LHMmTMnuryRI0dG+2NZuHAhq1atiu59l5aWsnHjxrjLBliyZAn3338/u3fv5quvvqJz585ccsklVWr+5JNPWLNmDRdeeCEQ2Ttu3bp1zM83adKkIz4YorS0lGuuuYaNGzdiZuzfvx+I9PZ47733UlRUxKhRo2jfvv0Rt+Nbb73FDTfcwDHHRH6tatJ97sKFC1m4cCHdu3cHIv9r2rhxIwMGDKhRF74i32bVBnrw16EsmMwIhgb6AKaaKykpYfHixaxZswYzo7y8HDPj/vvvj/mUISDu/AqVO+lyd2bPns3QoUOrtKl4OtKh9u7dy4033khBQQFnnHEGM2bMYO/evTFr6Ny5c5UOrpJ11113MXjwYF555RUKCwsZNGgQAOPGjePcc8/ljTfeYOjQoTzxxBNVuvqNVVNtus+dNm0a119//WGv1aQLX5Fvs4SOoZtZYzNbAWwHFrn734KX7jWzVWY208yOPcIiGqz58+dz9dVX8/nnn1NYWMjmzZtp164dy5Yto3///rz88sscPHiQbdu2sXTpUgA6duzIpk2bKCwsBDjicd2hQ4fy6KOPRvd6N2zYwK5du+IuuyK8W7VqRVlZWZUTjpW7nj377LMpLi6OBvr+/fuTfsB0aWlp9BmqlXtHrOir/eabb2bkyJGsWrUqbve3EHmE3mOPPRY9GRyr+9yMjIzotqhs6NChPPXUU9HHzn3xxRds3749bhe+InK4hALd3cvdPQfIAnqbWRdgGtAR6AWcBNwW671mlmdmBWZWUFxcnMC6UjtUJz8/n8svv7zKvNGjR/Pcc88xevRosrKy6NKlC9dffz3nnnsuJ554Ik2bNuWRRx7hoosuon///px66qlx+yifOHEinTp1okePHtHlHDhwIO6yW7RowXXXXUfXrl257LLL6NWrV3RZFV3P5uTkUF5ezvz587ntttvo1q0bOTk5ca8YmTlzZpXLFiv+EFX4+c9/zrRp0+jXr1+VPt1feOEFunTpQk5ODh9//DFXX301J598Mv369aNLly5MmTLlsM965plnRp8p+txzzx1WS15eHtnZ2dGTohWGDBnCuHHj6NOnD127dmXMmDHs3LmT1atXR0/M3nvvvdx5550xP6OIJNF9rplNB3a5+4OV5g0CJrv7ES9sPhq7zy0rK+P444+npKSE3r17895770UfD3f88cfj7vzkJz+hffv2TJo0KSXLlppp6L9DEh4Nvfvcao+hm1kmsN/dd5hZU+AHwH1m1trdt1rkoOllwJpaV90AjRgxgh07drBv3z7uuuuuaOD+4Q9/YN68eezbt4/u3bvHPPab7LJFRJJR7R66mWUD84DGRA7RvOjuvzSzxUAmYEQuW7zB3cviL+no3EOXhk+/Q5IuR/0euruvArrHmP/9JGuLtY6kr46Qb7eaHjIUCbN6v1O0SZMmlJSU6IspNebulJSU0KRJk/ouRaRBqPeHRGdlZVFUVEQiV8CIHKpJkyZkZWXVdxkiDUK9B3pGRkb0jksREUlevR9yERGR1FCgi4iEhAJdRCQkFOgiIiGhQBcRCQkFuohISCjQRURCQoEuIhISCnQRkZBQoIuIhIQCXUQkJBToIiIhoUAXEQkJBbqISEgo0EVEQqLaQDezJmb2oZmtNLO1ZvaLYH47M/ubmW00sxfM7Dt1X66IiMSTyB76N8D33b0bkANcZGbnAfcBM929PfA1cG3dlSkiItWpNtA9oiyYzAgGB74PzA/mzwMuq5MKRUQkIQkdQzezxma2AtgOLAI+A3a4+4GgSRFwet2UKCIiiUgo0N293N1zgCygN3BOrGax3mtmeWZWYGYFehC0iEjdqdFVLu6+A1gKnAe0MLOKh0xnAVvivGeOu+e6e25mZmZtahURkSNI5CqXTDNrEYw3BX4ArAeWAGOCZtcAr9ZVkSIiUr1jqm9Ca2CemTUm8gfgRXd/3czWAc+b2T3AP4An67BOERGpRrWB7u6rgO4x5m8icjxdREQaAN0pKiISEgp0EZGQUKCLiISEAl1EJCQU6CIiIaFAFxEJCQW6iEhIKNBFREJCgS4iEhIKdBGRkFCgi4iEhAJdRCQkFOgiIiGhQBcRCQkFuohISCjQRURCQoEuIhISiTxT9AwzW2Jm681srZndEsyfYWZfmNmKYBhW9+WKiEg8iTxT9ADw3+7+kZk1B5ab2aLgtZnu/mDdlSciIolK5JmiW4GtwfhOM1sPnF7XhYmISM3U6Bi6mbUl8sDovwWzbjKzVWb2lJm1THFtIiJSAwkHupkdD7wM3Oru/wYeBb4H5BDZg/9tnPflmVmBmRUUFxenoGQREYkloUA3swwiYf6su/8ZwN23uXu5ux8E/gD0jvVed5/j7rnunpuZmZmqukVE5BCJXOViwJPAenf/n0rzW1dqdjmwJvXliYhIohK5yqUf8ENgtZmtCObdDow1sxzAgULg+jqpUEREEpLIVS7LAIvx0oLUlyMiIsnSnaIiIiGhQBcRCQkFuohISCjQRURCQoEuIhISCnQRkZBQoIuIhIQCXUQkJBToIiIhoUAXEQkJBbqISEgo0EVEDmEWGY42CnQRkZBQoIuIhIQCXUQkJBToIiIhoUAXEQkJBbqISEgk8pDoM8xsiZmtN7O1ZnZLMP8kM1tkZhuDny3rvlwREYknkT30A8B/u/s5wHnAT8ysEzAVeNvd2wNvB9MiIlJPqg10d9/q7h8F4zuB9cDpwKXAvKDZPOCyuipSRESqV6Nj6GbWFugO/A041d23QiT0gVNSXZyIiCQu4UA3s+OBl4Fb3f3fNXhfnpkVmFlBcXFxMjWKiEgCEgp0M8sgEubPuvufg9nbzKx18HprYHus97r7HHfPdffczMzMVNQsIiIxJHKViwFPAuvd/X8qvfQacE0wfg3waurLExGRRB2TQJt+wA+B1Wa2Iph3O/Ab4EUzuxb4F3BF3ZQoIiKJqDbQ3X0ZEK8jyQtSW46IiCRLd4qKiISEAl1EJCQU6CIiIaFAFxEJCQW6iEhIKNBFREJCgS4iEhIKdBGRkFCgi4iEhAJdRCQkFOgiIiGhQBcRCQkFuohISCjQRURCQoEuIhISCnQRkZBQoIuIhEQizxR9ysy2m9maSvNmmNkXZrYiGIbVbZkiIlKdRPbQ5wIXxZg/091zgmFBassSEZGaqjbQ3f0d4Ks01CIiIrVQm2PoN5nZquCQTMuUVSQiIklJNtAfBb4H5ABbgd/Ga2hmeWZWYGYFxcXFSa5ORKThMosM9S2pQHf3be5e7u4HgT8AvY/Qdo6757p7bmZmZrJ1iohINZIKdDNrXWnycmBNvLYiIpIex1TXwMzygUFAKzMrAqYDg8wsB3CgELi+DmsUEZEEVBvo7j42xuwn66AWERGpBd0pKiISEgp0EZGQUKCLiISEAl1EJCQU6CIiIaFAFxEJCQW6iEhIKNBFREJCgS4iEhIKdBGRkFCgi4iEhAJdRCQkFOgiIiGhQBcRCQkFuohISCjQRURCotoHXEhIVDzB1r1+6xD5Fqn84Oh0fPW0hy4iEhLVBrqZPWVm281sTaV5J5nZIjPbGPxsWbdliohIdRLZQ58LXHTIvKnA2+7eHng7mBYRkXpUbaC7+zvAV4fMvhSYF4zPAy5LcV0iIlJDyR5DP9XdtwIEP09JXUkiIpKMOj8pamZ5ZlZgZgXFxcV1vbqGzew/Q12vQ0S+dZIN9G1m1hog+Lk9XkN3n+Puue6em5mZmeTqRESkOskG+mvANcH4NcCrqSlHRESSlchli/nAB8DZZlZkZtcCvwEuNLONwIXBtIiI1KNq7xR197FxXrogxbWIiEgt6Nb/o1G67ycWkaOCbv0XEQkJBbqISEgo0EVEQkKBLiISEjop2pDp5KeI1ID20EVEQkKBLiISEgp0EZGQUKCLiISEAl1EJCQU6JXF60tcfYyLyFFAgS4iEhIKdBGRkFCgi4iEhAJdRCQkdOt/bVScKK3tbfm6xV9EUkB76CIiIVGrPXQzKwR2AuXAAXfPTUVRIiJSc6k45DLY3b9MwXJERKQWdMhFRCQkahvoDiw0s+VmlhergZnlmVmBmRUUFxfXcnVHCd1ZKnJUCNtXtbaB3s/dewAXAz8xs/MPbeDuc9w9191zMzMza7k6ERGJp1aB7u5bgp/bgVeA3qkoSkREai7pQDezZmbWvGIcGAKsSVVhIiJSM7W5yuVU4BWLHIA6BnjO3f9vSqoSEZEaSzrQ3X0T0C2FtYiISC3oskURkZBQoIuIhIQCXUQkJBToIiIhoUAXEQmJ8PaHrj7GExOvT/dU9fUuImmjPXQRkZBQoIuIhIQCXUQkJBToIiIhcfSfFE3m5KdO+IlICGkPXUQkJBToIiIhoUAXEQkJBbqISEgcPSdFdednw6ATytJAHSkivi2/ttpDFxEJiVoFupldZGafmNmnZjY1VUWJiEjN1eYh0Y2Bh4GLgU7AWDPrlKrCRESkZmqzh94b+NTdN7n7PuB54NLUlCUiIjVVm0A/HdhcaboomCciIvWgNle5WIx5h51DNrM8IC+YLDOzT4LxVsCXya051qpTOj92balafu2WVbW2ut8Wic5P/t+z7jXU2hpqXXCU11ZPX4tWZvHrOlIcJOC7iTSqTaAXAWdUms4CthzayN3nAHMOnW9mBe6eW4v11xnVVnMNtS5ouLU11LpAtSWjIdRVm0Mufwfam1k7M/sOcCXwWmrKEhGRmkp6D93dD5jZTcCbQGPgKXdfm7LKRESkRmp1p6i7LwAWJPn2ww7DNCCqreYaal3QcGtrqHWBaktGvddlHvZ7YUVEviV067+ISEikPdDNbLKZuUUuDYz1ermZrQiGtJ5kTaC2a8xsYzBck4Z6fmVmq4JtsdDM2sRpl/ZtVoPa0rrNgnU+YGYfB/W9YmYt4rQrNLPVwWcoaEB1pb1LDTO7wszWmtlBM4t7pUa6t1kNa0vrdjOzk8xsUfC7vcjMWsZpl77vp7unbSBymeObwOdAqzhtytJZU6K1AScBm4KfLYPxlnVc0wmVxm8GHmso2yyR2upjmwXrHQIcE4zfB9wXp11hvN/D+qqLyAUGnwFnAd8BVgKd0lDbOcDZwFIg9wjt0rrNEq2tPrYbcD8wNRifeoTfs7R9P9O9hz4T+DkxbkBqAKqrbSiwyN2/cvevgUXARXVZkLv/u9JksyPUlnYJ1pb2bRbUttDdDwSTfyVyj0S9S7CueulSw93Xu/sn1bdMvwRrq4/tdikwLxifB1xWx+urVtoC3cxGAl+4+8pqmjYxswIz+6uZpWUDJVhbvXR1YGb3mtlmYDxwd5xmad9mCdbWELqHmAD8b5zXHFhoZsuDO5rTKV5dDWGbHUl9brMjqY/tdqq7bwUIfp4Sp13avp8pfcCFmb0FnBbjpTuA24n8l7M6Z7r7FjM7C1hsZqvd/bMGUFtCXR2ksi53f9Xd7wDuMLNpwE3A9Bht077NEqytTrZZIrUFbe4ADgDPxllMv2C7nQIsMrOP3f2deq6rXrdZAlK+zVJUW9q/nzVYTJ18P2NJaaC7+w9izTezrkA7YKVFOjTIAj4ys97u/v8OWcaW4OcmM1sKdCdybKy+aysCBlWaziJyTK9O6orhOeANYgR6urdZDWqrk20G1dcWnIAdAVzgwYHMGMuo2G7bzewVIv9tr1U4paCuhLrUqIvaElxGyrdZimqrk+12pLrMbJuZtXb3rWbWGtgeZxl18v2MKV0H6w85SVBI7BOPLYFjg/FWwEbScEIowdpOAv4Z1NgyGD+pjmtpX2n8p8D8hrLNEqwt7dssWO9FwDog8whtmgHNK42/D1zUAOo6hsjJ43b85+Re57reZpXWv5T4Jx7Tvs1qUFvatxvwAFVPit4fo01av59p+YeI8SGjoQnkAk8E432B1cE/xmrg2oZSWzA9Afg0GH6chlpeBtYAq4D/A5zeULZZIrXVxzYL1vkpkeOpK4LhsWB+G2BBMH5WsM1WAmuJ/Ne+3usKpocBG4jsxdV5XcE6Lyeyl/sNsA14syFss0Rrq4/tBpwMvB2E9NsEOyv1+f3UnaIiIiGhO0VFREJCgS4iEhIKdBGRkFCgi4iEhAJdRCQkFOgiIiGhQBcRCQkFuohISPx/tXEmUvnqJloAAAAASUVORK5CYII=\n",
"text/plain": [
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plt.hist(\n",
" [means.flatten(), aggregates.flatten()], \n",
" color=['red', 'blue'], \n",
" bins=50\n",
");\n",
"plt.legend(['Mean Own Elasticities', 'Aggregate Elasticities']);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Marginal Costs and Markups\n",
"\n",
"To compute marginal costs, $c$, the `product_data` passed to [`Problem`](https://pyblp.readthedocs.io/en/stable/_api/pyblp.Problem.html#pyblp.Problem) must have had a `firm_ids` field. Since we included firm IDs when configuring the problem, we can use [`ProblemResults.compute_costs`](https://pyblp.readthedocs.io/en/stable/_api/pyblp.ProblemResults.compute_costs.html#pyblp.ProblemResults.compute_costs)."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAF4NJREFUeJzt3XuQnXWd5/H3l5AQBJRLWomJTOIWTOSaDD0kExgWiDICUSwkJgxKFCx2l0HdZZWLDGil1EKkkLG0xBiE4GISQSyC6LBZLioWl+mYAEHEEIjQkoWYEZarIea7f5ynM4fQne4+5/Qlv7xfVafOOb/n9s3TTz7969/znOdEZiJJKtdOQ12AJGlgGfSSVDiDXpIKZ9BLUuEMekkqnEEvSYUz6CWpcAa9JBXOoJekwu081AUAjBkzJidMmDDUZUjSdmX58uV/zMy23uYbFkE/YcIEOjo6hroMSdquRMTv+zKfQzeSVDiDXpIKZ9BLUuGGxRi9pMHz+uuv09nZyWuvvTbUpaiPRo8ezfjx4xk5cmRDyxv00g6ms7OTPfbYgwkTJhARQ12OepGZbNiwgc7OTiZOnNjQOhy6kXYwr732Gvvss48hv52ICPbZZ5+m/gIz6KUdkCG/fWn252XQS1LhHKOXdnATLrytpetbe9lJvc4TEXz0ox/l+9//PgCbNm1i7NixTJ06lZ/85CctqeOTn/wk5513HgceeGC/l73uuuvo6Ojgm9/85pum/exnP+OSSy7h5ZdfJjOZOXMmV1xxRb/Wv3LlSp555hlOPPHEftfWCINeO5RthVpfAkqtsdtuu7Fq1SpeffVVdt11V5YtW8a4ceP6tY5Nmzax8849R9iCBQuaLfNNVq1axbnnnsttt93GpEmT2LRpE/Pnz+/3elauXElHR8egBb1DN5KGxAknnMBtt9V+8S5atIjTTjtty7QHHniA6dOnM2XKFKZPn85jjz0G1Hras2bN4gMf+ADHH388mzdv5pxzzuGggw5i5syZnHjiidx0000AHHPMMVturbL77rtz8cUXc9hhhzFt2jSeffZZAG699VamTp3KlClTeO9737ulvSeXX345F198MZMmTQJg55135pxzzgHg97//PTNmzODQQw9lxowZPPXUUwDceOONHHzwwRx22GEcffTRbNy4kUsvvZQlS5YwefJklixZws9//nMmT57M5MmTmTJlCi+++GKrdjNg0EsaInPmzGHx4sW89tprPPTQQ0ydOnXLtEmTJvGLX/yCFStWMG/ePD7/+c9vmXbvvfeycOFC7rzzTm6++WbWrl3Lww8/zIIFC7j33nu73dbLL7/MtGnTePDBBzn66KP57ne/C8BRRx3Ffffdx4oVK5gzZw6XX375NmtetWoVhx9+eLfTzj33XM444wweeughTj/9dD796U8DMG/ePG6//XYefPBBli5dyqhRo5g3bx6zZ89m5cqVzJ49myuuuIJvfetbrFy5kl/+8pfsuuuu/dqXvXHoRtu1noZiHIYZ/g499FDWrl3LokWL3jSE8cILLzB37lxWr15NRPD6669vmfa+972PvffeG4B77rmHWbNmsdNOO7Hvvvty7LHHdrutUaNGMXPmTAAOP/xwli1bBtQ+UzB79mzWrVvHxo0bG75OHWq/gG6++WYAPvaxj3H++ecDcOSRR/Lxj3+cj3zkI5xyyindLnvkkUdy3nnncfrpp3PKKacwfvz4huvojj16SUPmgx/8IJ/97GffMGwDcMkll3DssceyatUqbr311jdcQ77bbrtteZ2ZfdrOyJEjt1yiOGLECDZt2gTApz71Kc4991wefvhhvvOd7/R6rfpBBx3E8uXL+7TNru1dffXVfOlLX+Lpp59m8uTJbNiw4U3zXnjhhSxYsIBXX32VadOm8dvf/rZP2+grg17SkDnzzDO59NJLOeSQQ97Q/sILL2w5OXvdddf1uPxRRx3Fj370IzZv3syzzz7L3Xff3a/t129n4cKFvc7/uc99jq985Sv87ne/A2Dz5s1ceeWVAEyfPp3FixcDcMMNN3DUUUcBsGbNGqZOncq8efMYM2YMTz/9NHvssccbxuHXrFnDIYccwgUXXEB7e3vLg77XoZuI+B4wE3guMw+u2r4GfADYCKwBPpGZz1fTLgLOAv4CfDozb29pxZJaaiiHucaPH89nPvOZN7Wff/75zJ07lyuvvJLjjjuux+U//OEPc8cdd3DwwQdzwAEHMHXqVN72trf1eftf/OIXmTVrFuPGjWPatGk8+eST25z/0EMP5aqrruK0007jlVdeISI46aTa/vvGN77BmWeeyde+9jXa2tq49tprgdovh9WrV5OZzJgxg8MOO4z99tuPyy67jMmTJ3PRRRdxzz33cNdddzFixAgOPPBATjjhhD7/G/oievvTJyKOBl4Crq8L+uOBOzNzU0R8FSAzL4iIA4FFwBHAO4H/AxyQmX/Z1jba29vTLx5RI/o7Ru/llfDoo4/ynve8Z6jLaJmXXnqJ3XffnQ0bNnDEEUfwq1/9in333Xeoy2q57n5uEbE8M9t7W7bXHn1m/iIiJmzV9r/r3t4HnFq9PhlYnJl/Bp6MiMephX73p8IlqUkzZ87k+eefZ+PGjVxyySVFhnyzWnHVzZnAkur1OGrB36WzapOkAdHfcfkdUVNBHxEXA5uAG7qaupmt27GhiDgbOBtgv/32a6YMaUCVeAlnZnpjs+1IX68u6knDV91ExFxqJ2lPz/+oohN4V91s44Fnuls+M+dnZntmtre19fol5pJaZPTo0WzYsKHp8NDg6Lof/ejRoxteR0M9+oh4P3AB8J8z85W6SUuBH0TEldROxu4PPNBwdZJabvz48XR2drJ+/fqhLkV91PUNU43qy+WVi4BjgDER0Ql8AbgI2AVYVv35d19m/tfMfCQifgj8htqQzj/1dsWNpME1cuTIpj4Bqu1PX666Oa2b5mu2Mf+XgS83U5QkqXW8141UafV92aXhwlsgSFLhDHpJKpxDNyqSwzDSfzDopQaV+EEqlcmhG0kqnEEvSYUz6CWpcAa9JBXOoJekwhn0klQ4g16SCmfQS1LhDHpJKpxBL0mFM+glqXDe60YaJN4bR0PFHr0kFc6gl6TCGfSSVDiDXpIKZ9BLUuEMekkqnEEvSYXrNegj4nsR8VxErKpr2zsilkXE6up5r6o9IuIbEfF4RDwUEX8zkMVLknrXlx79dcD7t2q7ELgjM/cH7qjeA5wA7F89zga+3ZoyJUmN6vWTsZn5i4iYsFXzycAx1euFwN3ABVX79ZmZwH0RsWdEjM3Mda0qWBruevoErDRUGh2jf0dXeFfPb6/axwFP183XWbW9SUScHREdEdGxfv36BsuQJPWm1Sdjo5u27G7GzJyfme2Z2d7W1tbiMiRJXRoN+mcjYixA9fxc1d4JvKtuvvHAM42XJ0lqVqNBvxSYW72eC9xS135GdfXNNOAFx+claWj1ejI2IhZRO/E6JiI6gS8AlwE/jIizgKeAWdXsPwVOBB4HXgE+MQA1qwD9vWWvJzilxvXlqpvTepg0o5t5E/inZouSJLWOn4yVpML5DVMaVhyikVrPHr0kFc6gl6TCGfSSVDiDXpIKZ9BLUuEMekkqnEEvSYUz6CWpcAa9JBXOoJekwhn0klQ4g16SCmfQS1LhDHpJKpxBL0mFM+glqXAGvSQVzqCXpMIZ9JJUOINekgpn0EtS4Qx6SSpcU0EfEf8jIh6JiFURsSgiRkfExIi4PyJWR8SSiBjVqmIlSf3XcNBHxDjg00B7Zh4MjADmAF8Fvp6Z+wN/As5qRaGSpMbs3ILld42I14G3AOuA44B/rKYvBL4IfLvJ7UjFmnDhbd22r73spEGuRKVquEefmX8ArgCeohbwLwDLgeczc1M1WycwrrvlI+LsiOiIiI7169c3WoYkqRfNDN3sBZwMTATeCewGnNDNrNnd8pk5PzPbM7O9ra2t0TIkSb1o5mTse4EnM3N9Zr4O3AxMB/aMiK4hofHAM03WKElqQjNB/xQwLSLeEhEBzAB+A9wFnFrNMxe4pbkSJUnNaPhkbGbeHxE3Ab8GNgErgPnAbcDiiPhS1XZNKwrV9qmnE42SBk9TV91k5heAL2zV/ARwRDPrldQzr9JRf/nJWEkqnEEvSYUz6CWpcAa9JBXOoJekwhn0klS4Zm9qJgFeLz8Q3KdqFXv0klQ4g16SCmfQS1LhDHpJKpxBL0mF86obqRDe7Ew9sUcvSYUz6CWpcAa9JBXOoJekwhn0klQ4g16SCmfQS1LhDHpJKpxBL0mFM+glqXAGvSQVrqmgj4g9I+KmiPhtRDwaEX8XEXtHxLKIWF0979WqYiVJ/ddsj/5fgH/NzEnAYcCjwIXAHZm5P3BH9V6SNEQaDvqIeCtwNHANQGZuzMzngZOBhdVsC4EPNVukJKlxzfTo3w2sB66NiBURsSAidgPekZnrAKrnt3e3cEScHREdEdGxfv36JsqQJG1LM0G/M/A3wLczcwrwMv0YpsnM+ZnZnpntbW1tTZQhSdqWZr54pBPozMz7q/c3UQv6ZyNibGaui4ixwHPNFimpcX4hiRru0Wfm/wWejoi/rppmAL8BlgJzq7a5wC1NVShJakqzXyX4KeCGiBgFPAF8gtovjx9GxFnAU8CsJrchSWpCU0GfmSuB9m4mzWhmvRq+ehoGkDR8+clYSSqcQS9JhTPoJalwBr0kFc6gl6TCGfSSVDiDXpIKZ9BLUuEMekkqnEEvSYUz6CWpcAa9JBXOoJekwhn0klQ4g16SCmfQS1LhDHpJKlyzXyWoQvlNUlI57NFLUuEMekkqnEM30g6qp+G5tZedNMiVaKDZo5ekwhn0klS4poM+IkZExIqI+En1fmJE3B8RqyNiSUSMar5MSVKjWtGj/wzwaN37rwJfz8z9gT8BZ7VgG5KkBjUV9BExHjgJWFC9D+A44KZqloXAh5rZhiSpOc326K8Czgc2V+/3AZ7PzE3V+05gXJPbkCQ1oeGgj4iZwHOZuby+uZtZs4flz46IjojoWL9+faNlSJJ60UyP/kjggxGxFlhMbcjmKmDPiOi6Pn888Ex3C2fm/Mxsz8z2tra2JsqQJG1Lwx+YysyLgIsAIuIY4LOZeXpE3AicSi385wK3tKBOSYNkW/c58sNU26eBuI7+AuC8iHic2pj9NQOwDUlSH7XkFgiZeTdwd/X6CeCIVqxXktQ8PxkrSYUz6CWpcAa9JBXOoJekwhn0klQ4g16SCmfQS1LhDHpJKpxBL0mFM+glqXAtuQWCpB1DTzc882Znw5tBvwPY1t0IpVbwF8Dw5tCNJBXOoJekwhn0klQ4g16SCmfQS1LhDHpJKpxBL0mFM+glqXAGvSQVzqCXpMIZ9JJUOO91I2nQeW+cwdVwjz4i3hURd0XEoxHxSER8pmrfOyKWRcTq6nmv1pUrSeqvZnr0m4D/mZm/jog9gOURsQz4OHBHZl4WERcCFwIXNF/qjsdej7Z33jl1eGg46DNzHbCuev1iRDwKjANOBo6pZlsI3I1BPyj8TyWpOy05GRsRE4ApwP3AO6pfAl2/DN7eim1IkhrTdNBHxO7Aj4D/npn/rx/LnR0RHRHRsX79+mbLkCT1oKmgj4iR1EL+hsy8uWp+NiLGVtPHAs91t2xmzs/M9sxsb2tra6YMSdI2NDxGHxEBXAM8mplX1k1aCswFLqueb2mqQkk7DC9AGBjNXHVzJPAx4OGIWFm1fZ5awP8wIs4CngJmNVeiJKkZzVx1cw8QPUye0eh6JUmt5S0QJKlwBr0kFc6gl6TCGfSSVDiDXpIKZ9BLUuG8H/0w4M3IpMb4Aau+sUcvSYWzRz+IWtVz9y8A7Wg85ptjj16SCmfQS1LhDHpJKpxBL0mFM+glqXAGvSQVzssrJe3wSv/glUHfhNIPDml75f/NN3LoRpIKZ49+APgpPml42lH/b9qjl6TC2aPvgx21FyCpDAa9JDVgezrh69CNJBXOoJekwg1Y0EfE+yPisYh4PCIuHKjtSJK2LTKz9SuNGAH8Dngf0An8G3BaZv6mu/nb29uzo6Oj5XX0dwzNk66SBspAjN1HxPLMbO9tvoE6GXsE8HhmPlEVsxg4Geg26JvRSDgb6JIG21CevB2ooZtxwNN17zurNknSIBuoHn100/aGMaKIOBs4u3r7UkQ8NkC19NcY4I9DXcQ2WF/jhnNtYH3NGM61wTbqi682td6/6stMAxX0ncC76t6PB56pnyEz5wPzB2j7DYuIjr6MeQ0V62vccK4NrK8Zw7k2GPr6Bmro5t+A/SNiYkSMAuYASwdoW5KkbRiQHn1mboqIc4HbgRHA9zLzkYHYliRp2wbsFgiZ+VPgpwO1/gE07IaTtmJ9jRvOtYH1NWM41wZDXN+AXEcvSRo+vAWCJBWu+KDv7VYMEbFLRCyppt8fERPqpl1UtT8WEf/Q13UOdG0R8b6IWB4RD1fPx9Utc3e1zpXV4+1DUN+EiHi1roar65Y5vKr78Yj4RkR0dynuQNd3el1tKyNic0RMrqYN5v47OiJ+HRGbIuLUrabNjYjV1WNuXXtL9l+jtUXE5Ii4NyIeiYiHImJ23bTrIuLJun03uZHamqmvmvaXuhqW1rVPrI6D1dVxMWowa4uIY7c67l6LiA9V01q277qVmcU+qJ0IXgO8GxgFPAgcuNU85wBXV6/nAEuq1wdW8+8CTKzWM6Iv6xyE2qYA76xeHwz8oW6Zu4H2Id53E4BVPaz3AeDvqH3W4mfACYNd31bzHAI8MUT7bwJwKHA9cGpd+97AE9XzXtXrvVq1/5qs7QBg/+r1O4F1wJ7V++vq5x2KfVdNe6mH9f4QmFO9vhr4b4Nd21Y/438H3tLKfdfTo/Qe/ZZbMWTmRqDrVgz1TgYWVq9vAmZUvaSTgcWZ+efMfBJ4vFpfX9Y5oLVl5orM7PpcwiPA6IjYpYEaBqS+nlYYEWOBt2bmvVk7uq8HPjTE9Z0GLGqwhqbqy8y1mfkQsHmrZf8BWJaZ/56ZfwKWAe9v4f5ruLbM/F1mrq5ePwM8B7Q1UMOA1NeT6ud+HLXjAGrHxaDuu62cCvwsM19poIZ+Kz3o+3Irhi3zZOYm4AVgn20s26rbOzRTW70PAysy8891bddWf/5d0sTQSLP1TYyIFRHx84j4+7r5O3tZ52DV12U2bw76wdp//V22VfuvJcdwRBxBrVe7pq75y9WQzteb6Hw0W9/oiOiIiPu6hkao/dyfr46DRtbZqtq6zOHNx10r9l23Sg/6Xm/FsI15+tveX83UVpsYcRDwVeC/1E0/PTMPAf6+enysgdqarW8dsF9mTgHOA34QEW/t4zoHo77axIipwCuZuapu+mDuv/4uO5jH3rZXUPvr4vvAJzKzq+d6ETAJ+FtqQxMXNFBbK+rbL2ufQv1H4KqI+E8tWGerauvad4dQ+5xRl1btu26VHvS93oqhfp6I2Bl4G7Wxs56W7cs6B7o2ImI88GPgjMzc0qPKzD9Uzy8CP6D2p2YjGq6vGu7aUNWxnFqP74Bq/vG9rHPA66ub/qZe1SDvv/4u26r919QxXP3Svg3458y8r6s9M9dlzZ+Baxmafdc1pETW7p57N7VzWn8E9qyOg36vs1W1VT4C/DgzX6+ruVX7rlulB31fbsWwFOi6quFU4M5q/HMpMCdqV25MBPandiKsVbd3aLi2iNiT2n+0izLzV10zR8TOETGmej0SmAmsojHN1NcWte8kICLeTW3fPZGZ64AXI2JaNSRyBnDLYNdX1bUTMIvaGCtV22Dvv57cDhwfEXtFxF7A8cDtLdx/DddWzf9j4PrMvHGraWOr56A2/j3o+67aZ7tUr8cARwK/qX7ud1E7DqB2XAzqvqvzpvNCLdx33Ruos7zD5QGcSO1LUNYAF1dt84APVq9HAzdSO9n6APDuumUvrpZ7jLqrG7pb52DWBvwz8DKwsu7xdmA3YDnwELWTtP8CjBiC+j5cbf9B4NfAB+rW2U7tIF4DfJPqQ3tD8LM9Brhvq/UN9v77W2o9xJeBDcAjdcueWdX9OLXhkZbuv0ZrAz4KvL7VsTe5mnYn8HBV3/8Cdh/sfQdMr2p4sHo+q26d766Og8er42KXIfi5TgD+AOy01Tpbtu+6e/jJWEkqXOlDN5K0wzPoJalwBr0kFc6gl6TCGfSSVDiDXpIKZ9BLUuEMekkq3P8HFYRJLRECGdsAAAAASUVORK5CYII=\n",
"text/plain": [
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"costs = results.compute_costs()\n",
"plt.hist(costs, bins=50);\n",
"plt.legend([\"Marginal Costs\"]);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Other methods that compute supplyside outputs often compute marginal costs themselves. For example, [`ProblemResults.compute_markups`](https://pyblp.readthedocs.io/en/stable/_api/pyblp.ProblemResults.compute_markups.html#pyblp.ProblemResults.compute_markups) will compute marginal costs when estimating markups, $\\mathscr{M}$, but computation can be sped up if we just use our precomputed values."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAFWhJREFUeJzt3X+Q3XV97/HnuxIbiSA/stDc/OhCDQLFNPGuBEEtSKEkVaKFeMkUGgRdYMSxtlxNqVaGTmeYNsDQseqNJZPE0QC9/ApKf3C1lMsFhA2haTCm5ccWt2RIDN6oE+I15H3/ON/Ek7jJnpwfe7KfPB8zO+ecz/fX+5PdvPa7n/M9n29kJpKkcv1StwuQJHWWQS9JhTPoJalwBr0kFc6gl6TCGfSSVDiDXpIKZ9BLUuEMekkq3GHdLgBg4sSJ2dvb2+0yJGlMWb169Q8ys2ek9Q6KoO/t7WVgYKDbZUjSmBIR/9HIeg7dSFLhDHpJKpxBL0mFOyjG6CUdOn72s58xNDTE9u3bu13KmDF+/HimTJnCuHHjmtreoJc0qoaGhjjiiCPo7e0lIrpdzkEvM9myZQtDQ0OccMIJTe3DoRtJo2r79u0ce+yxhnyDIoJjjz22pb+ADHpJo86QPzCt/nsZ9JJUOMfoJXVV76JvtnV/gzf9zojrRASXXnopX/3qVwHYsWMHkyZNYvbs2XzjG99o+FgPP/wwixcvPqBtusGg74B9/eA28gMoqfMmTJjAunXreO2113jTm97EQw89xOTJkw9oHzt27OhQde3n0I2kQ9KcOXP45jdrJ2UrV65kwYIFu5c9+eSTnHnmmcyaNYszzzyTDRs2ALBs2TLmz5/PBz7wAc4///w99vfUU08xa9YsXnjhBW644QYWL168e9lpp53G4OAgg4ODnHzyySxcuJAZM2Zw8cUXs23bNgAWLVrEqaeeyowZM7juuuva2leDXtIh6ZJLLuGOO+5g+/btrF27ltmzZ+9edvLJJ/PII4+wZs0abrzxRq6//vrdyx5//HGWL1/Ot7/97d1tjz32GFdffTX3338/J5544n6Pu2HDBvr7+1m7di1HHnkkX/ziF3n11Ve59957efbZZ1m7di2f/exn29pXg17SIWnGjBkMDg6ycuVK5s6du8eyrVu3Mn/+fE477TQ+9alP8eyzz+5edt5553HMMcfsfr1+/Xr6+/t54IEHmDZt2ojHnTp1KmeddRYAl156KY8++ihHHnkk48eP56Mf/Sj33HMPhx9+eJt6WWPQSzpkXXjhhVx33XV7DNsAfO5zn+Occ85h3bp1PPDAA3tcwz5hwoQ91p00aRLjx49nzZo1u9sOO+wwdu7cuft1/fZ7XyoZERx22GE8+eSTXHTRRdx3331ccMEFbenf7nraujdJGkOuuOIK3vKWt/D2t7+dhx9+eHf71q1bd785u2zZsv3u46ijjuL222/n/PPPZ8KECZx99tn09vbuvhLn6aef5sUXX9y9/ksvvcTjjz/Ou971LlauXMm73/1ufvKTn7Bt2zbmzp3LGWecwVvf+ta29nPEoI+IqcAK4FeAncCSzLwtIo4B7gR6gUHgw5n5w6j9uroNmAtsAy7PzKfbWrWkYnTzarQpU6bwyU9+8hfaP/3pT7Nw4UJuueUW3ve+9424n+OPP54HHniAOXPmsHTpUi666CJWrFjBzJkzeec738lJJ520e91TTjmF5cuXc9VVVzF9+nSuueYatm7dyrx589i+fTuZya233trWfkZm7n+FiEnApMx8OiKOAFYDHwQuB17NzJsiYhFwdGZ+JiLmAp+gFvSzgdsyc/Y+dg9AX19flnTjES+vlPZt/fr1nHLKKd0uoysGBwd5//vfz7p16w542+H+3SJidWb2jbTtiGP0mblx1xl5Zv4YWA9MBuYBy6vVllMLf6r2FVnzBHBU9ctCktQFBzRGHxG9wCzgO8DxmbkRar8MIuK4arXJwPfrNhuq2ja2Wmy3eIYuqR16e3ubOptvVcNX3UTEm4G7gT/IzB/tb9Vh2n5hfCgi+iNiICIGNm/e3GgZkgow0pCx9tTqv1dDZ/QRMY5ayH8tM++pml+JiEnV2fwkYFPVPgRMrdt8CvDy3vvMzCXAEqiN0TdZf1e1e44O6VAwfvx4tmzZ4lTFDdo1H/348eOb3kcjV90EcDuwPjNvqVu0ClgI3FQ93l/Xfm1E3EHtzditu4Z4JGnKlCkMDQ3hX/KN23WHqWY1ckZ/FnAZ8K8R8UzVdj21gL8rIq4EXgLmV8sepHbFzXPULq/8SNPVSSrOuHHjmr5TkpozYtBn5qMMP+4OcO4w6yfw8RbrkiS1iVMgSFLhDHpJKpxBL0mFM+glqXAGvSQVzqCXpMI5H/0ocs4cSd3gGb0kFc6gl6TCGfSSVDiDXpIKZ9BLUuEMekkqnEEvSYUz6CWpcAa9JBVuxKCPiKURsSki1tW13RkRz1Rfg7vuPBURvRHxWt2yL3eyeEnSyBqZAmEZ8AVgxa6GzPxvu55HxM3A1rr1n8/Mme0qUJLUmkZuJfhIRPQOt6y6cfiHgfe1tyxJUru0Okb/HuCVzPz3urYTImJNRPxzRLxnXxtGRH9EDETEgHeDl6TOaTXoFwAr615vBKZl5izgD4GvR8SRw22YmUsysy8z+3p6elosQ5K0L00HfUQcBvwucOeutsz8aWZuqZ6vBp4HTmq1SElS81qZj/63gO9l5tCuhojoAV7NzNcj4kRgOvBCizUWz3nqJXVSI5dXrgQeB94WEUMRcWW16BL2HLYBeC+wNiL+BfifwNWZ+Wo7C5YkHZhGrrpZsI/2y4dpuxu4u/WyJEnt4idjJalwBr0kFc6gl6TCGfSSVDiDXpIKZ9BLUuEMekkqnEEvSYUz6CWpcAa9JBXOoJekwrUye2Vx9jWLpCSNZZ7RS1LhDHpJKpxBL0mFa+TGI0sjYlNErKtruyEi/jMinqm+5tYt++OIeC4iNkTEb3eqcElSYxp5M3YZ8AVgxV7tt2bm4vqGiDiV2p2nfh34L8D/ioiTMvP1NtR6yPEWg5LaYcQz+sx8BGj0doDzgDuqm4S/CDwHnN5CfZKkFrUyRn9tRKythnaOrtomA9+vW2eoapMkdUmzQf8l4NeAmcBG4OaqPYZZN4fbQUT0R8RARAxs3ry5yTIkSSNpKugz85XMfD0zdwJf4efDM0PA1LpVpwAv72MfSzKzLzP7enp6milDktSApoI+IibVvfwQsOuKnFXAJRHxyxFxAjAdeLK1EiVJrRjxqpuIWAmcDUyMiCHg88DZETGT2rDMIHAVQGY+GxF3Ad8FdgAf94obSequEYM+MxcM03z7ftb/c+DPWylKktQ+fjJWkgpn0EtS4Qx6SSqcQS9JhTPoJalwBr0kFc6gl6TCGfSSVDiDXpIKZ9BLUuEMekkqnEEvSYUz6CWpcAa9JBXOoJekwhn0klS4EYM+IpZGxKaIWFfX9pcR8b2IWBsR90bEUVV7b0S8FhHPVF9f7mTxkqSRNXJGvwy4YK+2h4DTMnMG8G/AH9ctez4zZ1ZfV7enTElSsxq5leAjEdG7V9s/1r18Ari4vWV1Vu+ib3a7BEkaNe0Yo78C+Lu61ydExJqI+OeIeE8b9i9JasGIZ/T7ExF/AuwAvlY1bQSmZeaWiPivwH0R8euZ+aNhtu0H+gGmTZvWShmSpP1o+ow+IhYC7wd+LzMTIDN/mplbquergeeBk4bbPjOXZGZfZvb19PQ0W4YkaQRNBX1EXAB8BrgwM7fVtfdExBuq5ycC04EX2lGoJKk5Iw7dRMRK4GxgYkQMAZ+ndpXNLwMPRQTAE9UVNu8FboyIHcDrwNWZ+WqHapckNaCRq24WDNN8+z7WvRu4u9WitH/7u2po8KbfGcVKJI0FfjJWkgpn0EtS4Qx6SSqcQS9JhTPoJalwBr0kFc6gl6TCGfSSVDiDXpIK19LslTr47OtTs35iVjp0eUYvSYUz6CWpcAa9JBXOoJekwhn0klQ4g16SCtdQ0EfE0ojYFBHr6tqOiYiHIuLfq8ejq/aIiL+KiOciYm1EvKNTxUuSRtbodfTLgC8AK+raFgHfysybImJR9fozwBxq94qdDswGvlQ9qou8vl46dDV0Rp+ZjwB73/t1HrC8er4c+GBd+4qseQI4KiImtaNYSdKBa2WM/vjM3AhQPR5XtU8Gvl+33lDVtoeI6I+IgYgY2Lx5cwtlSJL2pxNvxsYwbfkLDZlLMrMvM/t6eno6UIYkCVoL+ld2DclUj5uq9iFgat16U4CXWziOJKkFrQT9KmBh9XwhcH9d++9XV9+cAWzdNcQjSRp9DV11ExErgbOBiRExBHweuAm4KyKuBF4C5lerPwjMBZ4DtgEfaXPNkqQD0FDQZ+aCfSw6d5h1E/h4K0VJktrHT8ZKUuEMekkqnEEvSYUz6CWpcAa9JBXOoJekwhn0klQ4g16SCmfQS1LhDHpJKpxBL0mFa/RWgiqUtxiUyucZvSQVzqCXpMIZ9JJUuKbH6CPibcCddU0nAn8KHAV8DNh1x+/rM/PBpiuUJLWk6aDPzA3ATICIeAPwn8C91O4odWtmLm5LhZKklrRr6OZc4PnM/I827U+S1CbtCvpLgJV1r6+NiLURsTQijm7TMSRJTWg56CPijcCFwN9WTV8Cfo3asM5G4OZ9bNcfEQMRMbB58+bhVpEktUE7zujnAE9n5isAmflKZr6emTuBrwCnD7dRZi7JzL7M7Ovp6WlDGZKk4bQj6BdQN2wTEZPqln0IWNeGY0iSmtTSFAgRcThwHnBVXfNfRMRMIIHBvZZJkkZZS0GfmduAY/dqu6yliiRJbeUnYyWpcAa9JBXOoJekwhn0klS4om88sq+bakjSoaTooFfzvPOUVA6HbiSpcAa9JBXOoJekwhn0klQ4g16SCmfQS1LhDHpJKpxBL0mFM+glqXAGvSQVruUpECJiEPgx8DqwIzP7IuIY4E6gl9pdpj6cmT9s9ViSpAPXrrluzsnMH9S9XgR8KzNviohF1evPtOlY6iLnwJHGnk4N3cwDllfPlwMf7NBxJEkjaEfQJ/CPEbE6IvqrtuMzcyNA9Xjc3htFRH9EDETEwObNm9tQhiRpOO0YujkrM1+OiOOAhyLie41slJlLgCUAfX192YY6JEnDaPmMPjNfrh43AfcCpwOvRMQkgOpxU6vHkSQ1p6Wgj4gJEXHErufA+cA6YBWwsFptIXB/K8eRJDWv1aGb44F7I2LXvr6emX8fEU8Bd0XElcBLwPwWjyNJalJLQZ+ZLwC/MUz7FuDcVvYtSWoPPxkrSYXz5uBqCz9IJR28PKOXpMIZ9JJUOINekgpn0EtS4Qx6SSqcQS9JhTPoJalwBr0kFc6gl6TCGfSSVDiDXpIKZ9BLUuGc1Ewd5WRnUvc1HfQRMRVYAfwKsBNYkpm3RcQNwMeAXXf8vj4zH2y1UB0a/MUgtV8rZ/Q7gD/KzKer2wmujoiHqmW3Zubi1suTJLWq6aDPzI3Axur5jyNiPTC5XYVJktqjLW/GRkQvMAv4TtV0bUSsjYilEXF0O44hSWpOy2/GRsSbgbuBP8jMH0XEl4A/A7J6vBm4Ypjt+oF+gGnTprVUw77GdSVJLZ7RR8Q4aiH/tcy8ByAzX8nM1zNzJ/AV4PThts3MJZnZl5l9PT09rZQhSdqPpoM+IgK4HVifmbfUtU+qW+1DwLrmy5MktaqVoZuzgMuAf42IZ6q264EFETGT2tDNIHBVSxWqSA63SaOnlatuHgVimEVeMy9JBxGnQJCkwhn0klQ4g16SCuekZlLFeXZUKs/oJalwBr0kFc6gl6TCOUavMeFAP2DluLr0c57RS1LhDHpJKpxBL0mFc4xehxwnVNOhxqCXRuAHqTTWOXQjSYXzjF5FcnhG+jmDXjpIOWSkdulY0EfEBcBtwBuAv8nMmzp1LKkb2vVXg8GtTutI0EfEG4C/Bs4DhoCnImJVZn63E8eT1F7+NVGWTp3Rnw48l5kvAETEHcA8wKCX9nKgfxm08/2HAw1ufwGMTZ0K+snA9+teDwGzO3QsSYVo15xGB+MvpG7WFJnZ/p1GzAd+OzM/Wr2+DDg9Mz9Rt04/0F+9fBuwoe2F7Gki8IMOH6Nb7NvYU2q/wL6Npl/NzJ6RVurUGf0QMLXu9RTg5foVMnMJsKRDx/8FETGQmX2jdbzRZN/GnlL7BfbtYNSpD0w9BUyPiBMi4o3AJcCqDh1LkrQfHTmjz8wdEXEt8A/ULq9cmpnPduJYkqT969h19Jn5IPBgp/bfhFEbJuoC+zb2lNovsG8HnY68GStJOng4qZkkFa64oI+ICyJiQ0Q8FxGLhln+hxHx3YhYGxHfiohf7UadzRipb3XrXRwRGRFj4uqARvoVER+uvm/PRsTXR7vGZjXw8zgtIv4pItZUP5Nzu1HngYqIpRGxKSLW7WN5RMRfVf1eGxHvGO0am9VA336v6tPaiHgsIn5jtGs8YJlZzBe1N36fB04E3gj8C3DqXuucAxxePb8GuLPbdberb9V6RwCPAE8Afd2uu03fs+nAGuDo6vVx3a67jX1bAlxTPT8VGOx23Q327b3AO4B1+1g+F/g7IIAzgO90u+Y29u3Mup/FOWOhb6Wd0e+eeiEz/x+wa+qF3TLznzJzW/XyCWrX+I8FI/at8mfAXwDbR7O4FjTSr48Bf52ZPwTIzE2jXGOzGulbAkdWz9/CXp83OVhl5iPAq/tZZR6wImueAI6KiEmjU11rRupbZj6262eRMZIhpQX9cFMvTN7P+ldSO+sYC0bsW0TMAqZm5jdGs7AWNfI9Owk4KSL+T0Q8Uc2MOhY00rcbgEsjYojaVWqfoAwH+n9xrBoTGVLafPQxTNuwlxVFxKVAH/CbHa2offbbt4j4JeBW4PLRKqhNGvmeHUZt+OZsamdP/zsiTsvM/9vh2lrVSN8WAMsy8+aIeBfw1apvOztfXkc1/H9xrIqIc6gF/bu7XctISjujH3HqBYCI+C3gT4ALM/Ono1Rbq0bq2xHAacDDETFIbVx01Rh4Q7aR79kQcH9m/iwzX6Q2L9L0UaqvFY307UrgLoDMfBwYT20+lbGuof+LY1VEzAD+BpiXmVu6Xc9ISgv6EadeqIY3/ge1kB8rY70wQt8yc2tmTszM3szspTZ2eGFmDnSn3IY1Ml3GfdTeRCciJlIbynlhVKtsTiN9ewk4FyAiTqEW9JtHtcrOWAX8fnX1zRnA1szc2O2i2iEipgH3AJdl5r91u55GFDV0k/uYeiEibgQGMnMV8JfAm4G/jQiAlzLzwq4V3aAG+zbmNNivfwDOj4jvAq8D/30snEU12Lc/Ar4SEZ+iNrRxeVaXcxzMImIltaG0idX7C58HxgFk5pepvd8wF3gO2AZ8pDuVHrgG+vanwLHAF6sM2ZEH+URnfjJWkgpX2tCNJGkvBr0kFc6gl6TCGfSSVDiDXpIKZ9BLUuEMekkqnEEvSYX7/1YNSrFocF3fAAAAAElFTkSuQmCC\n",
"text/plain": [
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"markups = results.compute_markups(costs=costs)\n",
"plt.hist(markups, bins=50);\n",
"plt.legend([\"Markups\"]);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Mergers\n",
"\n",
"Before computing postmerger outputs, we'll supplement our premerger markups with some other outputs. We'll compute HerfindahlHirschman Indices, $\\text{HHI}$, with [`ProblemResults.compute_hhi`](https://pyblp.readthedocs.io/en/stable/_api/pyblp.ProblemResults.compute_hhi.html#pyblp.ProblemResults.compute_hhi); populationnormalized gross expected profits, $\\pi$, with [`ProblemResults.compute_profits`](https://pyblp.readthedocs.io/en/stable/_api/pyblp.ProblemResults.compute_profits.html#pyblp.ProblemResults.compute_profits); and populationnormalized consumer surpluses, $\\text{CS}$, with [`ProblemResults.compute_consumer_surpluses`](https://pyblp.readthedocs.io/en/stable/_api/pyblp.ProblemResults.compute_consumer_surpluses.html#pyblp.ProblemResults.compute_consumer_surpluses)."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"hhi = results.compute_hhi()\n",
"profits = results.compute_profits(costs=costs)\n",
"cs = results.compute_consumer_surpluses()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To compute postmerger outputs, we'll create a new set of firm IDs that represent a merger of firms ``2`` and ``1``."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"product_data['merger_ids'] = product_data['firm_ids'].replace(2, 1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can use [`ProblemResults.compute_approximate_prices`](https://pyblp.readthedocs.io/en/stable/_api/pyblp.ProblemResults.compute_approximate_prices.html#pyblp.ProblemResults.compute_approximate_prices) or [`ProblemResults.compute_prices`](https://pyblp.readthedocs.io/en/stable/_api/pyblp.ProblemResults.compute_prices.html#pyblp.ProblemResults.compute_prices) to estimate postmerger prices. The first method, which is discussed, for example, in [Nevo (1997)](https://pyblp.readthedocs.io/en/stable/references.html#nevo1997), assumes that shares and their price derivatives are unaffected by the merger. The second method does not make these assumptions and iterates over the $\\zeta$markup equation from [Morrow and Skerlos (2011)](https://pyblp.readthedocs.io/en/stable/references.html#morrowandskerlos2011) to solve the full system of $J_t$ equations and $J_t$ unknowns in each market $t$. We'll use the latter, since it is fast enough for this example problem."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"changed_prices = results.compute_prices(\n",
" firm_ids=product_data['merger_ids'], \n",
" costs=costs\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If the problem was configured with more than two columns of firm IDs, we could estimate postmerger prices for the other mergers with the `firms_index` argument, which is by default `1`.\n",
"\n",
"We'll compute postmerger shares with [`ProblemResults.compute_shares`](https://pyblp.readthedocs.io/en/stable/_api/pyblp.ProblemResults.compute_shares.html#pyblp.ProblemResults.compute_shares)."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"changed_shares = results.compute_shares(changed_prices)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Postmerger prices and shares are used to compute other postmerger outputs. For example, $\\text{HHI}$ increases."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAD8CAYAAABXe05zAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAEKNJREFUeJzt3X2QXXV9x/H31yRkKQkSk62Dhu1i61OIawhLJJPo8FAVE4TigMpYmlpxx6FxEqdMjePUCQ5/BMYCdlCmodLYYskwaKZOQOUpPsAU8iAx5AFioOm4hOcUdBsTEvz2j727LnHv7t3k3uz+dt+vmTt77rm/e+73/O7JJ+f87rnnRmYiSSrHG4a7AEnS0BjcklQYg1uSCmNwS1JhDG5JKozBLUmFMbglqTAGtyQVxuCWpMKMb8RCp02blq2trY1YtCSNSps2bXoxM5traduQ4G5tbWXjxo2NWLQkjUoR8T+1tnWoRJIKY3BLUmEMbkkqTEPGuCWNLAcPHqSzs5P9+/cPdyljXlNTE9OnT2fChAlHvAyDWxoDOjs7mTx5Mq2trUTEcJczZmUmL730Ep2dnZx66qlHvJyahkoi4qSIuDMiHo+IHREx94hfUdIxt3//fqZOnWpoD7OIYOrUqUd95FPrHvfXgR9m5iURcRzwR0f1qpKOOUN7ZKjH+zBocEfEicAHgL8GyMxXgVeP+pUlSUeklj3utwEvAP8aEe8FNgFLMvP/GlqZpIZpXXZXXZe3e8XCQdtMmjSJrq6u3vurVq1i48aN3HTTTSxfvpxJkyZx1VVX/b7Gyhf5pk2b9gfP7fHss8+ydOlSNmzYwMSJE2ltbeXGG29kz549fO1rX2Pt2rX1WcERppbgHg/MBj6fmY9ExNeBZcA/9G0UER1AB0BLS0u969QIVi0EavnHPNrYF8dOZnLxxRezaNEiVq9eDcDmzZt57rnnhrmyxqvlw8lOoDMzH6ncv5PuIH+dzFyZme2Z2d7cXNPX7SXpiK1bt44JEybwuc99rnferFmzeP/73w9AV1cXl1xyCe9617v41Kc+RWYC8NWvfpUzzzyTmTNn0tHR0Tv/7LPP5otf/CJz5szhHe94Bz/72c8A2LdvHx//+Mdpa2vjE5/4BO973/t6L+lxzz33MHfuXGbPns2ll17ae1SwbNkyZsyYQVtb2+uOIupl0D3uzHw2In4VEe/MzCeA84Dtda9E0qj229/+llmzZvXe37t3LxdeeGHv/RtuuIHbbrut9/6ePXsGXN7WrVs544wzqj7+6KOPsm3bNt7ylrcwb948HnroIebPn8/ixYv5yle+AsDll1/O2rVr+ehHPwrAoUOHWL9+PXfffTdXX3019913H9/85jeZMmUKW7ZsYevWrb3r8OKLL3LNNddw3333ccIJJ3Dttddy/fXXs3jxYtasWcPjjz9ORPDyyy8PvbMGUetZJZ8HvlM5o+Qp4NN1r0TSqHb88cezefPm3vs9Y9w9vvCFL/zBGPfRmDNnDtOnTwe698R3797N/PnzWbduHddddx379u1j7969nHbaab3B/bGPfQyAM844g927dwPw4IMPsmTJEgBmzpxJW1sbAA8//DDbt29n3rx5ALz66qvMnTuXE088kaamJq644goWLlzIBRdccFTr0Z+agjszNwPtdX91STpCp512GnfeeWfVxydOnNg7PW7cOA4dOsT+/fu58sor2bhxI6eccgrLly9/3TnVPc/paQ/0DqUcLjP54Ac/yO233/4Hj61fv57777+f1atXc9NNN/HAAw8c0TpW47VKJBXp3HPP5cCBA9xyyy298zZs2MBPfvKTqs/pCelp06bR1dU1YPD3mD9/PnfccQcA27dv57HHHgPgrLPO4qGHHmLXrl1A91j4zp076erq4pVXXmHBggXceOONrzvKqBe/8i6NQaPhLJeIYM2aNSxdupQVK1bQ1NTUezrg008/3e9zTjrpJD772c/ynve8h9bWVs4888xBX+fKK69k0aJFtLW1cfrpp9PW1sYb3/hGmpubWbVqFZdddhkHDhwA4JprrmHy5MlcdNFF7N+/n8zkhhtuqOt6A0S1w4Cj0d7env6QwtjhKXC/N1L7YseOHbz73e8e1hpK9dprr3Hw4EGampp48sknOe+889i5cyfHHXfcES+zv/cjIjZlZk1D0u5xS9IA9u3bxznnnMPBgwfJTG6++eajCu16MLglaQCTJ08ecT/F6IeT0hjRiGFRDV093geDWxoDmpqaeOmllwzvYdZzPe6mpqajWo5DJdIYMH36dDo7O3nhhReGu5Qxr+cXcI6GwS2NARMmTDiqX1zRyOJQiSQVxuCWpMIY3JJUGINbkgpjcEtSYQxuSSqMwS1JhTG4JakwBrckFcbglqTCGNySVBiDW5IKY3BLUmEMbkkqjMEtSYUxuCWpMAa3JBWmpl/AiYjdwG+A14BDmdneyKIkSdUN5afLzsnMFxtWiSSpJg6VSFJhat3jTuCeiEjgnzNz5eENIqID6ABoaWmpX4UalVqX3TWk9rtXLGxQJapVtffM9+bYq3WPe15mzgY+AvxtRHzg8AaZuTIz2zOzvbm5ua5FSpJ+r6bgzsw9lb/PA2uAOY0sSpJU3aDBHREnRMTknmngQ8DWRhcmSepfLWPcbwbWRERP+//IzB82tCpJUlWDBndmPgW89xjUIkmqgacDSlJhDG5JKozBLUmFMbglqTAGtyQVxuCWpMIY3JJUGINbkgpjcEtSYQxuSSqMwS1JhTG4JakwBrckFcbglqTCGNySVBiDW5IKY3BLUmEMbkkqjMEtSYUxuCWpMAa3JBXG4JakwhjcklQYg1uSCmNwS1Jhag7uiBgXEY9GxNpGFiRJGthQ9riXADsaVYgkqTY1BXdETAcWAv/S2HIkSYOpdY/7RuDvgd81sBZJUg3GD9YgIi4Ans/MTRFx9gDtOoAOgJaWlroVqMZpXXbXcJdw1Kqtw+4VC+vSXhqJatnjngdcGBG7gdXAuRFx2+GNMnNlZrZnZntzc3Ody5Qk9Rg0uDPzS5k5PTNbgU8CD2TmXza8MklSvzyPW5IKM+gYd1+Z+WPgxw2pRJJUE/e4JakwBrckFcbglqTCGNySVBiDW5IKY3BLUmEMbkkqjMEtSYUxuCWpMAa3JBXG4JakwhjcklQYg1uSCmNwS1JhDG5JKozBLUmFMbglqTAGtyQVxuCWpMIY3JJUGINbkgpjcEtSYQxuSSqMwS1JhTG4JakwBrckFWbQ4I6IpohYHxG/iIhtEXH1sShMktS/8TW0OQCcm5ldETEBeDAifpCZDze4NklSPwYN7sxMoKtyd0Lllo0sSpJUXS173ETEOGAT8GfANzLzkX7adAAdAC0tLfWsUYdpXXbXkNrvXrGwQZWMXEPto5L6tFqtY/F9Hqtq+nAyM1/LzFnAdGBORMzsp83KzGzPzPbm5uZ61ylJqhjSWSWZ+TLwY+D8hlQjSRpULWeVNEfESZXp44E/Bx5vdGGSpP7VMsZ9MvDtyjj3G4A7MnNtY8uSJFVTy1klW4DTj0EtkqQa+M1JSSqMwS1JhTG4JakwBrckFcbglqTCGNySVBiDW5IKY3BLUmEMbkkqjMEtSYUxuCWpMAa3JBXG4JakwhjcklQYg1uSCmNwS1JhDG5JKozBLUmFMbglqTAGtyQVxuCWpMIY3JJUGINbkgpjcEtSYQxuSSrMoMEdEadExLqI2BER2yJiybEoTJLUv/E1tDkE/F1m/jwiJgObIuLezNze4NokSf0YdI87M5/JzJ9Xpn8D7ADe2ujCJEn9G9IYd0S0AqcDjzSiGEnS4GoZKgEgIiYB3wWWZuav+3m8A+gAaGlpqVuBY1nrsruGu4SjUnr9x0I9+6jasnavWFi312jk6w7UF41eh2qGq08HU9Med0RMoDu0v5OZ3+uvTWauzMz2zGxvbm6uZ42SpD5qOaskgG8BOzLz+saXJEkaSC173POAy4FzI2Jz5bagwXVJkqoYdIw7Mx8E4hjUIkmqgd+clKTCGNySVBiDW5IKY3BLUmEMbkkqjMEtSYUxuCWpMAa3JBXG4JakwhjcklQYg1uSCmNwS1JhDG5JKozBLUmFMbglqTAGtyQVxuCWpMIY3JJUGINbkgpjcEtSYQxuSSqMwS1JhTG4JakwBrckFcbglqTCGNySVJhBgzsibo2I5yNi67EoSJI0sFr2uFcB5ze4DklSjQYN7sz8KbD3GNQiSarB+HotKCI6gA6AlpaWI15O67K7htR+94qFR/xajTDU+qHx63AkNWlgQ+3TkfgeNLqmasuv5/Y+XHlxLNZtIHX7cDIzV2Zme2a2Nzc312uxkqTDeFaJJBXG4JakwtRyOuDtwH8B74yIzoj4TOPLkiRVM+iHk5l52bEoRJJUG4dKJKkwBrckFcbglqTCGNySVBiDW5IKY3BLUmEMbkkqjMEtSYUxuCWpMAa3JBXG4JakwhjcklQYg1uSCmNwS1JhDG5JKozBLUmFMbglqTAGtyQVxuCWpMIY3JJUGINbkgpjcEtSYQxuSSqMwS1JhTG4JakwNQV3RJwfEU9ExK6IWNbooiRJ1Q0a3BExDvgG8BFgBnBZRMxodGGSpP7Vssc9B9iVmU9l5qvAauCixpYlSaqmluB+K/CrPvc7K/MkScMgMnPgBhGXAh/OzCsq9y8H5mTm5w9r1wF0VO6+E3ii/uWOKNOAF4e7iBHIfumf/VKdfdPtTzKzuZaG42to0wmc0uf+dGDP4Y0ycyWwsqbyRoGI2JiZ7cNdx0hjv/TPfqnOvhm6WoZKNgBvj4hTI+I44JPA9xtbliSpmkH3uDPzUEQsBn4EjANuzcxtDa9MktSvWoZKyMy7gbsbXEtpxsyw0BDZL/2zX6qzb4Zo0A8nJUkji195l6TCGNx9RMStEfF8RGztM+9NEXFvRPyy8ndKZX5ExD9VLgOwJSJm93nOokr7X0bEouFYl3qq0i/LI+LpiNhcuS3o89iXKv3yRER8uM/8UXXphIg4JSLWRcSOiNgWEUsq88f0NjNAv4z5baZuMtNb5QZ8AJgNbO0z7zpgWWV6GXBtZXoB8AMggLOARyrz3wQ8Vfk7pTI9ZbjXrQH9shy4qp+2M4BfABOBU4En6f5Qe1xl+m3AcZU2M4Z73Y6yX04GZlemJwM7K+s/preZAfplzG8z9bq5x91HZv4U2HvY7IuAb1emvw38RZ/5/5bdHgZOioiTgQ8D92bm3sz8X+Be4PzGV984VfqlmouA1Zl5IDP/G9hF92UTRt2lEzLzmcz8eWX6N8AOur9VPKa3mQH6pZoxs83Ui8E9uDdn5jPQvUECf1yZX+1SAGPpEgGLK4f8t/YMBzBG+yUiWoHTgUdwm+l1WL+A20xdGNxHLvqZlwPMH21uBv4UmAU8A/xjZf6Y65eImAR8F1iamb8eqGk/80Zt3/TTL24zdWJwD+65yuEslb/PV+ZXuxRATZcIKF1mPpeZr2Xm74Bb6D6shTHWLxExge5w+k5mfq8ye8xvM/31i9tM/Rjcg/s+0PMp/yLgP/vM/6vKmQJnAa9UDot/BHwoIqZUDgU/VJk3qvQEU8XFQM8ZJ98HPhkREyPiVODtwHpG4aUTIiKAbwE7MvP6Pg+N6W2mWr+4zdTRcH86OpJuwO10H8IdpPt/+88AU4H7gV9W/r6p0jbo/oGJJ4HHgPY+y/kbuj9g2QV8erjXq0H98u+V9d5C9z+mk/u0/3KlX54APtJn/gK6zzB4EvjycK9XHfplPt2H7luAzZXbgrG+zQzQL2N+m6nXzW9OSlJhHCqRpMIY3JJUGINbkgpjcEtSYQxuSSqMwS1JhTG4JakwBrckFeb/AYDW1SQB+8IsAAAAAElFTkSuQmCC\n",
"text/plain": [
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"changed_cs = results.compute_consumer_surpluses(changed_prices)\n",
"plt.hist(changed_cs  cs, bins=50);\n",
"plt.legend([\"Consumer Surplus Changes\"]);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Bootstrapping Results\n",
"\n",
"Postestimation outputs can be informative, but they don't mean much without a sense sampletosample variability. One way to estimate confidence intervals for postestimation outputs is with a standard bootstrap procedure:\n",
"\n",
"1. Construct a large number of bootstrap samples by sampling with replacement from the original product data.\n",
"2. Initialize and solve a [`Problem`](https://pyblp.readthedocs.io/en/stable/_api/pyblp.Problem.html#pyblp.Problem) for each bootstrap sample.\n",
"3. Compute the desired postestimation output for each bootstrapped [`ProblemResults`](https://pyblp.readthedocs.io/en/stable/_api/pyblp.ProblemResults.html#pyblp.ProblemResults) and from the resulting empirical distribution, construct boostrap confidence intervals.\n",
"\n",
"Although appealing because of its simplicity, the computational resources required for this procedure are often prohibatively expensive. Furthermore, human oversight of the optimization routine is often required to determine whether the routine ran into any problems and if it successfully converged. Human oversight of estimation for each bootstrapped problem is usually not feasible.\n",
"\n",
"A more reasonable alternative is a parametric bootstrap procedure:\n",
"\n",
"1. Construct a large number of draws from the estimated joint distribution of parameters.\n",
"2. Compute the implied mean utility, $\\delta$, and shares, $s$, for each draw. If a supply side was estimated, also computed the implied marginal costs, $c$, and prices, $p$.\n",
"3. Compute the desired postestimation output under each of these parametric bootstrap samples. Again, from the resulting empirical distribution, construct boostrap confidence intervals.\n",
"\n",
"Compared to the standard bootstrap procedure, the parametric bootstrap requires far fewer computational resources, and is simple enough to not require human oversight of each bootstrap iteration. The primary complication to this procedure is that when supply is estimated, equilibrium prices and shares need to be computed for each parametric bootstrap sample by iterating over the $\\zeta$markup equation from [Morrow and Skerlos (2011)](https://pyblp.readthedocs.io/en/stable/references.html#morrowandskerlos2011). Although nontrivial, this fixed point iteration problem is much less demanding than the full optimization routine required to solve the BLP problem from the start.\n",
"\n",
"An empirical distribution of results computed according to this parametric bootstrap procedure can be created with the [`ProblemResults.bootstrap`](https://pyblp.readthedocs.io/en/stable/_api/pyblp.ProblemResults.bootstrap.html#pyblp.ProblemResults.bootstrap) method, which returns a [`BootstrappedResults`](https://pyblp.readthedocs.io/en/stable/_api/pyblp.BootstrappedResults.html#pyblp.BootstrappedResults) class that can be used just like [`ProblemResults`](https://pyblp.readthedocs.io/en/stable/_api/pyblp.ProblemResults.html#pyblp.ProblemResults) to compute various postestimation outputs. The difference is that [`BootstrappedResults`](https://pyblp.readthedocs.io/en/stable/_api/pyblp.BootstrappedResults.html#pyblp.BootstrappedResults) methods return arrays with an extra first dimension, along which bootstrapped results are stacked.\n",
"\n",
"We'll construct 90% parametric bootstrap confidence intervals for estimated mean own elasticities in each market of the fake cereal problem. Usually, bootstrapped confidence intervals should be based on thousands of draws, but we'll only use a few for the sake of speed in this example."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Bootstrapped Problem Results Summary:\n",
"======================\n",
"Computation Bootstrap\n",
" Time Draws \n",
" \n",
" 00:00:24 100 \n",
"======================"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"bootstrapped_results = results.bootstrap(draws=100, seed=0)\n",
"bootstrapped_results"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
"
],
"text/plain": [
" market_ids weights nodes0 nodes1 nodes2 nodes3 nodes4 \\\n",
"0 1971 0.005 0.548814 0.457760 0.564690 0.395537 0.392173 \n",
"1 1971 0.005 0.715189 0.376918 0.839746 0.844017 0.041157 \n",
"2 1971 0.005 0.602763 0.702335 0.376884 0.150442 0.923301 \n",
"3 1971 0.005 0.544883 0.207324 0.499676 0.306309 0.406235 \n",
"4 1971 0.005 0.423655 0.074280 0.081302 0.094570 0.944282 \n",
"\n",
" income \n",
"0 9.728478 \n",
"1 7.908957 \n",
"2 11.079404 \n",
"3 17.641671 \n",
"4 12.423995 "
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"blp_agent_data.head()"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/xpython",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.6"
}
},
"nbformat": 4,
"nbformat_minor": 2
}PK+N₯dέrrKpyblpstable/_downloads/5b262cb63645dabfa5773dc21f5fa2b9/logit_nested.ipynb{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Logit and Nested Logit Tutorial"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"'0.7.0'"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import pyblp\n",
"import numpy as np\n",
"import pandas as pd\n",
"\n",
"pyblp.options.digits = 2\n",
"pyblp.options.verbose = False\n",
"pyblp.__version__"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this tutorial, we'll use data from [Nevo (2000)](https://pyblp.readthedocs.io/en/stable/references.html#nevo2000) to solve the paper's fake cereal problem. Locations of CSV files that contain the data are in the [`data`](https://pyblp.readthedocs.io/en/stable/_api/pyblp.data.html#modulepyblp.data) module.\n",
"\n",
"We will compare two simple models, the plain (IIA) logit model and the nested logit (GEV) model using the fake cereal dataset of [Nevo (2000)](https://pyblp.readthedocs.io/en/stable/references.html#nevo2000).\n",
"\n",
"## Theory of Plain Logit\n",
"\n",
"Let's start with the plain logit model under independence of irrelevant alternatives (IIA). In this model (indirect) utility is given by\n",
"\n",
"$$U_{jti} = \\alpha p_{jt} + x_{jt} \\beta^x + \\xi_{jt} + \\epsilon_{jti},$$\n",
"\n",
"where $\\epsilon_{jti}$ is distributed IID with the Type I Extreme Value (Gumbel) distribution. It is common to normalize the mean utility of the outside good to zero so that $U_{0ti} = \\epsilon_{0ti}$. This gives us aggregate marketshares\n",
"\n",
"$$s_{jt} = \\frac{\\exp(\\alpha p_{jt} + x_{jt} \\beta^x + \\xi_{jt})}{1 + \\sum_k \\exp(\\alpha p_{jt} + x_{kt} \\beta^x + \\xi_{kt})}.$$\n",
"\n",
"If we take logs we get\n",
"\n",
"$$\\log s_{jt} = \\alpha p_{jt} + x_{jt} \\beta^x + \\xi_{jt}  0  \\log \\sum_k \\exp(\\alpha p_{jt} + x_{kt} \\beta^x + \\xi_{kt})$$\n",
"\n",
"and\n",
"\n",
"$$\\log s_{0t} = 0  \\log \\sum_k \\exp(\\alpha p_{jt} + x_{kt} \\beta^x + \\xi_{kt}).$$\n",
"\n",
"By differencing the above we get a linear estimating equation:\n",
"\n",
"$$\\log s_{jt}  \\log s_{0t} = \\alpha p_{jt} + x_{jt} \\beta^x + \\xi_{jt}.$$\n",
"\n",
"Because the left hand side is data, we can estimate this model using linear IV GMM."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Application of Plain Logit\n",
"\n",
"A Logit [`Problem`](https://pyblp.readthedocs.io/en/stable/_api/pyblp.Problem.html#pyblp.Problem) can be created by simply excluding the formulation for the nonlinear parameters, $X_2$, along with any agent information. In other words, it requires only specifying the _linear component_ of demand.\n",
"\n",
"We'll set up and solve a simple version of the fake data cereal problem from [Nevo (2000)](https://pyblp.readthedocs.io/en/stable/references.html#nevo2000). Since we won't include any nonlinear characteristics or parameters, we don't have to worry about configuring an optimization routine.\n",
"\n",
"There are some important reserved variable names:\n",
"\n",
" `market_ids` are the unique market identifiers which we subscript with $t$.\n",
" `shares` specifies the marketshares which need to be between zero and one, and within a market ID, $\\sum_{j} s_{jt} \\leq 1$.\n",
" `prices` are prices $p_{jt}$. These have some special properties and are _always_ treated as endogenous.\n",
" `demand_instruments0`, `demand_instruments1`, and so on are numbered demand instruments. These represent only the _excluded_ instruments. The exogenous regressors in $X_1$ will be automatically added to the set of instruments.\n",
"\n",
"We begin with two steps:\n",
"\n",
"1. Load the product data which at a minimum consists of `market_ids`, `shares`, `prices`, and at least a single column of demand instruments, `demand_instruments0`.\n",
"2. Define a [`Formulation`](https://pyblp.readthedocs.io/en/stable/_api/pyblp.Formulation.html#pyblp.Formulation) for the $X_1$ (linear) demand model.\n",
"\n",
"  This and all other formulas are similar to R and [patsy](https://patsy.readthedocs.io/en/stable/) formulas.\n",
"  It includes a constant by default. To exclude the constant, specify either a `0` or a `1`.\n",
"  To efficiently include fixed effects, use the `absorb` option and specify which categorical variables you would like to absorb.\n",
"  Some model reduction may happen automatically. The constant will be excluded if you include fixed effects and some precautions are taken against collinearity. However, you will have to make sure that differentlynamed variables are not collinear.\n",
" \n",
"3. Combine the [`Formulation`](https://pyblp.readthedocs.io/en/stable/_api/pyblp.Formulation.html#pyblp.Formulation) and product data to construct a [`Problem`](https://pyblp.readthedocs.io/en/stable/_api/pyblp.Problem.html#pyblp.Problem).\n",
"4. Use [`Problem.solve`](https://pyblp.readthedocs.io/en/stable/_api/pyblp.Problem.solve.html#pyblp.Problem.solve) to estimate paramters."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Loading the Data\n",
"\n",
"The `product_data` argument of [`Problem`](https://pyblp.readthedocs.io/en/stable/_api/pyblp.Problem.html#pyblp.Problem) should be a structured arraylike object with fields that store data. Product data can be a structured [NumPy](https://www.numpy.org/) array, a [pandas](https://pandas.pydata.org/) DataFrame, or other similar objects."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"