Replication of the classic BLP (1995) paper

Code is a summer self-exercise as preparation for graduate school:

Some helpful sources:

  1. Chris Conlon and Jeff Gortmaker's pyblp package, with their working paper: https://jeffgortmaker.com/files/pyblp.pdf
  1. Mark Ponder's Blog: https://mark-ponder.com/tutorials/static-discrete-choice-models/random-coefficients-blp/
    • Note that code on that website has quite a few fatal typos in it
  1. Kohei Kawaguchi's online empirical IO assignments: https://kohei-kawaguchi.github.io/EmpiricalIO/assignment4.html
  1. I also try to follow the procedure of Aviv Nevo's Practioners Guide (JEMS 2000)

Note: Code is run on a fairly budget laptop (i5-8250U, 8GB RAM), so performance should be better in faster computers

In [1]:
import pandas as pd
import numpy as np
from scipy.optimize import minimize
import scipy
from numba import jit, njit, prange
import time
import multiprocessing as mp
import pickle
from sklearn.linear_model import LinearRegression
import pyblp

R code for cleaning

Dataframe cleaning is a lot more concise in R...

library(magrittr)

df <- cars %>%
  dplyr::mutate(ln_hpwt = log(hpwt),
                ln_space = log(space),
                ln_mpg = log(mpg),
                ln_mpd = log(mpd),
                ln_price = log(price),
                trend = market + 70,
                cons = 1) %>%
  dplyr::group_by(model_year) %>%
  dplyr::mutate(s_0 = log(1-sum(share))) %>%
  dplyr::ungroup() %>%
  dplyr::mutate(s_i = log(share)) %>%
  dplyr::mutate(dif = s_i - s_0, 
                dif_2 = log(share) - log(share_out),
                ln_price = log(price)) %>%
  dplyr::arrange(market, firmid)
In [2]:
# Read the dataframe
df = pd.read_csv("../data/blp_1995_data.csv")
df = df.drop(df.columns[0], axis = 1)

# python code for cleaning
df[["ln_hpwt", "ln_space", "ln_mpg", "ln_mpd", "ln_price"]] = \
    df[["hpwt", "space", "mpg", "mpd", "price"]].apply(lambda x: np.log(x))

# instrument change
df["trend"] = df.market.map(lambda x: x + 70)  # use with non pyblp instruments
# df["trend"] = df.market

df["cons"] = 1

df["s_0"] = np.log(1 - df.share.groupby(df["model_year"]).transform("sum"))

df["s_i"] = np.log(df.share)
df["dif"] = df.s_i - df.s_0
df["dif_2"] = np.log(df.share) - np.log(df.share_out)
df["ln_price"] = np.log(df.price)

df.head()
Out[2]:
modelvec newmodv model_year id firmid market hpwt space air mpd ... ln_space ln_mpg ln_mpd ln_price trend cons s_0 s_i dif dif_2
0 AMGREM AMGREM71 71 129 15 1 0.528997 1.1502 0 1.888146 ... 0.139936 0.528862 0.635595 1.596515 71 1 -0.171483 -6.858013 -6.686531 -6.730300
1 AMHORN AMHORN71 71 130 15 1 0.494324 1.2780 0 1.935989 ... 0.245296 0.553885 0.660618 1.707662 71 1 -0.171483 -7.308233 -7.136750 -7.180520
2 AMJAVL AMJAVL71 71 132 15 1 0.467613 1.4592 0 1.716799 ... 0.377888 0.433729 0.540462 1.961311 71 1 -0.171483 -7.983628 -7.812146 -7.855915
3 AMMATA AMMATA71 71 134 15 1 0.426540 1.6068 0 1.687871 ... 0.474245 0.416735 0.523468 1.922716 71 1 -0.171483 -7.557843 -7.386360 -7.430130
4 AMAMBS AMAMBS71 71 136 15 1 0.452489 1.6458 0 1.504286 ... 0.498227 0.301585 0.408318 2.189237 71 1 -0.171483 -7.724201 -7.552718 -7.596488

5 rows × 30 columns

In [3]:
# this is here because we may want to use their instruments
product_data = pd.read_csv(pyblp.data.BLP_PRODUCTS_LOCATION)
In [4]:
# demand variables
X = df[["cons", "hpwt", "air", "mpd", "space"]].as_matrix()

# suppy variables
W = df[["cons", "ln_hpwt", "air", "ln_mpg", "ln_space", "trend"]].as_matrix()

# price
p = df.price.values

# initial delta_0 estimate: log(share) - log(share outside good)
delta_0 = df.dif_2.values

# number of goods per market
J = df.groupby("year").sum().cons.values

# number of draws per market
N = 500

# number of markets
T = len(J)

# Estimated log income means for years 1971 - 1990
incomeMeans = [2.01156, 2.06526, 2.07843, 2.05775, 2.02915, 2.05346, 2.06745,
               2.09805, 2.10404, 2.07208, 2.06019, 2.06561, 2.07672, 2.10437, 2.12608, 2.16426,
               2.18071, 2.18856, 2.21250, 2.18377]

# standard deviation of log incomes, assuming empirically given in 1995
sigma_v = 1.72

# number of terms that have the random coefficient
#  according to table 4, this is constant, hp/wt, air, mp$, size, price
k = 5

# markets for itj
markets = df.market.values

# unique markets
marks = np.unique(df.market)

# firms
firms = np.reshape(df.firmid.values, (-1,1))
In [5]:
# define a class so I can repeatedly update the delta value
class delta:
    def __init__(self, delta):
        self.delta = delta

# initialize a delta object using the delta_0 values
d = delta(delta_0)

Even though the orignal BLP (1995) paper did not use different draws in each market, I decided that this is a better idea than using the same draws for all markets.

In [6]:
# set seed
np.random.seed(135567)

m_t = np.repeat(incomeMeans, N)

# matrix of simulated values
#  number of rows = number of simulations
#  treat last column is price/income
# note that BLP treat it as the same individuals in the market across all years (Nevo 2000)

# different draws for each market
V = np.reshape(np.random.standard_normal((k + 1) * N * T), (T * N, k + 1))

# # same draws for each market
# V = pd.read_csv("../v_ik.csv", header = None).as_matrix()
# V = np.reshape(np.random.standard_normal((k + 1) * N), (N, k + 1))

# income if we have different draws per market
y_it = np.exp(m_t + sigma_v * V[:, k]).reshape(T,N).T

# draws for income if same draws in every market
# y_it = np.exp(m_t + sigma_v * np.tile(V[:, k], len(incomeMeans))).reshape(T,N).T

If same draws are used in each market, you must remove the "N * t" indexing and just have "i"

In [7]:
# the loops that calculate utility in a separate function so that it can be
#  run in parallel. 
@jit(nopython = True, parallel = True)
def util_iter(out, x, v, p, y, delta, theta_2, J, T, N):
    # first iterate over the individuals 
    for i in prange(N):
        # iterator through t and j
        tj = 0
        
        # iterate over the markets
        for t in prange(T):
            # market size of market t
            mktSize = J[t]
            
            # income for individual i in market t
            y_im = y[i, t]
            
            # iterate over goods in a particular market
            for j in prange(mktSize):
                out[tj, i] = delta[tj] + \
                v[N * t + i, 0] * theta_2[0] * x[tj, 0] + \
                v[N * t + i, 1] * theta_2[1] * x[tj, 1] + \
                v[N * t + i, 2] * theta_2[2] * x[tj, 2] + \
                v[N * t + i, 3] * theta_2[3] * x[tj, 3] + \
                v[N * t + i, 4] * theta_2[4] * x[tj, 4] - \
                theta_2[5] / y_im * p[tj]
                
                tj += 1
    return out

# computes indirect utility given parameters
#  x: matrix of demand characteristics
#  v: monte carlo draws of N simulations
#  p: price vector
#  y: income of individuals
#  delta: guess for the mean utility
#  theta_2: non-linear params (sigma - can think of as stdev's)
#  J: vector of number of goods per market
#  T: numer of markets
#  N: number of simulations

@jit
def compute_indirect_utility(x, v, p, y, delta, theta_2, J, T, N):
    # make sure theta_2 are positive
    theta_2 = np.abs(theta_2)
    
    # output matrix
    out = np.zeros((sum(J), N))
    
    # call the iteration function to calculate utilities
    out = util_iter(out, x, v, p, y, delta, theta_2, J, T, N)
     
    return out
In [8]:
# computes the implied shares of goods in each market given inputs
# same inputs as above function

@jit
def compute_share(x, v, p, y, delta, theta_2, J, T, N):
    q = np.zeros((np.sum(J), N))
    
    # obtain vector of indirect utilities
    u = compute_indirect_utility(x, v, p, y, delta, theta_2, J, T, N)
    
    # exponentiate the utilities
    exp_u = np.exp(u)
    
    # pointer to first good in the market
    first_good = 0
            
    for t in range(T):
        # market size of market t
        mktSize = J[t]

        # calculate the numerator of the share eq
        numer = exp_u[first_good:first_good + mktSize,:]

        # calculate the denom of the share eq
        denom = 1 + numer.sum(axis = 0)    
          
        # calculate the quantity each indv purchases of each good in each market
        q[first_good:first_good + mktSize,:] = numer/denom
        
        first_good += mktSize
    
    # to obtain shares, assume that each simulation carries the same weight.
    # this averages each row, which is essentially computing the sahres for each
    #  good in each market. 
    s = np.matmul(q, np.repeat(1/N, N))
    
    return [q,s]
In [9]:
@jit
def solve_delta(s, x, v, p, y, delta, theta_2, J, T, N, tol):
    # define the tolerance variable
    eps = 10
    
    # renaming delta as delta^r
    delta_old = delta
    
    while eps > tol:
        # Aviv's step 1: obtain predicted shares and quantities
        q_s = compute_share(x, v, p, y, delta_old, 
                            theta_2, J, T, N)
        
        # extract the shares
        sigma_jt = q_s[1]
        
        # step 2: use contraction mapping to find delta
        delta_new = delta_old + np.log(s/sigma_jt)
        
        # update tolerance
        eps = np.max(np.abs(delta_new - delta_old))
        
        delta_old = delta_new.copy()
    
    return delta_old
        
In [10]:
# calculates the marginal costs given probabilities and shares
#  q_s: output of compute share, a list of probabilities matrix(q) and shares vector(s)
#  firms: vector of firms operating in each market (length is JxT)
#  marks: vector of unique markets (length T)
#  markets: vector indicating observation in which market (length JxT)

@jit
def calc_mc(q_s, firms, p, y, alpha, J, T, N, marks, markets):
    
    # declare output vector
    out = np.zeros((np.sum(J)))
    
    # make sure the value of alpha is positive
    alpha = np.abs(alpha)
    
    # read in quantities
    q = q_s[0]
    
    # read in shares
    s = q_s[1].reshape((-1,1))
    
    # reshape some vectors into column vectors
    p = p.reshape((-1,1))
    
    # iterate over markets
    for m in marks:
        # obtain list of firms operating in that market/year
        firm_yr = firms[markets == m]
        
        # obtain list of prices of goods in that market/year
        price = p[markets == m]
        
        # J_t x J_t block matrix of 1's indicating goods belonging to same firm
        #  in that market/year
        # Also known as the ownership matrix
        same_firm = np.equal(firm_yr, np.transpose(firm_yr))
        
        # obtain matrix of probabilities for all simulations for goods in that 
        #  market/year
        yr = q[markets == m,:]
        
        # obtain number of goods in that market
        nobs = np.size(yr, 0)
        
        # this is the omega matrix initializing        
        grad = np.zeros((nobs, nobs))
        
        # calculate the omega matix by iterating through all individuals
        #  Omega matrix is cross-price deriv element-wise product with
        #  ownership matrix
        for i in range(N):
            yr_i = yr[:, i].reshape((-1, 1))
            grad += alpha / y[i, m - 1] * same_firm * 1/N * \
            (yr_i @ yr_i.T - np.diag(yr[:,i]))
        
        # Omega matrix actually requires negative cross price derivs
        subMatrix = -grad
        
        # now obtain the marginal costs
        b = np.linalg.inv(subMatrix) @ s[markets == m]
        mc = price - b
        mc[mc < 0] = .001
        
        # update entries in the output vector
        out[markets == m] = mc.flatten()
        
    return out
In [11]:
# generate instruments for BLP paper
# To be honest: I have no idea what's going on here.
# Additionally, there is discussion that the instruments BLP created were wrong.
#  However, this code below generates the instruments that correctly matches
#  the instruments used in the BLP paper.
# One can use packages to generate more "optimal" instruments instead, like Hausmann 
#  instruments
# A source of instruments exists in the Conlon and Gortmaker's pyblp package
def gen_inst(inx, firms, marks, markets):
    totMarket = np.zeros((np.size(inx, axis = 0), np.size(inx, axis = 1)))
    totFirm = np.zeros((np.size(inx, axis = 0), np.size(inx, axis = 1)))
    
    for m in marks:
        sub = inx[markets == m, :]
        firminfo = firms[markets == m]
        same_firm = np.equal(firminfo, np.transpose(firminfo))
        
        z_1 = np.zeros((np.size(sub,axis = 0), np.size(sub, axis = 1)))
        
        for i in range(np.size(sub, axis = 1)):
            z_1[:,i] = (sub[:,i].reshape((-1,1)) * same_firm).T.sum(axis = 0)
        
        totFirm[markets == m,:] = z_1
        
        z_1 = np.zeros((np.size(sub,axis = 0), np.size(sub, axis = 1)))
        
        for i in range(np.size(sub, axis = 1)):
            z_1[:,i] = (sub[:,i].reshape((-1,1)) * (same_firm + np.logical_not(same_firm))).sum(axis = 0)
            
        totMarket[markets == m, :] = z_1
        
    return [totFirm, totMarket]

            
In [12]:
# If you're looking at the four steps Aviv lists in his appendix, start here
# This is the objective function that we optimize the non-linear parameters over
def objective(theta_2, s, X, V, p, y, J, T, N, marks, markets, tol, 
              Z, Z_s, W, weigh, firms):
    
    obs = np.sum(J)
    
    # Aviv's step 1 & 2:
    d.delta = solve_delta(s, X, V, p, y, d.delta, theta_2, J, T, N, tol)
    
    # obtain the actual implied quantities and shares from converged delta
    q_s = compute_share(X, V, p, y, d.delta, theta_2, J, T, N)
    
    # calculate marginal costs
    mc = calc_mc(q_s, firms, p, y, theta_2[5], J, T, N, marks, markets).reshape((-1,1))
    
    # since we are using both demand and supply side variables,
    #  we want to stack the estimated delta and estimated mc
    y2 = np.vstack((d.delta.reshape((-1,1)), np.log(mc)))
    
    # create characteristics matrix that includes both supply and demand side
    #  with demand characteristics on the bottom left and supply on the top right
    x = scipy.linalg.block_diag(X,W)
    
    # create matrix of supply and demand instruments, again with
    #  demand instruments on the right and supply on the left (top/down changed)    
    z = scipy.linalg.block_diag(Z,Z_s)
    
    # get linear parameters (this FOC is from Aviv's appendix)
    b = np.linalg.inv(x.T @ z @ weigh @ z.T @ x) @ (x.T @ z @ weigh @ z.T @ y2)
    
    # Step 3: get the error term xi (also called omega)
    xi_w = y2 - x @ b
    
    # computeo g_bar in GMM methods
    g = z.T @ xi_w / obs
    
    obj = float(obs**2 * g.T @ weigh @ g)
   
    print([theta_2, obj])
    
    return obj

Table 1

In [13]:
# create table
table_1 = pd.DataFrame()

# calculate weighted means of the following column with quantity weights
table_1[[
    "P", "Dom", "Japan", "Eu", "HP_Wt", "Size", "Air", "MPG", "MPD", "drop"
   ]] = df[["price", "domestic", "japan", 
    "european", "hpwt", "space", 
    "air", "mpg", "mpd", "quantity"]] \
    .groupby(df.year) \
    .apply(lambda x: pd.Series(np.average(x, weights=x["quantity"], axis = 0)))

# count number of models per year/market
table_1.insert(0, "num_models", df.groupby("year").sum().cons.values)

# mean quantity per year/market
table_1.insert(1, "Q", df.quantity.groupby(df.year).mean()) 

# delete the extraneous weighted quantity column
table_1.drop('drop', axis=1, inplace=True)
In [14]:
# round the table to 3 digits to emulate the paper
np.round(table_1, 3)
Out[14]:
num_models Q P Dom Japan Eu HP_Wt Size Air MPG MPD
year
71 92 86.892 7.868 0.866 0.057 0.077 0.490 1.496 0.000 1.662 1.849
72 89 98.623 7.979 0.892 0.042 0.066 0.391 1.510 0.014 1.619 1.875
73 86 92.785 7.535 0.932 0.040 0.028 0.364 1.529 0.022 1.589 1.818
74 72 105.119 7.506 0.887 0.050 0.064 0.347 1.510 0.026 1.567 1.452
75 93 84.775 7.821 0.853 0.083 0.064 0.337 1.479 0.054 1.584 1.503
76 99 93.382 7.787 0.876 0.081 0.043 0.338 1.508 0.059 1.759 1.696
77 95 97.727 7.651 0.837 0.112 0.051 0.340 1.467 0.032 1.947 1.835
78 95 99.444 7.645 0.855 0.107 0.039 0.346 1.405 0.034 1.982 1.929
79 102 82.742 7.599 0.803 0.158 0.038 0.348 1.343 0.047 2.061 1.657
80 103 71.567 7.718 0.773 0.191 0.036 0.350 1.296 0.078 2.215 1.466
81 116 62.030 8.349 0.741 0.213 0.046 0.349 1.286 0.094 2.363 1.559
82 110 61.893 8.831 0.714 0.235 0.051 0.347 1.277 0.134 2.440 1.817
83 115 67.878 8.821 0.734 0.215 0.051 0.351 1.276 0.126 2.601 2.087
84 113 85.933 8.870 0.783 0.179 0.038 0.361 1.293 0.129 2.469 2.117
85 136 78.143 8.938 0.761 0.191 0.048 0.372 1.265 0.140 2.261 2.024
86 130 83.756 9.382 0.733 0.216 0.050 0.379 1.249 0.176 2.416 2.856
87 143 67.667 9.965 0.702 0.245 0.052 0.395 1.246 0.229 2.327 2.789
88 150 67.078 10.069 0.717 0.237 0.045 0.396 1.251 0.237 2.334 2.919
89 147 62.914 10.321 0.690 0.261 0.049 0.406 1.259 0.289 2.310 2.806
90 131 66.377 10.337 0.682 0.276 0.043 0.419 1.270 0.308 2.270 2.852

Table 3

In [15]:
# table 3 column 1
t3c1 = LinearRegression().fit(df[["hpwt", "air", "mpd", "space", "price"]], df.dif_2)
np.hstack((t3c1.intercept_, t3c1.coef_))
Out[15]:
array([-10.07300751,  -0.12309488,  -0.03441478,   0.26546568,
         2.34191416,  -0.08860627])
In [16]:
# table 3 column 3
t3c3 = LinearRegression().fit(
    df[["ln_hpwt", "air", "ln_mpg", "ln_space", "trend"]],
    df.ln_price)
np.hstack((t3c3.intercept_, t3c3.coef_))
Out[16]:
array([ 1.88192125,  0.52033668,  0.6797513 , -0.47064027,  0.12482708,
        0.01283075])
In [17]:
# ====================================================================
# IV Reg in table 3 column 2
# ====================================================================

# close to BLP original instruments
tempDemand = gen_inst(X, firms, marks, markets)
tempSupply = gen_inst(W, firms, marks, markets)

Z = np.hstack((X, tempDemand[0], tempDemand[1]))
baseData = np.hstack((p.reshape((-1,1)), X))

Z_s = np.hstack((W, tempSupply[0], tempSupply[1], df.mpd.values.reshape((-1,1))))

# # pyblp instruments
# Z = np.hstack((product_data[["demand_instruments0", "demand_instruments1", 
#                              "demand_instruments2", "demand_instruments3",
#                             "demand_instruments4", "demand_instruments5"]], X))

# Z_s = np.hstack((product_data[["supply_instruments0", "supply_instruments1", 
#                              "supply_instruments2", "supply_instruments3",
#                             "supply_instruments4", "supply_instruments5"]], W))

baseData = np.hstack((X, p.reshape((-1,1))))

zxw1 = Z.T @ baseData

bx1 = np.linalg.inv(zxw1.T @ zxw1)@ zxw1.T @ Z.T @ delta_0

e = delta_0 - baseData @ bx1

g_ind = e.reshape((-1,1)) * Z

demean = g_ind - g_ind.mean(axis=0).reshape((1,-1))

vg = demean.T @ demean / demean.shape[0]

w0 = np.linalg.inv(vg)

t3c2 = np.linalg.inv(zxw1.T @ w0 @ zxw1) @ zxw1.T @ w0 @ Z.T @ delta_0
t3c2
Out[17]:
array([-9.27629294,  1.94935115,  1.28739148,  0.05456147,  2.35760254,
       -0.21578695])
In [18]:
# obtain block-diag matrix of supply and demand instruments
z = scipy.linalg.block_diag(Z,Z_s)

# Recommended initial weighting matrix from Aviv's appendix
w1 = np.linalg.inv(z.T @ z)
In [19]:
# # run objective function once to see if it is working
# initial parameter guess (from BLP(1995))
theta_2 = [3.612, 4.628, 1.818, 1.050, 2.056, 43.501]

t0 = time.time()

obj = objective(theta_2, 
                df.share.values, 
                X, V, p, y_it, J, T, N, marks, markets, 1e-4, 
              Z, Z_s, W, w1, firms)


time.time() - t0
[[3.612, 4.628, 1.818, 1.05, 2.056, 43.501], 1264.5134198310147]
Out[19]:
12.652115106582642
In [20]:
# constraints for the minimization problem
#  not needed because utility considers abosolute value of params before calculating

# Step 4: search over parameters that minimize the objective function
t1 = time.time()

# set bounds for optimization (not entirely needed but oh well)
bnds = ((0,np.inf), (0,np.inf), (0,np.inf), 
        (0,np.inf), (0,np.inf), (5,np.inf))

res_init_wt = minimize(objective,
                      theta_2, 
                      args = (df.share.values, X, V, p, y_it, 
                              J, T, N, marks, markets, 1e-4, 
                              Z, Z_s, W, w1, firms), 
                      bounds = bnds,
                      method = "L-BFGS-B",
                      options = {'maxiter': 1000, 'maxfun': 1000, 'eps': 1e-3},
                      tol = 1e-4)
    

    
time.time() - t1
[array([ 3.612,  4.628,  1.818,  1.05 ,  2.056, 43.501]), 1264.5154952422406]
[array([ 3.613,  4.628,  1.818,  1.05 ,  2.056, 43.501]), 1264.5495867005272]
[array([ 3.612,  4.629,  1.818,  1.05 ,  2.056, 43.501]), 1264.548759632199]
[array([ 3.612,  4.628,  1.819,  1.05 ,  2.056, 43.501]), 1264.673516194809]
[array([ 3.612,  4.628,  1.818,  1.051,  2.056, 43.501]), 1265.684936696364]
[array([ 3.612,  4.628,  1.818,  1.05 ,  2.057, 43.501]), 1264.5471762736206]
[array([ 3.612,  4.628,  1.818,  1.05 ,  2.056, 43.502]), 1264.577495000419]
[array([ 3.51951909,  4.50950564,  1.7714523 ,  1.02311602,  2.0033586 ,
       42.51522832]), 1198.2776516959893]
[array([ 3.52051909,  4.50950564,  1.7714523 ,  1.02311602,  2.0033586 ,
       42.51522832]), 1198.2844462916214]
[array([ 3.51951909,  4.51050564,  1.7714523 ,  1.02311602,  2.0033586 ,
       42.51522832]), 1198.26665249537]
[array([ 3.51951909,  4.50950564,  1.7724523 ,  1.02311602,  2.0033586 ,
       42.51522832]), 1198.3715560402566]
[array([ 3.51951909,  4.50950564,  1.7714523 ,  1.02411602,  2.0033586 ,
       42.51522832]), 1199.3384615923162]
[array([ 3.51951909,  4.50950564,  1.7714523 ,  1.02311602,  2.0043586 ,
       42.51522832]), 1198.2405888372955]
[array([ 3.51951909,  4.50950564,  1.7714523 ,  1.02311602,  2.0033586 ,
       42.51622832]), 1198.2695914920744]
[array([ 3.62822966,  4.71009164,  1.55533262,  0.        ,  2.19643851,
       43.89476902]), 613.1277572623852]
[array([ 3.62922966,  4.71009164,  1.55533262,  0.        ,  2.19643851,
       43.89476902]), 613.1696522584633]
[array([ 3.62822966,  4.71109164,  1.55533262,  0.        ,  2.19643851,
       43.89476902]), 613.1740015556185]
[array([ 3.62822966,  4.71009164,  1.55633262,  0.        ,  2.19643851,
       43.89476902]), 613.2594959848929]
[array([3.62822966e+00, 4.71009164e+00, 1.55533262e+00, 1.00000000e-03,
       2.19643851e+00, 4.38947690e+01]), 613.3784168537053]
[array([ 3.62822966,  4.71009164,  1.55533262,  0.        ,  2.19743851,
       43.89476902]), 613.326291759063]
[array([ 3.62822966,  4.71009164,  1.55533262,  0.        ,  2.19643851,
       43.89576902]), 613.3045726820188]
[array([ 3.49296059,  4.57478362,  1.46183984,  0.        ,  2.08351271,
       42.22087485]), 570.8216310465589]
[array([ 3.49396059,  4.57478362,  1.46183984,  0.        ,  2.08351271,
       42.22087485]), 570.8175196680352]
[array([ 3.49296059,  4.57578362,  1.46183984,  0.        ,  2.08351271,
       42.22087485]), 570.7930180473637]
[array([ 3.49296059,  4.57478362,  1.46283984,  0.        ,  2.08351271,
       42.22087485]), 570.8467897691277]
[array([3.49296059e+00, 4.57478362e+00, 1.46183984e+00, 1.00000000e-03,
       2.08351271e+00, 4.22208749e+01]), 570.9210219919003]
[array([ 3.49296059,  4.57478362,  1.46183984,  0.        ,  2.08451271,
       42.22087485]), 570.8711849027618]
[array([ 3.49296059,  4.57478362,  1.46183984,  0.        ,  2.08351271,
       42.22187485]), 570.8508699696847]
[array([ 3.51067955,  4.6925611 ,  1.41868331,  0.        ,  1.98295798,
       41.93278984]), 561.5728301851626]
[array([ 3.51167955,  4.6925611 ,  1.41868331,  0.        ,  1.98295798,
       41.93278984]), 561.5682181194901]
[array([ 3.51067955,  4.6935611 ,  1.41868331,  0.        ,  1.98295798,
       41.93278984]), 561.543091556416]
[array([ 3.51067955,  4.6925611 ,  1.41968331,  0.        ,  1.98295798,
       41.93278984]), 561.5968075886202]
[array([3.51067955e+00, 4.69256110e+00, 1.41868331e+00, 1.00000000e-03,
       1.98295798e+00, 4.19327898e+01]), 561.6459706369062]
[array([ 3.51067955,  4.6925611 ,  1.41868331,  0.        ,  1.98395798,
       41.93278984]), 561.6077261728315]
[array([ 3.51067955,  4.6925611 ,  1.41868331,  0.        ,  1.98295798,
       41.93378984]), 561.6055931838306]
[array([ 3.58155541,  5.16367102,  1.24605719,  0.        ,  1.58073907,
       40.78044978]), 523.0548053801655]
[array([ 3.58255541,  5.16367102,  1.24605719,  0.        ,  1.58073907,
       40.78044978]), 523.0626028463197]
[array([ 3.58155541,  5.16467102,  1.24605719,  0.        ,  1.58073907,
       40.78044978]), 523.0418920350324]
[array([ 3.58155541,  5.16367102,  1.24705719,  0.        ,  1.58073907,
       40.78044978]), 523.0974054969456]
[array([3.58155541e+00, 5.16367102e+00, 1.24605719e+00, 1.00000000e-03,
       1.58073907e+00, 4.07804498e+01]), 523.1339700969833]
[array([ 3.58155541,  5.16367102,  1.24605719,  0.        ,  1.58173907,
       40.78044978]), 523.1056708840731]
[array([ 3.58155541,  5.16367102,  1.24605719,  0.        ,  1.58073907,
       40.78144978]), 523.1033430769804]
[array([ 3.8601008 ,  7.01515496,  0.5676285 ,  0.        ,  0.        ,
       36.25169969]), 398.12078023237416]
[array([ 3.8611008 ,  7.01515496,  0.5676285 ,  0.        ,  0.        ,
       36.25169969]), 398.1306707757394]
[array([ 3.8601008 ,  7.01615496,  0.5676285 ,  0.        ,  0.        ,
       36.25169969]), 398.13124564649866]
[array([ 3.8601008 ,  7.01515496,  0.5686285 ,  0.        ,  0.        ,
       36.25169969]), 398.1530784471046]
[array([3.86010080e+00, 7.01515496e+00, 5.67628503e-01, 1.00000000e-03,
       0.00000000e+00, 3.62516997e+01]), 398.1955752093014]
[array([3.86010080e+00, 7.01515496e+00, 5.67628503e-01, 0.00000000e+00,
       1.00000000e-03, 3.62516997e+01]), 398.17325765552835]
[array([ 3.8601008 ,  7.01515496,  0.5676285 ,  0.        ,  0.        ,
       36.25269969]), 398.1621786699984]
[array([ 3.94318709,  8.01505285,  0.089834  ,  0.        ,  0.        ,
       32.73127883]), 374.8132941580275]
[array([ 3.94418709,  8.01505285,  0.089834  ,  0.        ,  0.        ,
       32.73127883]), 374.8180228781935]
[array([ 3.94318709,  8.01605285,  0.089834  ,  0.        ,  0.        ,
       32.73127883]), 374.8201592951552]
[array([ 3.94318709,  8.01505285,  0.090834  ,  0.        ,  0.        ,
       32.73127883]), 374.8188251333163]
[array([3.94318709e+00, 8.01505285e+00, 8.98339984e-02, 1.00000000e-03,
       0.00000000e+00, 3.27312788e+01]), 374.872910238051]
[array([3.94318709e+00, 8.01505285e+00, 8.98339984e-02, 0.00000000e+00,
       1.00000000e-03, 3.27312788e+01]), 374.84543212183735]
[array([ 3.94318709,  8.01505285,  0.089834  ,  0.        ,  0.        ,
       32.73227883]), 374.8381715767841]
[array([ 4.12762049,  9.43679208,  0.        ,  0.        ,  0.        ,
       27.16844373]), 356.81003052758865]
[array([ 4.12862049,  9.43679208,  0.        ,  0.        ,  0.        ,
       27.16844373]), 356.8124400246089]
[array([ 4.12762049,  9.43779208,  0.        ,  0.        ,  0.        ,
       27.16844373]), 356.82018739584834]
[array([4.12762049e+00, 9.43679208e+00, 1.00000000e-03, 0.00000000e+00,
       0.00000000e+00, 2.71684437e+01]), 356.84048631372895]
[array([4.12762049e+00, 9.43679208e+00, 0.00000000e+00, 1.00000000e-03,
       0.00000000e+00, 2.71684437e+01]), 356.8152952905704]
[array([4.12762049e+00, 9.43679208e+00, 0.00000000e+00, 0.00000000e+00,
       1.00000000e-03, 2.71684437e+01]), 356.82745519704815]
[array([ 4.12762049,  9.43679208,  0.        ,  0.        ,  0.        ,
       27.16944373]), 356.82722997427794]
[array([ 4.47368382, 11.63538738,  0.        ,  0.        ,  0.        ,
       18.09171868]), 384.1842210302302]
[array([ 4.47468382, 11.63538738,  0.        ,  0.        ,  0.        ,
       18.09171868]), 384.1981207298126]
[array([ 4.47368382, 11.63638738,  0.        ,  0.        ,  0.        ,
       18.09171868]), 384.20040004553914]
[array([4.47368382e+00, 1.16353874e+01, 1.00000000e-03, 0.00000000e+00,
       0.00000000e+00, 1.80917187e+01]), 384.168733860721]
[array([4.47368382e+00, 1.16353874e+01, 0.00000000e+00, 1.00000000e-03,
       0.00000000e+00, 1.80917187e+01]), 384.24227351121846]
[array([4.47368382e+00, 1.16353874e+01, 0.00000000e+00, 0.00000000e+00,
       1.00000000e-03, 1.80917187e+01]), 384.1985737325732]
[array([ 4.47368382, 11.63538738,  0.        ,  0.        ,  0.        ,
       18.09271868]), 384.1890790905862]
[array([ 4.20842446,  9.95015255,  0.        ,  0.        ,  0.        ,
       25.04907605]), 356.5449296848203]
[array([ 4.20942446,  9.95015255,  0.        ,  0.        ,  0.        ,
       25.04907605]), 356.55463626779976]
[array([ 4.20842446,  9.95115255,  0.        ,  0.        ,  0.        ,
       25.04907605]), 356.5807170657544]
[array([4.20842446e+00, 9.95015255e+00, 1.00000000e-03, 0.00000000e+00,
       0.00000000e+00, 2.50490761e+01]), 356.52041475063965]
[array([4.20842446e+00, 9.95015255e+00, 0.00000000e+00, 1.00000000e-03,
       0.00000000e+00, 2.50490761e+01]), 356.67007618761346]
[array([4.20842446e+00, 9.95015255e+00, 0.00000000e+00, 0.00000000e+00,
       1.00000000e-03, 2.50490761e+01]), 356.577233426469]
[array([ 4.20842446,  9.95015255,  0.        ,  0.        ,  0.        ,
       25.05007605]), 356.5723169121138]
[array([ 4.26958666, 10.3387257 ,  0.        ,  0.        ,  0.        ,
       23.44488284]), 359.3874955527839]
[array([ 4.27058666, 10.3387257 ,  0.        ,  0.        ,  0.        ,
       23.44488284]), 359.3938108195522]
[array([ 4.26958666, 10.3397257 ,  0.        ,  0.        ,  0.        ,
       23.44488284]), 359.39654309323436]
[array([4.26958666e+00, 1.03387257e+01, 1.00000000e-03, 0.00000000e+00,
       0.00000000e+00, 2.34448828e+01]), 359.3954104699332]
[array([4.26958666e+00, 1.03387257e+01, 0.00000000e+00, 1.00000000e-03,
       0.00000000e+00, 2.34448828e+01]), 359.4238664572987]
[array([4.26958666e+00, 1.03387257e+01, 0.00000000e+00, 0.00000000e+00,
       1.00000000e-03, 2.34448828e+01]), 359.40190059108886]
[array([ 4.26958666, 10.3387257 ,  0.        ,  0.        ,  0.        ,
       23.44588284]), 359.39767669111734]
[array([ 4.22266172, 10.04060415,  0.        ,  0.        ,  0.        ,
       24.67565384]), 358.33224782288966]
[array([ 4.22366172, 10.04060415,  0.        ,  0.        ,  0.        ,
       24.67565384]), 358.3406867774916]
[array([ 4.22266172, 10.04160415,  0.        ,  0.        ,  0.        ,
       24.67565384]), 358.3487410543649]
[array([4.22266172e+00, 1.00406042e+01, 1.00000000e-03, 0.00000000e+00,
       0.00000000e+00, 2.46756538e+01]), 358.35177825448204]
[array([4.22266172e+00, 1.00406042e+01, 0.00000000e+00, 1.00000000e-03,
       0.00000000e+00, 2.46756538e+01]), 358.3790822220071]
[array([4.22266172e+00, 1.00406042e+01, 0.00000000e+00, 0.00000000e+00,
       1.00000000e-03, 2.46756538e+01]), 358.3591721465296]
[array([ 4.22266172, 10.04060415,  0.        ,  0.        ,  0.        ,
       24.67665384]), 358.35409692507386]
[array([ 4.21068588,  9.96451971,  0.        ,  0.        ,  0.        ,
       24.98976239]), 358.0085286521671]
[array([ 4.21168588,  9.96451971,  0.        ,  0.        ,  0.        ,
       24.98976239]), 358.0151200241609]
[array([ 4.21068588,  9.96551971,  0.        ,  0.        ,  0.        ,
       24.98976239]), 358.02225107989517]
[array([4.21068588e+00, 9.96451971e+00, 1.00000000e-03, 0.00000000e+00,
       0.00000000e+00, 2.49897624e+01]), 358.0258270606038]
[array([4.21068588e+00, 9.96451971e+00, 0.00000000e+00, 1.00000000e-03,
       0.00000000e+00, 2.49897624e+01]), 358.05076188249205]
[array([4.21068588e+00, 9.96451971e+00, 0.00000000e+00, 0.00000000e+00,
       1.00000000e-03, 2.49897624e+01]), 358.0321134231381]
[array([ 4.21068588,  9.96451971,  0.        ,  0.        ,  0.        ,
       24.99076239]), 358.02718058165624]
[array([ 4.20842446,  9.95015255,  0.        ,  0.        ,  0.        ,
       25.04907605]), 356.5603186619384]
[array([ 4.20942446,  9.95015255,  0.        ,  0.        ,  0.        ,
       25.04907605]), 356.56456965421035]
[array([ 4.20842446,  9.95115255,  0.        ,  0.        ,  0.        ,
       25.04907605]), 356.58713968846683]
[array([4.20842446e+00, 9.95015255e+00, 1.00000000e-03, 0.00000000e+00,
       0.00000000e+00, 2.50490761e+01]), 356.52449369055813]
[array([4.20842446e+00, 9.95015255e+00, 0.00000000e+00, 1.00000000e-03,
       0.00000000e+00, 2.50490761e+01]), 356.67274294106176]
[array([4.20842446e+00, 9.95015255e+00, 0.00000000e+00, 0.00000000e+00,
       1.00000000e-03, 2.50490761e+01]), 356.57891732448246]
[array([ 4.20842446,  9.95015255,  0.        ,  0.        ,  0.        ,
       25.05007605]), 356.57356443235966]
[array([ 4.43109527, 11.11125083,  0.42598884,  0.        ,  0.        ,
       19.37425735]), 389.3016958581568]
[array([ 4.43209527, 11.11125083,  0.42598884,  0.        ,  0.        ,
       19.37425735]), 389.3141523565936]
[array([ 4.43109527, 11.11225083,  0.42598884,  0.        ,  0.        ,
       19.37425735]), 389.3156427312736]
[array([ 4.43109527, 11.11125083,  0.42698884,  0.        ,  0.        ,
       19.37425735]), 389.31795717995874]
[array([4.43109527e+00, 1.11112508e+01, 4.25988837e-01, 1.00000000e-03,
       0.00000000e+00, 1.93742573e+01]), 389.34803424382505]
[array([4.43109527e+00, 1.11112508e+01, 4.25988837e-01, 0.00000000e+00,
       1.00000000e-03, 1.93742573e+01]), 389.31204503683216]
[array([ 4.43109527, 11.11125083,  0.42598884,  0.        ,  0.        ,
       19.37525735]), 389.3031402895383]
[array([ 4.24666115, 10.14953457,  0.07315015,  0.        ,  0.        ,
       24.07460486]), 359.40938242275723]
[array([ 4.24766115, 10.14953457,  0.07315015,  0.        ,  0.        ,
       24.07460486]), 359.42000347343287]
[array([ 4.24666115, 10.15053457,  0.07315015,  0.        ,  0.        ,
       24.07460486]), 359.4307119662535]
[array([ 4.24666115, 10.14953457,  0.07415015,  0.        ,  0.        ,
       24.07460486]), 359.4445404897392]
[array([4.24666115e+00, 1.01495346e+01, 7.31501521e-02, 1.00000000e-03,
       0.00000000e+00, 2.40746049e+01]), 359.44648499683444]
[array([4.24666115e+00, 1.01495346e+01, 7.31501521e-02, 0.00000000e+00,
       1.00000000e-03, 2.40746049e+01]), 359.43957695393533]
[array([ 4.24666115, 10.14953457,  0.07315015,  0.        ,  0.        ,
       24.07560486]), 359.4372501352327]
[array([4.21296773e+00, 9.97384308e+00, 8.69168528e-03, 0.00000000e+00,
       0.00000000e+00, 2.49332896e+01]), 358.573462549793]
[array([4.21396773e+00, 9.97384308e+00, 8.69168528e-03, 0.00000000e+00,
       0.00000000e+00, 2.49332896e+01]), 358.54330127811295]
[array([4.21296773e+00, 9.97484308e+00, 8.69168528e-03, 0.00000000e+00,
       0.00000000e+00, 2.49332896e+01]), 358.1951268010661]
[array([4.21296773e+00, 9.97384308e+00, 9.69168528e-03, 0.00000000e+00,
       0.00000000e+00, 2.49332896e+01]), 357.5875991079219]
[array([4.21296773e+00, 9.97384308e+00, 8.69168528e-03, 1.00000000e-03,
       0.00000000e+00, 2.49332896e+01]), 358.2228949674263]
[array([4.21296773e+00, 9.97384308e+00, 8.69168528e-03, 0.00000000e+00,
       1.00000000e-03, 2.49332896e+01]), 358.5577891920151]
[array([4.21296773e+00, 9.97384308e+00, 8.69168528e-03, 0.00000000e+00,
       0.00000000e+00, 2.49342896e+01]), 358.5207632612908]
[array([4.20856524e+00, 9.95088665e+00, 2.69331646e-04, 0.00000000e+00,
       0.00000000e+00, 2.50454881e+01]), 356.5780465764649]
[array([4.20956524e+00, 9.95088665e+00, 2.69331646e-04, 0.00000000e+00,
       0.00000000e+00, 2.50454881e+01]), 356.5827458157032]
[array([4.20856524e+00, 9.95188665e+00, 2.69331646e-04, 0.00000000e+00,
       0.00000000e+00, 2.50454881e+01]), 356.60619595490016]
[array([4.20856524e+00, 9.95088665e+00, 1.26933165e-03, 0.00000000e+00,
       0.00000000e+00, 2.50454881e+01]), 356.5417362889475]
[array([4.20856524e+00, 9.95088665e+00, 2.69331646e-04, 1.00000000e-03,
       0.00000000e+00, 2.50454881e+01]), 356.69431565905916]
[array([4.20856524e+00, 9.95088665e+00, 2.69331646e-04, 0.00000000e+00,
       1.00000000e-03, 2.50454881e+01]), 356.5979307213764]
[array([4.20856524e+00, 9.95088665e+00, 2.69331646e-04, 0.00000000e+00,
       0.00000000e+00, 2.50464881e+01]), 356.5926205628889]
[array([4.20844220e+00, 9.95024508e+00, 3.39482283e-05, 0.00000000e+00,
       0.00000000e+00, 2.50486238e+01]), 356.5722111062827]
[array([4.20944220e+00, 9.95024508e+00, 3.39482283e-05, 0.00000000e+00,
       0.00000000e+00, 2.50486238e+01]), 356.5725928589345]
[array([4.20844220e+00, 9.95124508e+00, 3.39482283e-05, 0.00000000e+00,
       0.00000000e+00, 2.50486238e+01]), 356.5930590511405]
[array([4.20844220e+00, 9.95024508e+00, 1.03394823e-03, 0.00000000e+00,
       0.00000000e+00, 2.50486238e+01]), 356.5288047238047]
[array([4.20844220e+00, 9.95024508e+00, 3.39482283e-05, 1.00000000e-03,
       0.00000000e+00, 2.50486238e+01]), 356.6768595128356]
[array([4.20844220e+00, 9.95024508e+00, 3.39482283e-05, 0.00000000e+00,
       1.00000000e-03, 2.50486238e+01]), 356.5821891240662]
[array([4.20844220e+00, 9.95024508e+00, 3.39482283e-05, 0.00000000e+00,
       0.00000000e+00, 2.50496238e+01]), 356.5766176624098]
[array([4.20842539e+00, 9.95015742e+00, 1.78791688e-06, 0.00000000e+00,
       0.00000000e+00, 2.50490522e+01]), 356.5737211603468]
[array([4.20942539e+00, 9.95015742e+00, 1.78791688e-06, 0.00000000e+00,
       0.00000000e+00, 2.50490522e+01]), 356.57233668666726]
[array([4.20842539e+00, 9.95115742e+00, 1.78791688e-06, 0.00000000e+00,
       0.00000000e+00, 2.50490522e+01]), 356.59193005944485]
[array([4.20842539e+00, 9.95015742e+00, 1.00178792e-03, 0.00000000e+00,
       0.00000000e+00, 2.50490522e+01]), 356.52751738918613]
[array([4.20842539e+00, 9.95015742e+00, 1.78791688e-06, 1.00000000e-03,
       0.00000000e+00, 2.50490522e+01]), 356.6747318136642]
[array([4.20842539e+00, 9.95015742e+00, 1.78791688e-06, 0.00000000e+00,
       1.00000000e-03, 2.50490522e+01]), 356.58023525454445]
[array([4.20842539e+00, 9.95015742e+00, 1.78791688e-06, 0.00000000e+00,
       0.00000000e+00, 2.50500522e+01]), 356.5745878688081]
[array([4.20842446e+00, 9.95015256e+00, 5.39271489e-09, 0.00000000e+00,
       0.00000000e+00, 2.50490760e+01]), 356.56911202835266]
[array([4.20942446e+00, 9.95015256e+00, 5.39271489e-09, 0.00000000e+00,
       0.00000000e+00, 2.50490760e+01]), 356.56997539555795]
[array([4.20842446e+00, 9.95115256e+00, 5.39271489e-09, 0.00000000e+00,
       0.00000000e+00, 2.50490760e+01]), 356.5905681179008]
[array([4.20842446e+00, 9.95015256e+00, 1.00000539e-03, 0.00000000e+00,
       0.00000000e+00, 2.50490760e+01]), 356.52676318403127]
[array([4.20842446e+00, 9.95015256e+00, 5.39271489e-09, 1.00000000e-03,
       0.00000000e+00, 2.50490760e+01]), 356.6742088955648]
[array([4.20842446e+00, 9.95015256e+00, 5.39271489e-09, 0.00000000e+00,
       1.00000000e-03, 2.50490760e+01]), 356.57991974281435]
[array([4.20842446e+00, 9.95015256e+00, 5.39271489e-09, 0.00000000e+00,
       0.00000000e+00, 2.50500760e+01]), 356.5743395309874]
[array([4.20842446e+00, 9.95015255e+00, 7.54837539e-14, 0.00000000e+00,
       0.00000000e+00, 2.50490761e+01]), 356.5689418457048]
[array([4.20942446e+00, 9.95015255e+00, 7.54837539e-14, 0.00000000e+00,
       0.00000000e+00, 2.50490761e+01]), 356.5698894054771]
[array([4.20842446e+00, 9.95115255e+00, 7.54837539e-14, 0.00000000e+00,
       0.00000000e+00, 2.50490761e+01]), 356.5905243904706]
[array([4.20842446e+00, 9.95015255e+00, 1.00000000e-03, 0.00000000e+00,
       0.00000000e+00, 2.50490761e+01]), 356.5267445310712]
[array([4.20842446e+00, 9.95015255e+00, 7.54837539e-14, 1.00000000e-03,
       0.00000000e+00, 2.50490761e+01]), 356.6741990480748]
[array([4.20842446e+00, 9.95015255e+00, 7.54837539e-14, 0.00000000e+00,
       1.00000000e-03, 2.50490761e+01]), 356.5799169951647]
[array([4.20842446e+00, 9.95015255e+00, 7.54837539e-14, 0.00000000e+00,
       0.00000000e+00, 2.50500761e+01]), 356.57433877256926]
[array([4.20842446e+00, 9.95015255e+00, 1.50818994e-23, 0.00000000e+00,
       0.00000000e+00, 2.50490761e+01]), 356.56894230327714]
[array([4.20942446e+00, 9.95015255e+00, 1.50818994e-23, 0.00000000e+00,
       0.00000000e+00, 2.50490761e+01]), 356.56989071081847]
[array([4.20842446e+00, 9.95115255e+00, 1.50818994e-23, 0.00000000e+00,
       0.00000000e+00, 2.50490761e+01]), 356.59052584923467]
[array([4.20842446e+00, 9.95015255e+00, 1.00000000e-03, 0.00000000e+00,
       0.00000000e+00, 2.50490761e+01]), 356.5267459179233]
[array([4.20842446e+00, 9.95015255e+00, 1.50818994e-23, 1.00000000e-03,
       0.00000000e+00, 2.50490761e+01]), 356.674200155439]
[array([4.20842446e+00, 9.95015255e+00, 1.50818994e-23, 0.00000000e+00,
       1.00000000e-03, 2.50490761e+01]), 356.57991789103096]
[array([4.20842446e+00, 9.95015255e+00, 1.50818994e-23, 0.00000000e+00,
       0.00000000e+00, 2.50500761e+01]), 356.57433952291626]
[array([ 4.20842446,  9.95015255,  0.        ,  0.        ,  0.        ,
       25.04907605]), 356.568942923869]
[array([ 4.20942446,  9.95015255,  0.        ,  0.        ,  0.        ,
       25.04907605]), 356.5698911691769]
[array([ 4.20842446,  9.95115255,  0.        ,  0.        ,  0.        ,
       25.04907605]), 356.5905261815787]
[array([4.20842446e+00, 9.95015255e+00, 1.00000000e-03, 0.00000000e+00,
       0.00000000e+00, 2.50490761e+01]), 356.52674615661664]
[array([4.20842446e+00, 9.95015255e+00, 0.00000000e+00, 1.00000000e-03,
       0.00000000e+00, 2.50490761e+01]), 356.67420032358564]
[array([4.20842446e+00, 9.95015255e+00, 0.00000000e+00, 0.00000000e+00,
       1.00000000e-03, 2.50490761e+01]), 356.5799180092802]
[array([ 4.20842446,  9.95015255,  0.        ,  0.        ,  0.        ,
       25.05007605]), 356.5743396158783]
Out[20]:
425.7770085334778
In [21]:
# save the output
outfile = open("res_init_wt_bfgs.pickle", "wb")
pickle.dump(res_init_wt, outfile)
outfile.close()

# remember to write code that reads in the above output later and comment out the save
In [22]:
obs = np.sum(J)

# approximate optimal weighting matrix
params_2 = res_init_wt.x

# calculate mean utility given the optimal parameters (with id weighting matrix)
d.delta = solve_delta(df.share.values, X, V, p, y_it, 
                        d.delta, params_2, J, T, N, 1e-5)

# calculate probabilities and shares given the optimal params (w/ id weight matrix)
q_s = compute_share(X, V, p, y_it, d.delta, params_2, J, T, N)

# calculate marginal costs
mc = calc_mc(q_s, firms, p, y_it, params_2[5], J, T, N, marks, markets).reshape((-1,1))

y2 = np.vstack((d.delta.reshape((-1,1)), np.log(mc)))
x = scipy.linalg.block_diag(X,W)
z = scipy.linalg.block_diag(Z,Z_s)

# this is the first order condition that solves for the linear parameters
b = np.linalg.inv(x.T @ z @ w1 @ z.T @ x) @ (x.T @ z @ w1 @ z.T @ y2)

# obtain the error
xi_w = y2 - x @ b


# update weighting matrix
g_ind = z * xi_w
vg = g_ind.T @ g_ind / obs

# obtain optimal weighting matrix
weight = scipy.linalg.inv(vg)
In [23]:
# now search for optimal parameters with the optimal weighting matrix
t1 = time.time()

res = minimize(objective,
              theta_2, 
              args = (df.share.values, X, V, p, y_it, 
                      J, T, N, marks, markets, 1e-4, 
                      Z, Z_s, W, weight, firms), 
              bounds = bnds,
              method = "L-BFGS-B",
              options = {'maxiter': 1000, 'maxfun': 1000, 'eps': 1e-3},
              tol = 1e-4)

time.time() - t1
[array([ 3.612,  4.628,  1.818,  1.05 ,  2.056, 43.501]), 1697266.2310209193]
[array([ 3.613,  4.628,  1.818,  1.05 ,  2.056, 43.501]), 1697286.901920809]
[array([ 3.612,  4.629,  1.818,  1.05 ,  2.056, 43.501]), 1697293.6694848957]
[array([ 3.612,  4.628,  1.819,  1.05 ,  2.056, 43.501]), 1697418.1729533367]
[array([ 3.612,  4.628,  1.818,  1.051,  2.056, 43.501]), 1698388.379205081]
[array([ 3.612,  4.628,  1.818,  1.05 ,  2.057, 43.501]), 1697231.0449830429]
[array([ 3.612,  4.628,  1.818,  1.05 ,  2.056, 43.502]), 1697308.0323311102]
[array([ 3.61189735,  4.62786847,  1.81794833,  1.04997016,  3.05599939,
       43.49990579]), 1672981.7068637346]
[array([ 3.61289735,  4.62786847,  1.81794833,  1.04997016,  3.05599939,
       43.49990579]), 1673024.999877508]
[array([ 3.61189735,  4.62886847,  1.81794833,  1.04997016,  3.05599939,
       43.49990579]), 1673035.5912970707]
[array([ 3.61189735,  4.62786847,  1.81894833,  1.04997016,  3.05599939,
       43.49990579]), 1673170.7675236792]
[array([ 3.61189735,  4.62786847,  1.81794833,  1.05097016,  3.05599939,
       43.49990579]), 1674189.016456585]
[array([ 3.61189735,  4.62786847,  1.81794833,  1.04997016,  3.05699939,
       43.49990579]), 1673050.5403659765]
[array([ 3.61189735,  4.62786847,  1.81794833,  1.04997016,  3.05599939,
       43.50090579]), 1673100.632818547]
[array([ 3.61196515,  4.62795535,  1.81798246,  1.04998987,  2.39546463,
       43.50062855]), 1680405.9889845813]
[array([ 3.61296515,  4.62795535,  1.81798246,  1.04998987,  2.39546463,
       43.50062855]), 1680415.1050752148]
[array([ 3.61196515,  4.62895535,  1.81798246,  1.04998987,  2.39546463,
       43.50062855]), 1680419.0397919223]
[array([ 3.61196515,  4.62795535,  1.81898246,  1.04998987,  2.39546463,
       43.50062855]), 1680538.4359266688]
[array([ 3.61196515,  4.62795535,  1.81798246,  1.05098987,  2.39546463,
       43.50062855]), 1681531.71830681]
[array([ 3.61196515,  4.62795535,  1.81798246,  1.04998987,  2.39646463,
       43.50062855]), 1680366.244283367]
[array([ 3.61196515,  4.62795535,  1.81798246,  1.04998987,  2.39546463,
       43.50162855]), 1680434.1274675922]
[array([ 3.61191923,  4.62789651,  1.81795935,  1.04997652,  2.84283363,
       43.50013904]), 1671518.3752811428]
[array([ 3.61291923,  4.62789651,  1.81795935,  1.04997652,  2.84283363,
       43.50013904]), 1671548.3091079644]
[array([ 3.61191923,  4.62889651,  1.81795935,  1.04997652,  2.84283363,
       43.50013904]), 1671565.7520418093]
[array([ 3.61191923,  4.62789651,  1.81895935,  1.04997652,  2.84283363,
       43.50013904]), 1671701.595428446]
[array([ 3.61191923,  4.62789651,  1.81795935,  1.05097652,  2.84283363,
       43.50013904]), 1672719.31760592]
[array([ 3.61191923,  4.62789651,  1.81795935,  1.04997652,  2.84383363,
       43.50013904]), 1671573.7299728075]
[array([ 3.61191923,  4.62789651,  1.81795935,  1.04997652,  2.84283363,
       43.50113904]), 1671626.252477229]
[array([ 3.61194954,  4.62793534,  1.8179746 ,  1.04998533,  2.54757009,
       43.50046212]), 1675437.260468715]
[array([ 3.61294954,  4.62793534,  1.8179746 ,  1.04998533,  2.54757009,
       43.50046212]), 1675442.7356603765]
[array([ 3.61194954,  4.62893534,  1.8179746 ,  1.04998533,  2.54757009,
       43.50046212]), 1675445.3642414098]
[array([ 3.61194954,  4.62793534,  1.8189746 ,  1.04998533,  2.54757009,
       43.50046212]), 1675563.181139877]
[array([ 3.61194954,  4.62793534,  1.8179746 ,  1.05098533,  2.54757009,
       43.50046212]), 1676562.6789916446]
[array([ 3.61194954,  4.62793534,  1.8179746 ,  1.04998533,  2.54857009,
       43.50046212]), 1675399.8612780077]
[array([ 3.61194954,  4.62793534,  1.8179746 ,  1.04998533,  2.54757009,
       43.50146212]), 1675461.315627017]
[array([ 3.61192822,  4.62790803,  1.81796387,  1.04997913,  2.75525705,
       43.50023486]), 1671859.453815337]
[array([ 3.61292822,  4.62790803,  1.81796387,  1.04997913,  2.75525705,
       43.50023486]), 1671887.8371488573]
[array([ 3.61192822,  4.62890803,  1.81796387,  1.04997913,  2.75525705,
       43.50023486]), 1671905.4661293754]
[array([ 3.61192822,  4.62790803,  1.81896387,  1.04997913,  2.75525705,
       43.50023486]), 1672040.2498468019]
[array([ 3.61192822,  4.62790803,  1.81796387,  1.05097913,  2.75525705,
       43.50023486]), 1673056.098632778]
[array([ 3.61192822,  4.62790803,  1.81796387,  1.04997913,  2.75625705,
       43.50023486]), 1671906.5776074864]
[array([ 3.61192822,  4.62790803,  1.81796387,  1.04997913,  2.75525705,
       43.50123486]), 1671961.5054070854]
[array([ 3.61192106,  4.62789886,  1.81796027,  1.04997705,  2.82498865,
       43.50015856]), 1671544.0238932627]
[array([ 3.61292106,  4.62789886,  1.81796027,  1.04997705,  2.82498865,
       43.50015856]), 1671572.5550322891]
[array([ 3.61192106,  4.62889886,  1.81796027,  1.04997705,  2.82498865,
       43.50015856]), 1671588.8833299351]
[array([ 3.61192106,  4.62789886,  1.81896027,  1.04997705,  2.82498865,
       43.50015856]), 1671722.8133868838]
[array([ 3.61192106,  4.62789886,  1.81796027,  1.05097705,  2.82498865,
       43.50015856]), 1672738.8261884947]
[array([ 3.61192106,  4.62789886,  1.81796027,  1.04997705,  2.82598865,
       43.50015856]), 1671591.6590438234]
[array([ 3.61192106,  4.62789886,  1.81796027,  1.04997705,  2.82498865,
       43.50115856]), 1671644.3582246203]
[array([ 3.61191923,  4.62789651,  1.81795935,  1.04997652,  2.84283363,
       43.50013904]), 1671514.4397682103]
[array([ 3.61291923,  4.62789651,  1.81795935,  1.04997652,  2.84283363,
       43.50013904]), 1671548.8157322267]
[array([ 3.61191923,  4.62889651,  1.81795935,  1.04997652,  2.84283363,
       43.50013904]), 1671567.9642323365]
[array([ 3.61191923,  4.62789651,  1.81895935,  1.04997652,  2.84283363,
       43.50013904]), 1671704.124790952]
[array([ 3.61191923,  4.62789651,  1.81795935,  1.05097652,  2.84283363,
       43.50013904]), 1672721.089078341]
[array([ 3.61191923,  4.62789651,  1.81795935,  1.04997652,  2.84383363,
       43.50013904]), 1671574.7729091577]
[array([ 3.61191923,  4.62789651,  1.81795935,  1.04997652,  2.84283363,
       43.50113904]), 1671627.0278715107]
[array([ 3.52984239,  4.50966098,  1.28939331,  0.        ,  3.690032  ,
       43.28271412]), 1029104.0167804187]
[array([ 3.53084239,  4.50966098,  1.28939331,  0.        ,  3.690032  ,
       43.28271412]), 1029117.6661600292]
[array([ 3.52984239,  4.51066098,  1.28939331,  0.        ,  3.690032  ,
       43.28271412]), 1029126.6983883999]
[array([ 3.52984239,  4.50966098,  1.29039331,  0.        ,  3.690032  ,
       43.28271412]), 1029173.1048616328]
[array([3.52984239e+00, 4.50966098e+00, 1.28939331e+00, 1.00000000e-03,
       3.69003200e+00, 4.32827141e+01]), 1029238.6440764571]
[array([ 3.52984239,  4.50966098,  1.28939331,  0.        ,  3.691032  ,
       43.28271412]), 1029185.9791082488]
[array([ 3.52984239,  4.50966098,  1.28939331,  0.        ,  3.690032  ,
       43.28371412]), 1029200.4941830151]
[array([ 3.52256154,  4.50167989,  1.21658533,  0.        ,  3.05393313,
       43.23201332]), 1032944.9142303743]
[array([ 3.52356154,  4.50167989,  1.21658533,  0.        ,  3.05393313,
       43.23201332]), 1032944.5139521437]
[array([ 3.52256154,  4.50267989,  1.21658533,  0.        ,  3.05393313,
       43.23201332]), 1032951.3574971608]
[array([ 3.52256154,  4.50167989,  1.21758533,  0.        ,  3.05393313,
       43.23201332]), 1032995.3075386124]
[array([3.52256154e+00, 4.50167989e+00, 1.21658533e+00, 1.00000000e-03,
       3.05393313e+00, 4.32320133e+01]), 1033052.0000959481]
[array([ 3.52256154,  4.50167989,  1.21658533,  0.        ,  3.05493313,
       43.23201332]), 1032999.0192450467]
[array([ 3.52256154,  4.50167989,  1.21658533,  0.        ,  3.05393313,
       43.23301332]), 1033022.681896546]
[array([ 3.52823455,  4.5078985 ,  1.27331497,  0.        ,  3.54956095,
       43.27151776]), 1026349.4068227457]
[array([ 3.52923455,  4.5078985 ,  1.27331497,  0.        ,  3.54956095,
       43.27151776]), 1026356.7860688473]
[array([ 3.52823455,  4.5088985 ,  1.27331497,  0.        ,  3.54956095,
       43.27151776]), 1026363.2772060981]
[array([ 3.52823455,  4.5078985 ,  1.27431497,  0.        ,  3.54956095,
       43.27151776]), 1026406.0793182249]
[array([3.52823455e+00, 4.50789850e+00, 1.27331497e+00, 1.00000000e-03,
       3.54956095e+00, 4.32715178e+01]), 1026468.7934783362]
[array([ 3.52823455,  4.5078985 ,  1.27331497,  0.        ,  3.55056095,
       43.27151776]), 1026403.6902427077]
[array([ 3.52823455,  4.5078985 ,  1.27331497,  0.        ,  3.54956095,
       43.27251776]), 1026419.781168961]
[array([ 3.53441428,  4.50899222,  1.15081984,  0.        ,  3.33275548,
       43.19921612]), 1027607.52050824]
[array([ 3.53541428,  4.50899222,  1.15081984,  0.        ,  3.33275548,
       43.19921612]), 1027607.1656692743]
[array([ 3.53441428,  4.50999222,  1.15081984,  0.        ,  3.33275548,
       43.19921612]), 1027608.7743588301]
[array([ 3.53441428,  4.50899222,  1.15181984,  0.        ,  3.33275548,
       43.19921612]), 1027643.9168713412]
[array([3.53441428e+00, 4.50899222e+00, 1.15081984e+00, 1.00000000e-03,
       3.33275548e+00, 4.31992161e+01]), 1027702.8761864667]
[array([ 3.53441428,  4.50899222,  1.15081984,  0.        ,  3.33375548,
       43.19921612]), 1027653.6376705859]
[array([ 3.53441428,  4.50899222,  1.15081984,  0.        ,  3.33275548,
       43.20021612]), 1027672.2219677776]
[array([ 3.52954573,  4.50813056,  1.24732465,  0.        ,  3.5035604 ,
       43.25617721]), 1024479.8659670524]
[array([ 3.53054573,  4.50813056,  1.24732465,  0.        ,  3.5035604 ,
       43.25617721]), 1024489.5085782404]
[array([ 3.52954573,  4.50913056,  1.24732465,  0.        ,  3.5035604 ,
       43.25617721]), 1024499.4692165486]
[array([ 3.52954573,  4.50813056,  1.24832465,  0.        ,  3.5035604 ,
       43.25617721]), 1024542.7793915409]
[array([3.52954573e+00, 4.50813056e+00, 1.24732465e+00, 1.00000000e-03,
       3.50356040e+00, 4.32561772e+01]), 1024606.617693387]
[array([ 3.52954573,  4.50813056,  1.24732465,  0.        ,  3.5045604 ,
       43.25617721]), 1024540.3299844918]
[array([ 3.52954573,  4.50813056,  1.24732465,  0.        ,  3.5035604 ,
       43.25717721]), 1024556.8937515435]
[array([ 3.53048437,  4.50829669,  1.22871877,  0.        ,  3.47062965,
       43.24519526]), 1028048.1157681576]
[array([ 3.53148437,  4.50829669,  1.22871877,  0.        ,  3.47062965,
       43.24519526]), 1028058.2797621175]
[array([ 3.53048437,  4.50929669,  1.22871877,  0.        ,  3.47062965,
       43.24519526]), 1027966.1861238715]
[array([ 3.53048437,  4.50829669,  1.22971877,  0.        ,  3.47062965,
       43.24519526]), 1027924.4281771705]
[array([3.53048437e+00, 4.50829669e+00, 1.22871877e+00, 1.00000000e-03,
       3.47062965e+00, 4.32451953e+01]), 1027882.8343640415]
[array([ 3.53048437,  4.50829669,  1.22871877,  0.        ,  3.47162965,
       43.24519526]), 1028099.1473214473]
[array([ 3.53048437,  4.50829669,  1.22871877,  0.        ,  3.47062965,
       43.24619526]), 1028096.3643808393]
[array([ 3.52965857,  4.50815053,  1.24508787,  0.        ,  3.4996015 ,
       43.25485697]), 1024312.1994143518]
[array([ 3.53065857,  4.50815053,  1.24508787,  0.        ,  3.4996015 ,
       43.25485697]), 1024326.3816745884]
[array([ 3.52965857,  4.50915053,  1.24508787,  0.        ,  3.4996015 ,
       43.25485697]), 1024339.2866495214]
[array([ 3.52965857,  4.50815053,  1.24608787,  0.        ,  3.4996015 ,
       43.25485697]), 1024384.8792831269]
[array([3.52965857e+00, 4.50815053e+00, 1.24508787e+00, 1.00000000e-03,
       3.49960150e+00, 4.32548570e+01]), 1024450.2774225465]
[array([ 3.52965857,  4.50815053,  1.24508787,  0.        ,  3.5006015 ,
       43.25485697]), 1024384.5923546111]
[array([ 3.52965857,  4.50815053,  1.24508787,  0.        ,  3.4996015 ,
       43.25585697]), 1024401.4270352563]
[array([ 3.52975702,  4.50816796,  1.24313648,  0.        ,  3.49614771,
       43.25370518]), 1024247.6852064499]
[array([ 3.53075702,  4.50816796,  1.24313648,  0.        ,  3.49614771,
       43.25370518]), 1024233.5739130742]
[array([ 3.52975702,  4.50916796,  1.24313648,  0.        ,  3.49614771,
       43.25370518]), 1024231.2786427821]
[array([ 3.52975702,  4.50816796,  1.24413648,  0.        ,  3.49614771,
       43.25370518]), 1024262.8272522852]
[array([3.52975702e+00, 4.50816796e+00, 1.24313648e+00, 1.00000000e-03,
       3.49614771e+00, 4.32537052e+01]), 1024320.4589923982]
[array([ 3.52975702,  4.50816796,  1.24313648,  0.        ,  3.49714771,
       43.25370518]), 1024258.1503032175]
[array([ 3.52975702,  4.50816796,  1.24313648,  0.        ,  3.49614771,
       43.25470518]), 1024272.3122774708]
[array([ 3.54048112,  4.52041967,  1.22474192,  0.        ,  3.48284256,
       43.23763903]), 1026325.2911288654]
[array([ 3.54148112,  4.52041967,  1.22474192,  0.        ,  3.48284256,
       43.23763903]), 1026363.2022073168]
[array([ 3.54048112,  4.52141967,  1.22474192,  0.        ,  3.48284256,
       43.23763903]), 1025967.4466040425]
[array([ 3.54048112,  4.52041967,  1.22574192,  0.        ,  3.48284256,
       43.23763903]), 1025579.4806804896]
[array([3.54048112e+00, 4.52041967e+00, 1.22474192e+00, 1.00000000e-03,
       3.48284256e+00, 4.32376390e+01]), 1024862.4469180055]
[array([ 3.54048112,  4.52041967,  1.22474192,  0.        ,  3.48384256,
       43.23763903]), 1026379.4652671984]
[array([ 3.54048112,  4.52041967,  1.22474192,  0.        ,  3.48284256,
       43.23863903]), 1026320.3568329145]
[array([ 3.53271193,  4.51154379,  1.23806805,  0.        ,  3.49248162,
       43.24927832]), 1023966.8031515223]
[array([ 3.53371193,  4.51154379,  1.23806805,  0.        ,  3.49248162,
       43.24927832]), 1023984.5474939578]
[array([ 3.53271193,  4.51254379,  1.23806805,  0.        ,  3.49248162,
       43.24927832]), 1023999.991661148]
[array([ 3.53271193,  4.51154379,  1.23906805,  0.        ,  3.49248162,
       43.24927832]), 1024047.6210730525]
[array([3.53271193e+00, 4.51154379e+00, 1.23806805e+00, 1.00000000e-03,
       3.49248162e+00, 4.32492783e+01]), 1024114.7737464573]
[array([ 3.53271193,  4.51154379,  1.23806805,  0.        ,  3.49348162,
       43.24927832]), 1024049.9395664802]
[array([ 3.53271193,  4.51154379,  1.23806805,  0.        ,  3.49248162,
       43.25027832]), 1024067.2158157988]
[array([ 3.53457687,  4.51367438,  1.23486921,  0.        ,  3.49016783,
       43.24648439]), 1023936.2722729219]
[array([ 3.53557687,  4.51367438,  1.23486921,  0.        ,  3.49016783,
       43.24648439]), 1023918.2476570662]
[array([ 3.53457687,  4.51467438,  1.23486921,  0.        ,  3.49016783,
       43.24648439]), 1023914.2409029156]
[array([ 3.53457687,  4.51367438,  1.23586921,  0.        ,  3.49016783,
       43.24648439]), 1023943.2397936321]
[array([3.53457687e+00, 4.51367438e+00, 1.23486921e+00, 1.00000000e-03,
       3.49016783e+00, 4.32464844e+01]), 1023999.243530603]
[array([ 3.53457687,  4.51367438,  1.23486921,  0.        ,  3.49116783,
       43.24648439]), 1023936.0925802401]
[array([ 3.53457687,  4.51367438,  1.23486921,  0.        ,  3.49016783,
       43.24748439]), 1023949.8726371753]
[array([ 3.56376867,  4.54759313,  1.19635256,  0.        ,  3.49502266,
       43.20977627]), 1028111.7004573392]
[array([ 3.56476867,  4.54759313,  1.19635256,  0.        ,  3.49502266,
       43.20977627]), 1028117.9948500774]
[array([ 3.56376867,  4.54859313,  1.19635256,  0.        ,  3.49502266,
       43.20977627]), 1028046.7856984123]
[array([ 3.56376867,  4.54759313,  1.19735256,  0.        ,  3.49502266,
       43.20977627]), 1028025.0891668407]
[array([3.56376867e+00, 4.54759313e+00, 1.19635256e+00, 1.00000000e-03,
       3.49502266e+00, 4.32097763e+01]), 1028011.8224846026]
[array([ 3.56376867,  4.54759313,  1.19635256,  0.        ,  3.49602266,
       43.20977627]), 1028143.6434923452]
[array([ 3.56376867,  4.54759313,  1.19635256,  0.        ,  3.49502266,
       43.21077627]), 1028145.2189025714]
[array([ 3.53651782,  4.51592962,  1.23230825,  0.        ,  3.49049063,
       43.24404368]), 1023801.9835210122]
[array([ 3.53751782,  4.51592962,  1.23230825,  0.        ,  3.49049063,
       43.24404368]), 1023806.900304563]
[array([ 3.53651782,  4.51692962,  1.23230825,  0.        ,  3.49049063,
       43.24404368]), 1023813.4610533181]
[array([ 3.53651782,  4.51592962,  1.23330825,  0.        ,  3.49049063,
       43.24404368]), 1023856.309194048]
[array([3.53651782e+00, 4.51592962e+00, 1.23230825e+00, 1.00000000e-03,
       3.49049063e+00, 4.32440437e+01]), 1023921.6995449837]
[array([ 3.53651782,  4.51592962,  1.23230825,  0.        ,  3.49149063,
       43.24404368]), 1023862.6576988089]
[array([ 3.53651782,  4.51592962,  1.23230825,  0.        ,  3.49049063,
       43.24504368]), 1023878.3324204369]
[array([ 3.53932794,  4.51919478,  1.22860049,  0.        ,  3.49095797,
       43.24051001]), 1023710.5895524976]
[array([ 3.54032794,  4.51919478,  1.22860049,  0.        ,  3.49095797,
       43.24051001]), 1023716.6117002342]
[array([ 3.53932794,  4.52019478,  1.22860049,  0.        ,  3.49095797,
       43.24051001]), 1023724.711126384]
[array([ 3.53932794,  4.51919478,  1.22960049,  0.        ,  3.49095797,
       43.24051001]), 1023766.1534342794]
[array([3.53932794e+00, 4.51919478e+00, 1.22860049e+00, 1.00000000e-03,
       3.49095797e+00, 4.32405100e+01]), 1023829.3700070523]
[array([ 3.53932794,  4.51919478,  1.22860049,  0.        ,  3.49195797,
       43.24051001]), 1023762.5958171049]
[array([ 3.53932794,  4.51919478,  1.22860049,  0.        ,  3.49095797,
       43.24151001]), 1023779.1027477634]
[array([ 3.55154831,  4.53339395,  1.21247652,  0.        ,  3.49299031,
       43.22514314]), 1027107.2131651404]
[array([ 3.55254831,  4.53339395,  1.21247652,  0.        ,  3.49299031,
       43.22514314]), 1027111.8935391384]
[array([ 3.55154831,  4.53439395,  1.21247652,  0.        ,  3.49299031,
       43.22514314]), 1026905.620868828]
[array([ 3.55154831,  4.53339395,  1.21347652,  0.        ,  3.49299031,
       43.22514314]), 1026764.6878732713]
[array([3.55154831e+00, 4.53339395e+00, 1.21247652e+00, 1.00000000e-03,
       3.49299031e+00, 4.32251431e+01]), 1026578.1190831733]
[array([ 3.55154831,  4.53339395,  1.21247652,  0.        ,  3.49399031,
       43.22514314]), 1027102.0311796132]
[array([ 3.55154831,  4.53339395,  1.21247652,  0.        ,  3.49299031,
       43.22614314]), 1027080.8653259344]
[array([ 3.54032964,  4.52035868,  1.22727881,  0.        ,  3.49112456,
       43.2392504 ]), 1023678.6902814157]
[array([ 3.54132964,  4.52035868,  1.22727881,  0.        ,  3.49112456,
       43.2392504 ]), 1023680.0125453514]
[array([ 3.54032964,  4.52135868,  1.22727881,  0.        ,  3.49112456,
       43.2392504 ]), 1023685.905110815]
[array([ 3.54032964,  4.52035868,  1.22827881,  0.        ,  3.49112456,
       43.2392504 ]), 1023726.6972281262]
[array([3.54032964e+00, 4.52035868e+00, 1.22727881e+00, 1.00000000e-03,
       3.49112456e+00, 4.32392504e+01]), 1023790.7237571832]
[array([ 3.54032964,  4.52035868,  1.22727881,  0.        ,  3.49212456,
       43.2392504 ]), 1023731.1833851463]
[array([ 3.54032964,  4.52035868,  1.22727881,  0.        ,  3.49112456,
       43.2402504 ]), 1023746.5567293942]
[array([ 3.54117147,  4.52133682,  1.22616807,  0.        ,  3.49126457,
       43.23819181]), 1023674.6615615797]
[array([ 3.54217147,  4.52133682,  1.22616807,  0.        ,  3.49126457,
       43.23819181]), 1023669.9394634217]
[array([ 3.54117147,  4.52233682,  1.22616807,  0.        ,  3.49126457,
       43.23819181]), 1023671.4261852528]
[array([ 3.54117147,  4.52133682,  1.22716807,  0.        ,  3.49126457,
       43.23819181]), 1023707.6707529848]
[array([3.54117147e+00, 4.52133682e+00, 1.22616807e+00, 1.00000000e-03,
       3.49126457e+00, 4.32381918e+01]), 1023767.7287819097]
[array([ 3.54117147,  4.52133682,  1.22616807,  0.        ,  3.49226457,
       43.23819181]), 1023699.6079096551]
[array([ 3.54117147,  4.52133682,  1.22616807,  0.        ,  3.49126457,
       43.23919181]), 1023715.6764809695]
[array([ 3.54635989,  4.52736539,  1.2193223 ,  0.        ,  3.49212744,
       43.23166748]), 1025614.407644709]
[array([ 3.54735989,  4.52736539,  1.2193223 ,  0.        ,  3.49212744,
       43.23166748]), 1025599.6131008514]
[array([ 3.54635989,  4.52836539,  1.2193223 ,  0.        ,  3.49212744,
       43.23166748]), 1024703.839306823]
[array([ 3.54635989,  4.52736539,  1.2203223 ,  0.        ,  3.49212744,
       43.23166748]), 1022989.6354692044]
[array([3.54635989e+00, 4.52736539e+00, 1.21932230e+00, 1.00000000e-03,
       3.49212744e+00, 4.32316675e+01]), 1023584.6202625274]
[array([ 3.54635989,  4.52736539,  1.2193223 ,  0.        ,  3.49312744,
       43.23166748]), 1025424.3447679423]
[array([ 3.54635989,  4.52736539,  1.2193223 ,  0.        ,  3.49212744,
       43.23266748]), 1025298.4512266654]
[array([ 3.54282907,  4.52326284,  1.22398097,  0.        ,  3.49154024,
       43.23610741]), 1023609.566511323]
[array([ 3.54382907,  4.52326284,  1.22398097,  0.        ,  3.49154024,
       43.23610741]), 1023606.8016021943]
[array([ 3.54282907,  4.52426284,  1.22398097,  0.        ,  3.49154024,
       43.23610741]), 1023609.5575378061]
[array([ 3.54282907,  4.52326284,  1.22498097,  0.        ,  3.49154024,
       43.23610741]), 1023645.9845143167]
[array([3.54282907e+00, 4.52326284e+00, 1.22398097e+00, 1.00000000e-03,
       3.49154024e+00, 4.32361074e+01]), 1023707.0644952692]
[array([ 3.54282907,  4.52326284,  1.22398097,  0.        ,  3.49254024,
       43.23610741]), 1023646.2952569533]
[array([ 3.54282907,  4.52326284,  1.22398097,  0.        ,  3.49154024,
       43.23710741]), 1023661.0462947906]
[array([ 3.54432864,  4.52500522,  1.2220024 ,  0.        ,  3.49178963,
       43.23422173]), 1021619.6694257593]
[array([ 3.54532864,  4.52500522,  1.2220024 ,  0.        ,  3.49178963,
       43.23422173]), 1020702.1044906578]
[array([ 3.54432864,  4.52600522,  1.2220024 ,  0.        ,  3.49178963,
       43.23422173]), 1023557.7835348978]
[array([ 3.54432864,  4.52500522,  1.2230024 ,  0.        ,  3.49178963,
       43.23422173]), 1023595.9656362459]
[array([3.54432864e+00, 4.52500522e+00, 1.22200240e+00, 1.00000000e-03,
       3.49178963e+00, 4.32342217e+01]), 1023657.2791487683]
[array([ 3.54432864,  4.52500522,  1.2220024 ,  0.        ,  3.49278963,
       43.23422173]), 1023589.780571027]
[array([ 3.54432864,  4.52500522,  1.2220024 ,  0.        ,  3.49178963,
       43.23522173]), 1023605.9714446875]
[array([ 3.54473288,  4.52547492,  1.22146903,  0.        ,  3.49185685,
       43.23371342]), 1023071.8080431024]
[array([ 3.54573288,  4.52547492,  1.22146903,  0.        ,  3.49185685,
       43.23371342]), 1023341.1480012906]
[array([ 3.54473288,  4.52647492,  1.22146903,  0.        ,  3.49185685,
       43.23371342]), 1023553.2412444146]
[array([ 3.54473288,  4.52547492,  1.22246903,  0.        ,  3.49185685,
       43.23371342]), 1023586.5234472909]
[array([3.54473288e+00, 4.52547492e+00, 1.22146903e+00, 1.00000000e-03,
       3.49185685e+00, 4.32337134e+01]), 1023645.1295174757]
[array([ 3.54473288,  4.52547492,  1.22146903,  0.        ,  3.49285685,
       43.23371342]), 1021957.6242148378]
[array([ 3.54473288,  4.52547492,  1.22146903,  0.        ,  3.49185685,
       43.23471342]), 1023592.3222085705]
[array([ 3.54432864,  4.52500522,  1.2220024 ,  0.        ,  3.49178963,
       43.23422173]), 1023588.9356674746]
[array([ 3.54532864,  4.52500522,  1.2220024 ,  0.        ,  3.49178963,
       43.23422173]), 1023579.130545111]
[array([ 3.54432864,  4.52600522,  1.2220024 ,  0.        ,  3.49178963,
       43.23422173]), 1023576.441132134]
[array([ 3.54432864,  4.52500522,  1.2230024 ,  0.        ,  3.49178963,
       43.23422173]), 1023605.3396744936]
[array([3.54432864e+00, 4.52500522e+00, 1.22200240e+00, 1.00000000e-03,
       3.49178963e+00, 4.32342217e+01]), 1023661.1422494783]
[array([ 3.54432864,  4.52500522,  1.2220024 ,  0.        ,  3.49278963,
       43.23422173]), 1023598.0284311832]
[array([ 3.54432864,  4.52500522,  1.2220024 ,  0.        ,  3.49178963,
       43.23522173]), 1023611.6189528179]
[array([ 3.57356456,  4.55921678,  1.17817244,  0.        ,  3.49391182,
       43.19422931]), 1028516.8399457252]
[array([ 3.57456456,  4.55921678,  1.17817244,  0.        ,  3.49391182,
       43.19422931]), 1028522.1245698994]
[array([ 3.57356456,  4.56021678,  1.17817244,  0.        ,  3.49391182,
       43.19422931]), 1028489.8707740743]
[array([ 3.57356456,  4.55921678,  1.17917244,  0.        ,  3.49391182,
       43.19422931]), 1028498.5119065555]
[array([3.57356456e+00, 4.55921678e+00, 1.17817244e+00, 1.00000000e-03,
       3.49391182e+00, 4.31942293e+01]), 1028524.5566041723]
[array([ 3.57356456,  4.55921678,  1.17817244,  0.        ,  3.49491182,
       43.19422931]), 1028556.8743293566]
[array([ 3.57356456,  4.55921678,  1.17817244,  0.        ,  3.49391182,
       43.19522931]), 1028565.421434999]
[array([ 3.54602233,  4.52698716,  1.21946325,  0.        ,  3.49191257,
       43.2319049 ]), 1025538.8034142067]
[array([ 3.54702233,  4.52698716,  1.21946325,  0.        ,  3.49191257,
       43.2319049 ]), 1025604.7802314135]
[array([ 3.54602233,  4.52798716,  1.21946325,  0.        ,  3.49191257,
       43.2319049 ]), 1024850.3708747209]
[array([ 3.54602233,  4.52698716,  1.22046325,  0.        ,  3.49191257,
       43.2319049 ]), 1023608.627894167]
[array([3.54602233e+00, 4.52698716e+00, 1.21946325e+00, 1.00000000e-03,
       3.49191257e+00, 4.32319049e+01]), 1023572.9351357062]
[array([ 3.54602233,  4.52698716,  1.21946325,  0.        ,  3.49291257,
       43.2319049 ]), 1025608.5401481852]
[array([ 3.54602233,  4.52698716,  1.21946325,  0.        ,  3.49191257,
       43.2329049 ]), 1025489.097084126]
[array([ 3.5443774 ,  4.52506228,  1.22192929,  0.        ,  3.49179317,
       43.23415503]), 1023574.1476076498]
[array([ 3.5453774 ,  4.52506228,  1.22192929,  0.        ,  3.49179317,
       43.23415503]), 1023567.0943754864]
[array([ 3.5443774 ,  4.52606228,  1.22192929,  0.        ,  3.49179317,
       43.23415503]), 1023566.4884147416]
[array([ 3.5443774 ,  4.52506228,  1.22292929,  0.        ,  3.49179317,
       43.23415503]), 1023598.247095055]
[array([3.54437740e+00, 4.52506228e+00, 1.22192929e+00, 1.00000000e-03,
       3.49179317e+00, 4.32341550e+01]), 1023656.1365285342]
[array([ 3.5443774 ,  4.52506228,  1.22192929,  0.        ,  3.49279317,
       43.23415503]), 1023594.0123696981]
[array([ 3.5443774 ,  4.52506228,  1.22192929,  0.        ,  3.49179317,
       43.23515503]), 1023608.1017897965]
[array([ 3.54442935,  4.52512308,  1.22185141,  0.        ,  3.49179694,
       43.23408397]), 1023588.4211776562]
[array([ 3.54542935,  4.52512308,  1.22185141,  0.        ,  3.49179694,
       43.23408397]), 1023573.2778621664]
[array([ 3.54442935,  4.52612308,  1.22185141,  0.        ,  3.49179694,
       43.23408397]), 1023570.4231933381]
[array([ 3.54442935,  4.52512308,  1.22285141,  0.        ,  3.49179694,
       43.23408397]), 1023600.0773880269]
[array([3.54442935e+00, 4.52512308e+00, 1.22185141e+00, 1.00000000e-03,
       3.49179694e+00, 4.32340840e+01]), 1023656.3866592862]
[array([ 3.54442935,  4.52512308,  1.22185141,  0.        ,  3.49279694,
       43.23408397]), 1023586.6825700551]
[array([ 3.54442935,  4.52512308,  1.22185141,  0.        ,  3.49179694,
       43.23508397]), 1023602.2121376498]
[array([ 3.54437983,  4.52506512,  1.22192566,  0.        ,  3.49179334,
       43.23415172]), 1023589.1374074619]
[array([ 3.54537983,  4.52506512,  1.22192566,  0.        ,  3.49179334,
       43.23415172]), 1023574.5790668075]
[array([ 3.54437983,  4.52606512,  1.22192566,  0.        ,  3.49179334,
       43.23415172]), 1023571.9242822158]
[array([ 3.54437983,  4.52506512,  1.22292566,  0.        ,  3.49179334,
       43.23415172]), 1023601.8468162068]
[array([3.54437983e+00, 4.52506512e+00, 1.22192566e+00, 1.00000000e-03,
       3.49179334e+00, 4.32341517e+01]), 1023658.3548652787]
[array([ 3.54437983,  4.52506512,  1.22192566,  0.        ,  3.49279334,
       43.23415172]), 1023588.7660867828]
[array([ 3.54437983,  4.52506512,  1.22192566,  0.        ,  3.49179334,
       43.23515172]), 1023604.3419187331]
[array([ 3.5443774 ,  4.52506228,  1.22192929,  0.        ,  3.49179317,
       43.23415503]), 1023588.6133227705]
[array([ 3.5453774 ,  4.52506228,  1.22192929,  0.        ,  3.49179317,
       43.23415503]), 1023574.2362814373]
[array([ 3.5443774 ,  4.52606228,  1.22192929,  0.        ,  3.49179317,
       43.23415503]), 1023571.7036171192]
[array([ 3.5443774 ,  4.52506228,  1.22292929,  0.        ,  3.49179317,
       43.23415503]), 1023601.7908459071]
[array([3.54437740e+00, 4.52506228e+00, 1.22192929e+00, 1.00000000e-03,
       3.49179317e+00, 4.32341550e+01]), 1023658.4020509939]
[array([ 3.5443774 ,  4.52506228,  1.22192929,  0.        ,  3.49279317,
       43.23415503]), 1023588.8571388279]
[array([ 3.5443774 ,  4.52506228,  1.22192929,  0.        ,  3.49179317,
       43.23515503]), 1023604.4451595861]
Out[23]:
645.4324803352356
In [24]:
outfile = open("res_bfgs.pickle", "wb")
pickle.dump(res, outfile)
outfile.close()
In [25]:
# obtain the linear parameters
params_3 = res.x

d.delta = solve_delta(df.share.values, X, V, p, y_it, 
                        d.delta, params_3, J, T, N, 1e-4)
q_s = compute_share(X, V, p, y_it, d.delta, params_3, J, T, N)
    
mc = calc_mc(q_s, firms, p, y_it, params_3[5], J, T, N, marks, markets).reshape((-1,1))
    
y2 = np.vstack((d.delta.reshape((-1,1)), np.log(mc)))
In [26]:
# linear parameters
#  this is the FOC on page 5 of Aviv's Appendix
# compare parameters to Table iv of BLP
#  first 5 are the demand side means
#  last 6 are the cost side params
b = np.linalg.inv(x.T @ z @ weight @ z.T @ x) @ (x.T @ z @ weight @ z.T @ y2)
b
Out[26]:
array([[-6.63020003],
       [ 3.14163877],
       [ 1.40442712],
       [ 0.42582024],
       [ 2.57036016],
       [ 0.9209081 ],
       [ 0.40041309],
       [ 0.52769782],
       [-0.66717413],
       [-0.37639983],
       [ 0.02182854]])
In [27]:
# non-linear parameters - these are the sigma estimates in table iv of BLP
res.x
Out[27]:
array([ 3.5443774 ,  4.52506228,  1.22192929,  0.        ,  3.49179317,
       43.23415503])
In [28]:
# average markup (compare to table 12 of Conlon and Gortmaker (2019))
np.mean((p.flatten() - mc.flatten())/p.flatten())
Out[28]:
0.31106243756696267
In [29]:
share_deriv = []
q = q_s[0]
s = q_s[1]
for m in marks:
    # obtain list of firms operating in that market/year
    firm_yr = firms[markets == m]

    # obtain list of prices of goods in that market/year
    price = p[markets == m]

    # J_t x J_t block matrix of 1's indicating goods belonging to same firm
    #  in that market/year
    # Also known as the ownership matrix
    same_firm = np.equal(firm_yr, np.transpose(firm_yr))

    # obtain matrix of probabilities for all simulations for goods in that 
    #  market/year
    yr = q[markets == m,:]

    # obtain number of goods in that market
    nobs = np.size(yr, 0)

    # this is the omega matrix initializing        
    deriv_m = np.zeros((nobs, nobs))

    # calculate the omega matix by iterating through all individuals
    #  Omega matrix is cross-price share deriv element-wise product with
    #  ownership matrix
    for i in range(N):
        yr_i = yr[:, i].reshape((-1, 1))
        deriv_m += params_3[5] / y_it[i, m - 1] * 1/N * \
        (yr_i @ yr_i.T - np.diag(yr[:,i]))
        
    share_deriv.append(deriv_m)     
In [30]:
# obtain the mazda 323 own-price deriv in 1990
mz323_deriv = (share_deriv[19][df.modelvec[markets == 20] == "MZ323"]).flatten()[df.modelvec[markets == 20] == "MZ323"]

# obtain the bmw 735i own-price deriv in 1990
bw735i_deriv = (share_deriv[19][df.modelvec[markets == 20] == "BW735i"]).flatten()[df.modelvec[markets == 20] == "BW735i"]
In [31]:
# obtain the mazda 323 mkt share in 1990
mz323_s = s[(markets == 20) & (df.modelvec == "MZ323")]

# obtain the bmw 735i market share in 1990
bw735i_s = s[(markets == 20) & (df.modelvec == "BW735i")]
In [32]:
# compare with table VI, first row first column of BLP
mz323_deriv/ mz323_s * 100
Out[32]:
array([-105.11981736])
In [33]:
# compare with table VI, last row last column of BLP
bw735i_deriv/bw735i_s * 100
Out[33]:
array([-7.93455796])
In [34]:
# now we attempt to create all of table VI

# extract all of the cars
table_vi_cars = ["MZ323", "NISENT", "FDESCO", "CVCAVA", "HDACCO", "FDTAUR", "BKCENT",
                "NIMAXI", "ACLEGE", "LNTOWN", "CDSEVI", "LXLS40", "BW735i"]
In [35]:
# obtain share derivs and shares relevant cars in 1990
deriv_1990 = share_deriv[19]
s_1990 = s[markets == 20]
In [36]:
# initialize the matrix dimensions of table VI
table_vi_mx = np.zeros((len(table_vi_cars), len(table_vi_cars)))

# iterate across cars to obtain table vi
row = 0
for car in table_vi_cars:
    col = 0
    for cars in table_vi_cars:
        table_vi_mx[row, col] = (deriv_1990[df.modelvec[markets == 20] == cars].flatten())[df.modelvec[markets == 20] == car] / \
            s_1990[df.modelvec[markets == 20] == car] * 100
        col += 1
    row += 1    
In [37]:
# store in dataframe for easier viewing
table_vi = pd.DataFrame(np.round(table_vi_mx, 3))

# renaming rows and columns
table_vi.index = table_vi_cars
table_vi.columns = table_vi_cars
In [38]:
table_vi
Out[38]:
MZ323 NISENT FDESCO CVCAVA HDACCO FDTAUR BKCENT NIMAXI ACLEGE LNTOWN CDSEVI LXLS40 BW735i
MZ323 -105.120 1.151 6.893 6.806 1.853 1.096 0.387 0.094 0.019 0.025 0.006 0.007 0.001
NISENT 0.535 -88.818 5.392 5.334 2.128 1.203 0.418 0.123 0.028 0.037 0.010 0.012 0.001
FDESCO 0.549 0.925 -85.786 5.490 2.169 1.217 0.419 0.119 0.026 0.038 0.009 0.011 0.001
CVCAVA 0.531 0.895 5.370 -84.471 2.225 1.414 0.451 0.115 0.024 0.046 0.009 0.011 0.001
HDACCO 0.102 0.253 1.501 1.574 -41.254 1.289 0.366 0.229 0.078 0.181 0.036 0.048 0.007
FDTAUR 0.080 0.190 1.122 1.332 1.716 -36.240 0.359 0.168 0.063 0.269 0.034 0.038 0.007
BKCENT 0.079 0.184 1.079 1.187 1.361 1.004 -38.302 0.380 0.146 0.278 0.062 0.070 0.009
NIMAXI 0.022 0.061 0.345 0.338 0.954 0.527 0.426 -22.780 0.215 0.345 0.095 0.120 0.015
ACLEGE 0.008 0.026 0.142 0.131 0.609 0.366 0.305 0.401 -15.909 0.377 0.106 0.134 0.022
LNTOWN 0.004 0.012 0.073 0.092 0.509 0.566 0.210 0.232 0.136 -13.393 0.101 0.129 0.027
CDSEVI 0.004 0.014 0.079 0.078 0.448 0.318 0.210 0.285 0.171 0.454 -12.468 0.138 0.028
LXLS40 0.004 0.014 0.078 0.076 0.486 0.292 0.190 0.293 0.174 0.466 0.111 -12.423 0.029
BW735i 0.002 0.006 0.035 0.038 0.264 0.221 0.102 0.146 0.114 0.381 0.090 0.115 -7.935

Citations

Berry, S., Levinsohn, J., & Pakes, A. (1995). Automobile prices in market equilibrium. Econometrica: Journal of the Econometric Society, 841-890.

Conlon, C., & Gortmaker, J. (2019). Best practices for differentiated products demand estimation with pyblp. Working paper. url: https://chrisconlon. github. io/site/pyblp. pdf.

Nevo, A. (2000). A practitioner's guide to estimation of random‐coefficients logit models of demand. Journal of economics & management strategy, 9(4), 513-548.