jeffgortmaker / pyblp

BLP Demand Estimation with Python
https://pyblp.readthedocs.io
MIT License
228 stars 82 forks source link

Beta Estimates for Unobserved Characteristics #115

Closed jessi96 closed 2 years ago

jessi96 commented 2 years ago

Hi,

thanks a lot for the great work on the PyBLP package. At the moment, I am working on the model and try to compute the coefficients and standard errors for the constant + the unobserved product characteristics. Contrary to the solution of the thread cited below, I do not use product fixed effects, but rather brand fixed effects. Therefore, I do not have a single number that I can align to each brand regarding the unobserved characteristics. Is there still a way to compute the constant + unobserved characteristics coefficients, e.g. take the mean / median value of the unobserved characteristic (e.g. sugar) per brand?

You should be able to run this just before the "Restricting Parameters" section:

# solve the problem with dummies instead of absorbing the product fixed effects
X1_formulation_with_d = pyblp.Formulation('0 + prices + C(product_ids)')
product_formulations_with_d = (X1_formulation_with_d, X2_formulation)
nevo_problem_with_d = pyblp.Problem(
    product_formulations_with_d,
    product_data,
    agent_formulation,
    agent_data
)
nevo_results_with_d = nevo_problem_with_d.solve(
    initial_sigma,
    initial_pi,
    optimization=pyblp.Optimization('bfgs', {'gtol': 0.5}),  # use a loose tolerance as in the original paper
    method='1s'
)

# collect the estimated fixed effects, their covariances, and the desired regressors
d = nevo_results_with_d.beta[1:]
V = nevo_results_with_d.parameter_covariances[-d.size:, -d.size:]
X = np.zeros((d.size, 3))
for i, label in enumerate(nevo_results_with_d.beta_labels[1:]):
    row = product_data[product_data['product_ids'] == label.split("'")[1]].iloc[0]
    X[i] = np.r_[1, row[['sugar', 'mushy']]]

# form the minimum distance estimates and their standard errors
W = np.linalg.inv(V)
beta = np.linalg.solve(X.T @ W @ X, X.T @ W @ d)
beta_se = np.sqrt(np.diag(np.linalg.inv(X.T @ W @ X)) / nevo_problem.N)
print(np.c_[beta, beta_se])

The above approximately replicates the original paper, giving:

[[-1.85991326  0.25721032]
 [ 0.14225657  0.01287912]
 [ 0.75831534  0.19939801]]

The point estimates and the standard errors for the last two covariates are a bit off, perhaps because the results are somewhat sensitive to the optimizer termination tolerance. (Above, I set a very loose tolerance to get a similar local (not global) optimum that was reported in the table you posted above.) Or maybe the reported estimates used a slightly different standard error formula -- not sure.

Hope this helps. I'll keep this issue open for now as a reminder to incorporate something like this into the notebook.

Originally posted by @jeffgortmaker in https://github.com/jeffgortmaker/pyblp/issues/93#issuecomment-888527383

jeffgortmaker commented 2 years ago

So to be clear, the cited solution computes coefficients on observed characteristics: a constant, sugar, and mushy. I'm guessing you're referring to them as unobserved characteristics because we're assuming there is unobserved preference heterogeneity for these same characteristics by including them in our X2 formula. Technically, the only unobserved characteristic in the model is xi. Hopefully that clears up some terminology!

The above solution uses a minimum distance procedure because if we tried to include a constant, sugar, and mushy in our X1 formula, these would be perfectly collinear with the product fixed effects, so we couldn't estimate their coefficients. Based on your description, I don't think you have the same issue. If there's characteristic variation within brand, you should be able to identify coefficients on these characteristics from this within-brand variation. Does that make sense?

jeffgortmaker commented 2 years ago

I'm going to close this for now, but feel free to keep commenting / re-open the issue if there's anything else.