Problem set 2 model solutions¶

ECON-L1300

Helena Rantakaulio, 29.2.2024

Exercise 1¶

Note that the actual code displayed in the html file is Matlab code.

In [14]:
from IPython.display import HTML
HTML(filename='contraction_mapping.html')
Out[14]:
Untitled
%ECON-L1300
%Problem set 2
%Exercise 1
%Author: Helena Rantakaulio
 
% Should use from the data:
% roof
% true market shares
% Should use from the model:
% true theta_2
 
clear;
boat_data=readmatrix("boat_data.xlsx");
roof = boat_data(:,6); %The variable with the random coef 4000x1
s = boat_data(:,2); %Market shares, is 4000x1
 
NS = 1000; %Number of draws/simulations
N = 4000; %total number of products x markets
Tol = 1*10^(-12); %Set the tolerance
Err = 1; %Initiate the error between new and old delta
 
%Set the initial guess for each delta as 0. We need N deltas.
de_1 = zeros(N,1);
 
%set the guess of theta_2 as the true value of theta_2
%the real theta2 = sigma = sqrt(sigma^2) = sqrt(16)
theta_2 = sqrt(16);
 
%Simulate \tilde{\beta_i}'s
%by drawing NS observations from the standard normal distribution and
%multiplying them by theta2
rng(0); %set seed as 0
beta_i = normrnd(0,1,NS,1)*theta_2; %NSx1 matrix
 
%Is exactly the same as
rng(0);
beta_i2 = normrnd(0,theta_2,NS,1);
 
%Calculate roof * \tilde{\beta_i} for each simulated value
%Each column contains all values of roof multiplied with one simulated value
%so one column corresponds with one simulated beta_i.
%When we iterate over the simulated beta_i's, we'll iterate over the
%columns of roof_i.
roof_mat = repmat(roof,1,NS); %4000x1000
beta_i_mat = repmat(beta_i',N,1); %4000x1000
roof_i = roof_mat.*beta_i_mat; %4000x1000
 
%Solve the contraction mapping
while (Err>=Tol)
 
de=de_1; %current guess of mean utility
s_pred = zeros(N,1); %matrix for predicted market shares
s_num = zeros(N,1); %matrix for the numerators in market shares
 
%Calculate the numerical integral:
 
for i=1:1:NS %iterate over number of draws for beta_i
%calculate the numerators
s_num= exp(de + roof_i(:,i)); %4000x1
 
%calculate the denominator for each market
s_den= reshape(s_num, 4, 1000); %each column is one market
s_den= sum(s_den); %sum by column, 1x1000
s_den= repmat(s_den, 4, 1); %4x1000
s_den= reshape(s_den, 4000, 1); %4000x1
 
%calculate the value inside the integral for each jt
s_value = s_num./(1+s_den); %4000x1
%add this to the sum of values evaluated for each jt using
%different simulated beta_i's
s_pred = s_pred + s_value; %4000x1
 
end
 
%To calculate the integral, take the mean of the sum of the logit
%market shares. You simulated NS values, so divide the sum by NS.
s_pred = s_pred/NS; %4000x1
 
%Update the guess for delta
de_1 = de + log(s) - log(s_pred);
 
%Calculate the maximum difference between deltas.
%The while loop continues until the maximum difference if within the
%tolerance.
Err = max(abs(de_1-de));
%disp(Err);
end
 
delta=de_1;
disp("Convergence achieved!");
Convergence achieved!
 
%Calculate means and standard deviations by product over all markets
d_means = mean((reshape(delta, 4, 1000))');
d_stdev = std((reshape(delta, 4, 1000))');
 
disp(d_means);
0.9179 0.7010 -0.8059 -0.8478
disp(d_stdev);
1.5866 1.6082 1.7084 1.7068
 
 

Note: when calculating the error, I followed K. Sudhir's example (Yale school of management) and took the maximum of the absolute values of the differences in mean utilities. Conlon and Gortmaker (2020) take the norm of vector subtraction.

Exercise 3¶

You were asked to estimate the model using the exogenous demand shifters and exogenous cost shifters as instruments. The possible instruments include

  • the firm's own cost shifter
  • competitors' cost shifters
  • competitors' qualities

Note that you also have many options when using the competitors' characteristics as instruments: You could take a sum of the characteristics of all competitors or you could take the closest competitor's characteristic as a single instrument and the sum of the other competitors' characteristics as another instrument. The "within-nest" instruments may be stronger and you could also just use those as instruments. Let's try this using the following instruments:

  • the firm's own cost shifter
  • the closest competitor's cost shifter (i.e. for a boat that has a roof this is the cost shifter of the other boat that has a roof)
  • the closest competitor's quality
  • sum of the other competitors' cost shifters
  • sum of the other competitors' qualities
In [1]:
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import pyblp
pyblp.options.verbose = True
pyblp.options.digits = 4
from tabulate import tabulate
import time

product_data = pd.read_excel('boat_data.xlsx')

Note that pyBLP recognizes the demand instruments by the variable names. So, you should name your instruments demand_instruments0, demand_instruments1 and so on.

In [2]:
# Own cost shifter
product_data = product_data.rename(columns={'cost_shifter': 'demand_instruments0'})

# Closest competitor's cost shifter
product_data['demand_instruments1']= product_data.groupby(['market_ids', 'roof'])['demand_instruments0'].transform("sum") - product_data['demand_instruments0']

# Closest competitor's quality
product_data['demand_instruments2']= product_data.groupby(['market_ids', 'roof'])['quality'].transform("sum") - product_data['quality']

# Sum of the other competitors' cost shifters
product_data['demand_instruments3']= product_data.groupby(['market_ids'])['demand_instruments0'].transform("sum") - product_data['demand_instruments0'] - product_data['demand_instruments1']

# Sum of the other competitors' qualities
product_data['demand_instruments4']= product_data.groupby(['market_ids'])['quality'].transform("sum") - product_data['quality'] - product_data['demand_instruments2']
In [3]:
product_data.head()
Out[3]:
prices shares length quality demand_instruments0 roof firm_ids market_ids demand_instruments1 demand_instruments2 demand_instruments3 demand_instruments4
0 5.858942 0.498415 9.336710 1.229133 0.802284 1 1 1 1.520421 0.752459 1.094955 1.098638
1 4.600115 0.102175 7.092352 0.752459 1.520421 1 2 1 0.802284 1.229133 1.094955 1.098638
2 1.000218 0.004851 6.159682 0.341548 0.334998 0 3 1 0.759957 0.757090 2.322705 1.981593
3 3.233179 0.216493 5.780847 0.757090 0.759957 0 4 1 0.334998 0.341548 2.322705 1.981593
4 3.827983 0.354432 9.336710 0.573911 0.164722 1 1 2 0.762669 0.879428 1.483692 2.171100
In [4]:
X1_formulation = pyblp.Formulation('0 + quality + length + prices + roof') #Linear demand parameters
X2_formulation = pyblp.Formulation('0 + roof') #Nonlinear demand model, i.e. random coef variabels. 0 for no random coefficient on the constant
product_formulations = (X1_formulation, X2_formulation)

n_sim = 5000
mc_integration = pyblp.Integration('monte_carlo', size = n_sim, specification_options={'seed': 0}) ## Monte carlo simulation for integration in the inner loop
problem = pyblp.Problem(product_formulations, product_data, integration = mc_integration) ## Formulate problem: combine Formulation, product_data, and the Integration configuration

Initializing the problem ...
Initialized the problem after 00:00:02.

Dimensions:
==========================================
 T     N     F      I      K1    K2    MD 
----  ----  ---  -------  ----  ----  ----
1000  4000   4   5000000   4     1     8  
==========================================

Formulations:
============================================================
       Column Indices:            0       1       2      3  
-----------------------------  -------  ------  ------  ----
 X1: Linear Characteristics    quality  length  prices  roof
X2: Nonlinear Characteristics   roof                        
============================================================

Estimate the model using the exogenous cost and demand shifters as instruments (not optimal instruments)¶

In [5]:
bfgs = pyblp.Optimization('bfgs', {'gtol': 1e-4}) ## Use Broyden–Fletcher–Goldfarb–Shanno algorithm within PyBlp
initial_sigma = np.ones((1,1)) ## An initial guess for the nonlinear parameters 
#Note that as meantioned in the tutorial, this is a simpler, non-default optimization configuration 
#and the tolerance is relatively loose.
In [6]:
# Solve in parallel

with pyblp.parallel(6):
    results = problem.solve(
    initial_sigma, # Initial sigmas
    optimization=bfgs, ## Artelys Knitro is best practice but commercial
    method='1s' ## 1 Step GMM, it should be okay if you didn't specify this
)

Starting a pool of 6 processes ...
Started the process pool after 00:00:00.
Solving the problem ...

Nonlinear Coefficient Initial Values:
==================
Sigma:     roof   
------  ----------
 roof   +1.000E+00
==================
Starting optimization ...

GMM   Computation  Optimization   Objective   Fixed Point  Contraction  Clipped  Objective    Objective    Gradient             
Step     Time       Iterations   Evaluations  Iterations   Evaluations  Shares     Value     Improvement     Norm       Theta   
----  -----------  ------------  -----------  -----------  -----------  -------  ----------  -----------  ----------  ----------
 1     00:03:48         0             1          7504         23273        0     +1.933E+02               +8.415E+01  +1.000E+00
 1     00:03:21         0             2          8837         27088        0     +1.041E+02  +8.913E+01   +8.319E+01  +2.010E+00
 1     00:03:10         0             3          10022        30705        0     +1.249E+02               +1.127E+02  +6.050E+00
 1     00:03:26         0             4          9574         29313        0     +1.214E+01  +9.200E+01   +2.271E+00  +3.968E+00
 1     00:03:51         1             5          9631         29488        0     +1.213E+01  +1.236E-02   +1.977E+00  +4.050E+00
 1     00:03:16         2             6          9605         29408        0     +1.211E+01  +2.084E-02   +1.322E+00  +4.037E+00
 1     00:04:41         3             7          9581         29358        0     +1.209E+01  +1.690E-02   +2.621E-03  +4.012E+00
 1     00:04:20         4             8          9578         29362        0     +1.209E+01  +6.654E-08   +3.522E-06  +4.012E+00

Optimization completed after 00:29:53.
Computing the Hessian and estimating standard errors ...
Computed results after 00:09:12.

Problem Results Summary:
======================================================================================
GMM   Objective    Gradient               Clipped  Weighting Matrix  Covariance Matrix
Step    Value        Norm      Hessian    Shares   Condition Number  Condition Number 
----  ----------  ----------  ----------  -------  ----------------  -----------------
 1    +1.209E+01  +3.522E-06  +5.162E+01     0        +4.491E+02        +1.162E+04    
======================================================================================

Cumulative Statistics:
===========================================================================
Computation  Optimizer  Optimization   Objective   Fixed Point  Contraction
   Time      Converged   Iterations   Evaluations  Iterations   Evaluations
-----------  ---------  ------------  -----------  -----------  -----------
 00:39:05       Yes          5             9          83910       257357   
===========================================================================

Nonlinear Coefficient Estimates (Robust SEs in Parentheses):
====================
Sigma:      roof    
------  ------------
 roof    +4.012E+00 
        (+4.019E-01)
====================

Beta Estimates (Robust SEs in Parentheses):
======================================================
  quality        length        prices         roof    
------------  ------------  ------------  ------------
 +2.008E+00    +5.285E-01    -2.042E+00    +3.812E+00 
(+7.854E-02)  (+3.284E-02)  (+8.561E-02)  (+1.432E-01)
======================================================
Terminating the pool of 6 processes ...
Terminated the process pool after 00:00:00.
In [7]:
results
Out[7]:
Problem Results Summary:
======================================================================================
GMM   Objective    Gradient               Clipped  Weighting Matrix  Covariance Matrix
Step    Value        Norm      Hessian    Shares   Condition Number  Condition Number 
----  ----------  ----------  ----------  -------  ----------------  -----------------
 1    +1.209E+01  +3.522E-06  +5.162E+01     0        +4.491E+02        +1.162E+04    
======================================================================================

Cumulative Statistics:
===========================================================================
Computation  Optimizer  Optimization   Objective   Fixed Point  Contraction
   Time      Converged   Iterations   Evaluations  Iterations   Evaluations
-----------  ---------  ------------  -----------  -----------  -----------
 00:39:05       Yes          5             9          83910       257357   
===========================================================================

Nonlinear Coefficient Estimates (Robust SEs in Parentheses):
====================
Sigma:      roof    
------  ------------
 roof    +4.012E+00 
        (+4.019E-01)
====================

Beta Estimates (Robust SEs in Parentheses):
======================================================
  quality        length        prices         roof    
------------  ------------  ------------  ------------
 +2.008E+00    +5.285E-01    -2.042E+00    +3.812E+00 
(+7.854E-02)  (+3.284E-02)  (+8.561E-02)  (+1.432E-01)
======================================================

Notice that the estimates are very close to the true parameter values

  • prices -2
  • length 0.5
  • quality 2
  • mean of roof 4 (beta estimates)
  • sigma for roof 4 (nonlinear coefficient estimates)

Estimate the demand model constructing approximations to the optimal instruments¶

In [8]:
#Optimal instruments
updated_problem = results.compute_optimal_instruments(method='approximate').to_problem()
with pyblp.parallel(6):
    updated_results = updated_problem.solve(
        results.sigma,
        results.pi,
        optimization=bfgs,
        method='1s'
)

Computing optimal instruments for theta ...
Computed optimal instruments after 00:00:03.

Optimal Instrument Results Summary:
=======================
Computation  Error Term
   Time        Draws   
-----------  ----------
 00:00:03        1     
=======================
Re-creating the problem ...
Re-created the problem after 00:00:00.

Dimensions:
==========================================
 T     N     F      I      K1    K2    MD 
----  ----  ---  -------  ----  ----  ----
1000  4000   4   5000000   4     1     5  
==========================================

Formulations:
============================================================
       Column Indices:            0       1       2      3  
-----------------------------  -------  ------  ------  ----
 X1: Linear Characteristics    quality  length  prices  roof
X2: Nonlinear Characteristics   roof                        
============================================================
Starting a pool of 6 processes ...
Started the process pool after 00:00:00.
Solving the problem ...

Nonlinear Coefficient Initial Values:
==================
Sigma:     roof   
------  ----------
 roof   +4.012E+00
==================
Starting optimization ...

GMM   Computation  Optimization   Objective   Fixed Point  Contraction  Clipped  Objective    Objective    Gradient             
Step     Time       Iterations   Evaluations  Iterations   Evaluations  Shares     Value     Improvement     Norm       Theta   
----  -----------  ------------  -----------  -----------  -----------  -------  ----------  -----------  ----------  ----------
 1     00:02:56         0             1          9578         29362        0     +3.220E-04               +3.200E-01  +4.012E+00
 1     00:02:59         0             2          9456         28925        0     +7.941E+00               +4.963E+01  +3.692E+00
 1     00:02:53         0             3          9578         29361        0     +1.689E-10  +3.220E-04   +2.317E-04  +4.010E+00
 1     00:02:52         1             4          9590         29388        0     +1.902E-18  +1.689E-10   +2.460E-08  +4.010E+00

Optimization completed after 00:11:40.
Computing the Hessian and estimating standard errors ...
Computed results after 00:08:41.

Problem Results Summary:
======================================================================================
GMM   Objective    Gradient               Clipped  Weighting Matrix  Covariance Matrix
Step    Value        Norm      Hessian    Shares   Condition Number  Condition Number 
----  ----------  ----------  ----------  -------  ----------------  -----------------
 1    +1.902E-18  +2.460E-08  +1.590E+02     0        +4.324E+04        +3.907E+03    
======================================================================================

Cumulative Statistics:
===========================================================================
Computation  Optimizer  Optimization   Objective   Fixed Point  Contraction
   Time      Converged   Iterations   Evaluations  Iterations   Evaluations
-----------  ---------  ------------  -----------  -----------  -----------
 00:20:21       Yes          2             5          47792       146424   
===========================================================================

Nonlinear Coefficient Estimates (Robust SEs in Parentheses):
====================
Sigma:      roof    
------  ------------
 roof    +4.010E+00 
        (+2.296E-01)
====================

Beta Estimates (Robust SEs in Parentheses):
======================================================
  quality        length        prices         roof    
------------  ------------  ------------  ------------
 +2.007E+00    +5.284E-01    -2.042E+00    +3.812E+00 
(+6.506E-02)  (+2.911E-02)  (+7.181E-02)  (+1.153E-01)
======================================================
Terminating the pool of 6 processes ...
Terminated the process pool after 00:00:00.
In [9]:
updated_results
Out[9]:
Problem Results Summary:
======================================================================================
GMM   Objective    Gradient               Clipped  Weighting Matrix  Covariance Matrix
Step    Value        Norm      Hessian    Shares   Condition Number  Condition Number 
----  ----------  ----------  ----------  -------  ----------------  -----------------
 1    +1.902E-18  +2.460E-08  +1.590E+02     0        +4.324E+04        +3.907E+03    
======================================================================================

Cumulative Statistics:
===========================================================================
Computation  Optimizer  Optimization   Objective   Fixed Point  Contraction
   Time      Converged   Iterations   Evaluations  Iterations   Evaluations
-----------  ---------  ------------  -----------  -----------  -----------
 00:20:21       Yes          2             5          47792       146424   
===========================================================================

Nonlinear Coefficient Estimates (Robust SEs in Parentheses):
====================
Sigma:      roof    
------  ------------
 roof    +4.010E+00 
        (+2.296E-01)
====================

Beta Estimates (Robust SEs in Parentheses):
======================================================
  quality        length        prices         roof    
------------  ------------  ------------  ------------
 +2.007E+00    +5.284E-01    -2.042E+00    +3.812E+00 
(+6.506E-02)  (+2.911E-02)  (+7.181E-02)  (+1.153E-01)
======================================================

Note that the estimates are very close to the estimates without approximations to the optimal instruments but the standard errors are smaller especially for the random coefficient. With different sets of instruments the point estimates may differ more.

Problem 4¶

In [10]:
elasticities_own = updated_results.compute_elasticities().reshape((1000,4,4)).mean(axis=0).diagonal()
elasticities_own

Computing elasticities with respect to prices ...
Finished after 00:00:02.
Out[10]:
array([-5.15695222, -5.16483513, -4.29922201, -4.17912948])

The logit elasticities in Problem set 1 were:

  • Firm 1: -5.600759
  • Firm 2: -5.289546
  • Firm 3: -4.072071
  • Firm 4: -3.934288

The nested logit elasticities in Problem set 1 were:

  • Firm 1: -5.695437
  • Firm 2: -5.484764
  • Firm 3: -4.065844
  • Firm 4: -3.933168

The BLP elasticities for firms 1 and 2 are smaller (closer to 0) and for firms 3 and 4 are larger (further from 0).

Let's compare in percentages:

In [11]:
elasticities_ps1_logit = np.array([-5.600759, -5.289546, -4.072071, -3.934288])
elasticities_ps1_nlogit = np.array([-5.695437, -5.484764, -4.065844, -3.933168])

print("Differences in % w.r.t to logit: ", (elasticities_own-elasticities_ps1_logit)*100/elasticities_own)
print("Differences in % w.r.t to nested logit: ", (elasticities_own-elasticities_ps1_nlogit)*100/elasticities_own)

Differences in % w.r.t to logit:  [-8.60598986 -2.41461483  5.28353748  5.85867184]
Differences in % w.r.t to nested logit:  [-10.44191923  -6.19436762   5.42837764   5.88547168]

Problem 5¶

In [12]:
diversions = updated_results.compute_diversion_ratios().reshape((1000,4,4)).mean(axis=0)
print(diversions[0,2], "and" ,diversions[2,0])

Computing diversion ratios with respect to prices ...
Finished after 00:00:02.

0.10439516533920497 and 0.15222955093717164

The diversion ratios calculated in Problem set 1 were

  • logit D_13: 0.1921416
  • logit D_31: 0.3445739
  • nested logit D_13: 0.09510307
  • nested logit D_31: 0.2321172

the BLP diversion ratio from product 1 to 3 is lower than the logit diversion ratio but a bit higher than the nested logit diversion ratio. the BLP diversion ratio from product 3 to 1 is lower than the logit and nested logit diversion ratios.

Let's compare in percentages:

In [13]:
D_13_logit = 0.1921416
D_31_logit = 0.3445739 
D_13_nlogit = 0.09510307
D_31_nlogit =  0.2321172 

print("Differences in % w.r.t to logit: ", (diversions[0,2]-D_13_logit)*100/diversions[0,2], "and" , (diversions[2,0]-D_31_logit)*100/diversions[2,0])
print("Differences in % w.r.t to nested logit: ", (diversions[0,2]-D_13_nlogit)*100/diversions[0,2], "and" , (diversions[2,0]-D_31_nlogit)*100/diversions[2,0])

Differences in % w.r.t to logit:  -84.05220143642265 and -126.35151840013832
Differences in % w.r.t to nested logit:  8.90088665410196 and -52.478410775710486