\n",
"\n",
"\n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
"

\n",
"

"
],
"text/plain": [
" market_ids city_ids quarter product_ids firm_ids brand_ids shares \\\n",
"0 C01Q1 1 1 F1B04 1 4 0.012417 \n",
"1 C01Q1 1 1 F1B06 1 6 0.007809 \n",
"2 C01Q1 1 1 F1B07 1 7 0.012995 \n",
"3 C01Q1 1 1 F1B09 1 9 0.005770 \n",
"4 C01Q1 1 1 F1B11 1 11 0.017934 \n",
"\n",
" prices sugar mushy ... demand_instruments10 \\\n",
"0 0.072088 2 1 ... 2.116358 \n",
"1 0.114178 18 1 ... -7.374091 \n",
"2 0.132391 4 1 ... 2.187872 \n",
"3 0.130344 3 0 ... 2.704576 \n",
"4 0.154823 12 0 ... 1.261242 \n",
"\n",
" demand_instruments11 demand_instruments12 demand_instruments13 \\\n",
"0 -0.154708 -0.005796 0.014538 \n",
"1 -0.576412 0.012991 0.076143 \n",
"2 -0.207346 0.003509 0.091781 \n",
"3 0.040748 -0.003724 0.094732 \n",
"4 0.034836 -0.000568 0.102451 \n",
"\n",
" demand_instruments14 demand_instruments15 demand_instruments16 \\\n",
"0 0.126244 0.067345 0.068423 \n",
"1 0.029736 0.087867 0.110501 \n",
"2 0.163773 0.111881 0.108226 \n",
"3 0.135274 0.088090 0.101767 \n",
"4 0.130640 0.084818 0.101075 \n",
"\n",
" demand_instruments17 demand_instruments18 demand_instruments19 \n",
"0 0.034800 0.126346 0.035484 \n",
"1 0.087784 0.049872 0.072579 \n",
"2 0.086439 0.122347 0.101842 \n",
"3 0.101777 0.110741 0.104332 \n",
"4 0.125169 0.133464 0.121111 \n",
"\n",
"[5 rows x 30 columns]"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"product_data = pd.read_csv(pyblp.data.NEVO_PRODUCTS_LOCATION)\n",
"product_data.head()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The product data contains `market_ids`, `product_ids`, `firm_ids`, `shares`, `prices`, a number of other firm IDs and product characteristics, and some pre-computed excluded `demand_instruments0`, `demand_instruments1`, and so on. The `product_ids` will be incorporated as fixed effects. \n",
"\n",
"For more information about the instruments and the example data as a whole, refer to the [`data`](https://pyblp.readthedocs.io/en/stable/_api/pyblp.data.html#module-pyblp.data) module."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Setting Up and Solving the Problem Without Demographics\n",
"\n",
"Formulations, product data, and an integration configuration are collectively used to initialize a [`Problem`](https://pyblp.readthedocs.io/en/stable/_api/pyblp.Problem.html#pyblp.Problem). Once initialized, [`Problem.solve`](https://pyblp.readthedocs.io/en/stable/_api/pyblp.Problem.solve.html#pyblp.Problem.solve) runs the estimation routine. The arguments to [`Problem.solve`](https://pyblp.readthedocs.io/en/stable/_api/pyblp.Problem.solve.html#pyblp.Problem.solve) configure how estimation is performed. For example, `optimization` and `iteration` arguments configure the optimization and iteration routines that are used by the outer and inner loops of estimation.\n",
"\n",
"We'll specify [`Formulation`](https://pyblp.readthedocs.io/en/stable/_api/pyblp.Formulation.html#pyblp.Formulation) configurations for $X_1$, the linear characteristics, and $X_2$, the nonlinear characteristics.\n",
"\n",
"- The formulation for $X_1$ consists of `prices` and fixed effects constructed from `product_ids`, which we will absorb using `absorb` argument of [`Formulation`](https://pyblp.readthedocs.io/en/stable/_api/pyblp.Formulation.html#pyblp.Formulation).\n",
"- If we were interested in reporting estimates for each fixed effect, we could replace the formulation for $X_1$ with `Formulation('prices + C(product_ids)')`.\n",
"- Because `sugar`, `mushy`, and the constant are collinear with `product_ids`, we can include them in $X_2$ but not in $X_1$."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(prices + Absorb[C(product_ids)], 1 + prices + sugar + mushy)"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"X1_formulation = pyblp.Formulation('0 + prices', absorb='C(product_ids)')\n",
"X2_formulation = pyblp.Formulation('1 + prices + sugar + mushy')\n",
"product_formulations = (X1_formulation, X2_formulation)\n",
"product_formulations"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We also need to specify an [`Integration`](https://pyblp.readthedocs.io/en/stable/_api/pyblp.Integration.html#pyblp.Integration) configuration. We consider two alternatives:\n",
"\n",
"1. Monte Carlo draws: we simulate 50 individuals from a random normal distribution.\n",
"2. Product rules: we construct nodes and weights according to a product rule that exactly integrates polynomials of degree $5 \\times 2 - 1 = 9$ or less."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Configured to construct nodes and weights with Monte Carlo simulation."
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"mc_integration = pyblp.Integration('monte_carlo', size=50, seed=0)\n",
"mc_integration"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Configured to construct nodes and weights according to the level-5 Gauss-Hermite product rule."
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pr_integration = pyblp.Integration('product', size=5)\n",
"pr_integration"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Dimensions:\n",
"============================================\n",
" T N F I K1 K2 MD ED \n",
"--- ---- --- ---- ---- ---- ---- ----\n",
"94 2256 5 4700 1 4 20 1 \n",
"============================================\n",
"\n",
"Formulations:\n",
"===========================================================\n",
" Column Indices: 0 1 2 3 \n",
"----------------------------- ------ ------ ----- -----\n",
" X1: Linear Characteristics prices \n",
"X2: Nonlinear Characteristics 1 prices sugar mushy\n",
"==========================================================="
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"mc_problem = pyblp.Problem(product_formulations, product_data, integration=mc_integration)\n",
"mc_problem"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Dimensions:\n",
"=============================================\n",
" T N F I K1 K2 MD ED \n",
"--- ---- --- ----- ---- ---- ---- ----\n",
"94 2256 5 58750 1 4 20 1 \n",
"=============================================\n",
"\n",
"Formulations:\n",
"===========================================================\n",
" Column Indices: 0 1 2 3 \n",
"----------------------------- ------ ------ ----- -----\n",
" X1: Linear Characteristics prices \n",
"X2: Nonlinear Characteristics 1 prices sugar mushy\n",
"==========================================================="
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pr_problem = pyblp.Problem(product_formulations, product_data, integration=pr_integration)\n",
"pr_problem"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As an illustration of how to configure the optimization routine, we'll use a simpler, non-default [`Optimization`](https://pyblp.readthedocs.io/en/stable/_api/pyblp.Optimization.html#pyblp.Optimization) configuration that doesn't support parameter bounds."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Configured to optimize using the BFGS algorithm implemented in SciPy with analytic gradients and options {}."
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"bfgs = pyblp.Optimization('bfgs')\n",
"bfgs"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We estimate three versions of the model:\n",
"\n",
"1. An unrestricted covariance matrix for random tastes using Monte Carlo integration.\n",
"2. An unrestricted covariance matrix for random tastes using the product rule.\n",
"3. A restricted diagonal matrix for random tastes using Monte Carlo Integration.\n",
"\n",
"Notice that the only thing that changes when we estimate the restricted covariance is our initial guess of $\\Sigma_0$. The lower diagonal in this initial guess is ignored because we are optimizing over the upper triangular (Cholesky) root of $\\Sigma$."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Problem Results Summary:\n",
"========================================================================================================================\n",
" Smallest Largest \n",
"Computation GMM Optimization Objective Fixed Point Contraction Objective Gradient Hessian Hessian \n",
" Time Step Iterations Evaluations Iterations Evaluations Value Infinity Norm Eigenvalue Eigenvalue\n",
"----------- ---- ------------ ----------- ----------- ----------- --------- ------------- ---------- ----------\n",
" 00:05:19 2 62 171 175814 541988 +3.0E+05 +8.0E-03 +1.5E+02 +1.2E+07 \n",
"========================================================================================================================\n",
"\n",
"Nonlinear Coefficient Estimates (Robust SEs in Parentheses):\n",
"======================================================\n",
"Sigma: 1 prices sugar mushy \n",
"------ ---------- ---------- ---------- ----------\n",
" 1 +9.1E-02 +1.3E-01 +1.1E-01 -5.7E-02 \n",
" (+1.8E+00) (+8.2E+00) (+9.5E+00) (+4.3E+00)\n",
" \n",
"prices -1.4E+00 +3.5E+00 +1.3E+01 \n",
" (+6.6E+01) (+7.4E+01) (+1.9E+01)\n",
" \n",
"sugar -1.1E-02 -1.4E-01 \n",
" (+1.7E-01) (+1.8E-01)\n",
" \n",
"mushy -2.6E-02 \n",
" (+7.5E-01)\n",
"======================================================\n",
"\n",
"Beta Estimates (Robust SEs in Parentheses):\n",
"==========\n",
" prices \n",
"----------\n",
" -3.4E+01 \n",
"(+1.5E+01)\n",
"=========="
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"results1 = mc_problem.solve(sigma=np.ones((4, 4)), optimization=bfgs)\n",
"results1"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Problem Results Summary:\n",
"========================================================================================================================\n",
" Smallest Largest \n",
"Computation GMM Optimization Objective Fixed Point Contraction Objective Gradient Hessian Hessian \n",
" Time Step Iterations Evaluations Iterations Evaluations Value Infinity Norm Eigenvalue Eigenvalue\n",
"----------- ---- ------------ ----------- ----------- ----------- --------- ------------- ---------- ----------\n",
" 00:10:49 2 69 153 100768 312910 +3.6E+05 +1.3E-02 +3.3E+00 +1.2E+07 \n",
"========================================================================================================================\n",
"\n",
"Nonlinear Coefficient Estimates (Robust SEs in Parentheses):\n",
"======================================================\n",
"Sigma: 1 prices sugar mushy \n",
"------ ---------- ---------- ---------- ----------\n",
" 1 -9.8E-10 -1.0E-09 +3.2E-01 +6.7E-01 \n",
" (+2.2E+09) (+5.1E+09) (+9.9E+01) (+5.2E+01)\n",
" \n",
"prices +1.4E-08 -6.2E+00 -1.1E+01 \n",
" (+2.8E+10) (+1.6E+03) (+9.3E+02)\n",
" \n",
"sugar +5.7E-02 +9.7E-02 \n",
" (+1.4E+01) (+7.9E+00)\n",
" \n",
"mushy -1.6E-01 \n",
" (+1.5E+01)\n",
"======================================================\n",
"\n",
"Beta Estimates (Robust SEs in Parentheses):\n",
"==========\n",
" prices \n",
"----------\n",
" -3.1E+01 \n",
"(+3.1E+01)\n",
"=========="
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"results2 = pr_problem.solve(sigma=np.ones((4, 4)), optimization=bfgs)\n",
"results2"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Problem Results Summary:\n",
"========================================================================================================================\n",
" Smallest Largest \n",
"Computation GMM Optimization Objective Fixed Point Contraction Objective Gradient Hessian Hessian \n",
" Time Step Iterations Evaluations Iterations Evaluations Value Infinity Norm Eigenvalue Eigenvalue\n",
"----------- ---- ------------ ----------- ----------- ----------- --------- ------------- ---------- ----------\n",
" 00:01:19 2 16 77 43049 134868 +4.0E+05 +3.9E-04 +2.4E+03 +1.4E+07 \n",
"========================================================================================================================\n",
"\n",
"Nonlinear Coefficient Estimates (Robust SEs in Parentheses):\n",
"======================================================\n",
"Sigma: 1 prices sugar mushy \n",
"------ ---------- ---------- ---------- ----------\n",
" 1 +5.2E-02 +0.0E+00 +0.0E+00 +0.0E+00 \n",
" (+1.1E+00) \n",
" \n",
"prices -4.3E-01 +0.0E+00 +0.0E+00 \n",
" (+8.0E+00) \n",
" \n",
"sugar +3.6E-02 +0.0E+00 \n",
" (+5.8E-02) \n",
" \n",
"mushy +5.0E-01 \n",
" (+1.7E+00)\n",
"======================================================\n",
"\n",
"Beta Estimates (Robust SEs in Parentheses):\n",
"==========\n",
" prices \n",
"----------\n",
" -3.0E+01 \n",
"(+1.4E+00)\n",
"=========="
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"results3 = mc_problem.solve(sigma=np.eye(4), optimization=bfgs)\n",
"results3"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We see that all three models give similar estimates of the price coefficient $\\hat{\\alpha} \\approx -30$. Note a few of the estimated terms on the diagonal of $\\Sigma$ are negative. Since the diagonal consists of standard deviations, negative values are unrealistic. When using another optimization routine that supports bounds (like the default L-BFGS-B routine), these diagonal elements are by default bounded from below by zero."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Adding Demographics to the Problem\n",
"\n",
"To add demographic data we need to make two changes:\n",
"\n",
"1. We need to load `agent_data`, which for this cereal problem contains pre-computed Monte Carlo draws and demographics.\n",
"2. We need to add an `agent_formulation` to the model.\n",
"\n",
"The `agent data` has several reserved column names.\n",
"\n",
"- `market_ids` are the index that link the `agent data` to the `market_ids` in `product data`.\n",
"- `weights` are the weights $w_{ti}$ attached to each agent. In each market, these should sum to one so that $\\sum_i w_{ti} = 1$. It is often the case the $w_{ti} = 1 / I_t$ where $I_t$ is the number of agents in market $t$, so that each agent gets equal weight. Other possibilities include quadrature nodes and weights.\n",
"- `nodes0`, `nodes1`, and so on are the nodes at which the unobserved agent tastes $\\mu_{jti}$ are evaluated. The nodes should be labeled from $0, \\ldots, K_2 - 1$ where $K_2$ is the number of random coefficients.\n",
"- Other fields are the realizations of the demographics $d$ themselves."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"market_ids | city_ids | quarter | product_ids | firm_ids | brand_ids | shares | prices | sugar | mushy | ... | demand_instruments10 | demand_instruments11 | demand_instruments12 | demand_instruments13 | demand_instruments14 | demand_instruments15 | demand_instruments16 | demand_instruments17 | demand_instruments18 | demand_instruments19 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|

0 | C01Q1 | 1 | 1 | F1B04 | 1 | 4 | 0.012417 | 0.072088 | 2 | 1 | ... | 2.116358 | -0.154708 | -0.005796 | 0.014538 | 0.126244 | 0.067345 | 0.068423 | 0.034800 | 0.126346 | 0.035484 |

1 | C01Q1 | 1 | 1 | F1B06 | 1 | 6 | 0.007809 | 0.114178 | 18 | 1 | ... | -7.374091 | -0.576412 | 0.012991 | 0.076143 | 0.029736 | 0.087867 | 0.110501 | 0.087784 | 0.049872 | 0.072579 |

2 | C01Q1 | 1 | 1 | F1B07 | 1 | 7 | 0.012995 | 0.132391 | 4 | 1 | ... | 2.187872 | -0.207346 | 0.003509 | 0.091781 | 0.163773 | 0.111881 | 0.108226 | 0.086439 | 0.122347 | 0.101842 |

3 | C01Q1 | 1 | 1 | F1B09 | 1 | 9 | 0.005770 | 0.130344 | 3 | 0 | ... | 2.704576 | 0.040748 | -0.003724 | 0.094732 | 0.135274 | 0.088090 | 0.101767 | 0.101777 | 0.110741 | 0.104332 |

4 | C01Q1 | 1 | 1 | F1B11 | 1 | 11 | 0.017934 | 0.154823 | 12 | 0 | ... | 1.261242 | 0.034836 | -0.000568 | 0.102451 | 0.130640 | 0.084818 | 0.101075 | 0.125169 | 0.133464 | 0.121111 |

5 rows \u00d7 30 columns

\n", "\n",
"\n",
"\n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
"

\n",
"

"
],
"text/plain": [
" market_ids city_ids quarter weights nodes0 nodes1 nodes2 \\\n",
"0 C01Q1 1 1 0.05 0.434101 -1.500838 -1.151079 \n",
"1 C01Q1 1 1 0.05 -0.726649 0.133182 -0.500750 \n",
"2 C01Q1 1 1 0.05 -0.623061 -0.138241 0.797441 \n",
"3 C01Q1 1 1 0.05 -0.041317 1.257136 -0.683054 \n",
"4 C01Q1 1 1 0.05 -0.466691 0.226968 1.044424 \n",
"\n",
" nodes3 income income_squared age child \n",
"0 0.161017 0.495123 8.331304 -0.230109 -0.230851 \n",
"1 0.129732 0.378762 6.121865 -2.532694 0.769149 \n",
"2 -0.795549 0.105015 1.030803 -0.006965 -0.230851 \n",
"3 0.259044 -1.485481 -25.583605 -0.827946 0.769149 \n",
"4 0.092019 -0.316597 -6.517009 -0.230109 -0.230851 "
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"agent_data = pd.read_csv(pyblp.data.NEVO_AGENTS_LOCATION)\n",
"agent_data.head()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The `agent formulation` tells us which columns of demographic information to interact with $X_2$."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"income + income_squared + age + child"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"agent_formulation = pyblp.Formulation('0 + income + income_squared + age + child')\n",
"agent_formulation"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This tells us to include demographic realizations for `income`, `income_squared`, `age`, and the presence of children, `child`, but to ignore other possible demographics when interacting demographics with $X_2$. We should also generally exclude the constant from the demographic formula.\n",
"\n",
"Now we configure the [`Problem`](https://pyblp.readthedocs.io/en/stable/_api/pyblp.Problem.html#pyblp.Problem) to include the `agent_formulation` and `agent_data`, which follow the `product_formulations` and `product_data`.\n",
"\n",
"When we display the class, it lists the demographic interactions in the table of formulations and reports $D = 4$, the dimension of the demographic draws."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Dimensions:\n",
"=================================================\n",
" T N F I K1 K2 D MD ED \n",
"--- ---- --- ---- ---- ---- --- ---- ----\n",
"94 2256 5 1880 1 4 4 20 1 \n",
"=================================================\n",
"\n",
"Formulations:\n",
"===================================================================\n",
" Column Indices: 0 1 2 3 \n",
"----------------------------- ------ -------------- ----- -----\n",
" X1: Linear Characteristics prices \n",
"X2: Nonlinear Characteristics 1 prices sugar mushy\n",
" d: Demographics income income_squared age child\n",
"==================================================================="
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"nevo_problem = pyblp.Problem(\n",
" product_formulations,\n",
" product_data,\n",
" agent_formulation,\n",
" agent_data\n",
")\n",
"nevo_problem"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The initialized problem can be solved with [`Problem.solve`](https://pyblp.readthedocs.io/en/stable/_api/pyblp.Problem.solve.html#pyblp.Problem.solve). We'll use the same starting values as [Nevo (2000)](https://pyblp.readthedocs.io/en/stable/references.html#nevo-2000). By passing a diagonal matrix as starting values for $\\Sigma$, we're choosing to ignore covariance terms. Similarly, zeros in the starting values for $\\Pi$ mean that those parameters will be fixed at zero.\n",
"\n",
"To replicate common estimates, we'll use the non-default BFGS optimization routine, and we'll configure [`Problem.solve`](https://pyblp.readthedocs.io/en/stable/_api/pyblp.Problem.solve.html#pyblp.Problem.solve) to use 1-step GMM instead of 2-step GMM."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Problem Results Summary:\n",
"========================================================================================================================\n",
" Smallest Largest \n",
"Computation GMM Optimization Objective Fixed Point Contraction Objective Gradient Hessian Hessian \n",
" Time Step Iterations Evaluations Iterations Evaluations Value Infinity Norm Eigenvalue Eigenvalue\n",
"----------- ---- ------------ ----------- ----------- ----------- --------- ------------- ---------- ----------\n",
" 00:01:18 1 51 58 46236 143503 +4.6E+00 +6.9E-06 +2.8E-05 +1.6E+04 \n",
"========================================================================================================================\n",
"\n",
"Nonlinear Coefficient Estimates (Robust SEs in Parentheses):\n",
"=====================================================================================================================\n",
"Sigma: 1 prices sugar mushy | Pi: income income_squared age child \n",
"------ ---------- ---------- ---------- ---------- | ------ ---------- -------------- ---------- ----------\n",
" 1 +5.6E-01 +0.0E+00 +0.0E+00 +0.0E+00 | 1 +2.3E+00 +0.0E+00 +1.3E+00 +0.0E+00 \n",
" (+1.6E-01) | (+1.2E+00) (+6.3E-01) \n",
" | \n",
"prices +3.3E+00 +0.0E+00 +0.0E+00 | prices +5.9E+02 -3.0E+01 +0.0E+00 +1.1E+01 \n",
" (+1.3E+00) | (+2.7E+02) (+1.4E+01) (+4.1E+00)\n",
" | \n",
"sugar -5.8E-03 +0.0E+00 | sugar -3.8E-01 +0.0E+00 +5.2E-02 +0.0E+00 \n",
" (+1.4E-02) | (+1.2E-01) (+2.6E-02) \n",
" | \n",
"mushy +9.3E-02 | mushy +7.5E-01 +0.0E+00 -1.4E+00 +0.0E+00 \n",
" (+1.9E-01) | (+8.0E-01) (+6.7E-01) \n",
"=====================================================================================================================\n",
"\n",
"Beta Estimates (Robust SEs in Parentheses):\n",
"==========\n",
" prices \n",
"----------\n",
" -6.3E+01 \n",
"(+1.5E+01)\n",
"=========="
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"initial_sigma = np.diag([0.3302, 2.4526, 0.0163, 0.2441])\n",
"initial_pi = np.array([\n",
" [ 5.4819, 0, 0.2037, 0 ],\n",
" [15.8935, -1.2000, 0, 2.6342],\n",
" [-0.2506, 0, 0.0511, 0 ],\n",
" [ 1.2650, 0, -0.8091, 0 ]\n",
"])\n",
"nevo_results = nevo_problem.solve(\n",
" initial_sigma,\n",
" initial_pi,\n",
" optimization=bfgs,\n",
" method='1s'\n",
")\n",
"nevo_results"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Results are similar to those in the original paper with a GMM objective value of $q(\\hat{\\theta}) = 4.56$ and a price coefficient of $\\hat{\\alpha} = -62.7$. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Restricting Parameters\n",
"\n",
"Because the interactions between `price`, `income`, and `income_squared` are potentially collinear, we might worry that $\\hat{\\Pi}_{21} = 588$ and $\\hat{\\Pi}_{22} = -30.2$ are pushing the price coefficient in opposite directions. Both are large in magnitude but statistically insignficant. One way of dealing with this is to restrict $\\Pi_{22} = 0$.\n",
"\n",
"There are two ways we can do this:\n",
"\n",
"1. Change the initial $\\Pi_0$ values to make this term zero.\n",
"2. Change the agent formula to drop `income_squared`.\n",
"\n",
"First, we'll change the initial $\\Pi_0$ values."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Problem Results Summary:\n",
"========================================================================================================================\n",
" Smallest Largest \n",
"Computation GMM Optimization Objective Fixed Point Contraction Objective Gradient Hessian Hessian \n",
" Time Step Iterations Evaluations Iterations Evaluations Value Infinity Norm Eigenvalue Eigenvalue\n",
"----------- ---- ------------ ----------- ----------- ----------- --------- ------------- ---------- ----------\n",
" 00:00:57 1 34 41 34457 106706 +1.5E+01 +5.2E-06 +4.7E-02 +1.7E+04 \n",
"========================================================================================================================\n",
"\n",
"Nonlinear Coefficient Estimates (Robust SEs in Parentheses):\n",
"=====================================================================================================================\n",
"Sigma: 1 prices sugar mushy | Pi: income income_squared age child \n",
"------ ---------- ---------- ---------- ---------- | ------ ---------- -------------- ---------- ----------\n",
" 1 +3.7E-01 +0.0E+00 +0.0E+00 +0.0E+00 | 1 +3.1E+00 +0.0E+00 +1.2E+00 +0.0E+00 \n",
" (+1.2E-01) | (+1.1E+00) (+1.0E+00) \n",
" | \n",
"prices +1.8E+00 +0.0E+00 +0.0E+00 | prices +4.2E+00 +0.0E+00 +0.0E+00 +1.2E+01 \n",
" (+9.2E-01) | (+4.6E+00) (+5.2E+00)\n",
" | \n",
"sugar -4.4E-03 +0.0E+00 | sugar -1.9E-01 +0.0E+00 +2.8E-02 +0.0E+00 \n",
" (+1.2E-02) | (+3.5E-02) (+3.2E-02) \n",
" | \n",
"mushy +8.6E-02 | mushy +1.5E+00 +0.0E+00 -1.5E+00 +0.0E+00 \n",
" (+1.9E-01) | (+6.5E-01) (+1.1E+00) \n",
"=====================================================================================================================\n",
"\n",
"Beta Estimates (Robust SEs in Parentheses):\n",
"==========\n",
" prices \n",
"----------\n",
" -3.2E+01 \n",
"(+2.3E+00)\n",
"=========="
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"restricted_pi = initial_pi.copy()\n",
"restricted_pi[1, 1] = 0\n",
"nevo_problem.solve(\n",
" initial_sigma,\n",
" restricted_pi,\n",
" optimization=bfgs,\n",
" method='1s'\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we'll drop both `income_squared` and the corresponding column in $\\Pi_0$."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Problem Results Summary:\n",
"========================================================================================================================\n",
" Smallest Largest \n",
"Computation GMM Optimization Objective Fixed Point Contraction Objective Gradient Hessian Hessian \n",
" Time Step Iterations Evaluations Iterations Evaluations Value Infinity Norm Eigenvalue Eigenvalue\n",
"----------- ---- ------------ ----------- ----------- ----------- --------- ------------- ---------- ----------\n",
" 00:00:54 1 34 41 34472 106726 +1.5E+01 +5.2E-06 +4.7E-02 +1.7E+04 \n",
"========================================================================================================================\n",
"\n",
"Nonlinear Coefficient Estimates (Robust SEs in Parentheses):\n",
"=====================================================================================================\n",
"Sigma: 1 prices sugar mushy | Pi: income age child \n",
"------ ---------- ---------- ---------- ---------- | ------ ---------- ---------- ----------\n",
" 1 +3.7E-01 +0.0E+00 +0.0E+00 +0.0E+00 | 1 +3.1E+00 +1.2E+00 +0.0E+00 \n",
" (+1.2E-01) | (+1.1E+00) (+1.0E+00) \n",
" | \n",
"prices +1.8E+00 +0.0E+00 +0.0E+00 | prices +4.2E+00 +0.0E+00 +1.2E+01 \n",
" (+9.2E-01) | (+4.6E+00) (+5.2E+00)\n",
" | \n",
"sugar -4.4E-03 +0.0E+00 | sugar -1.9E-01 +2.8E-02 +0.0E+00 \n",
" (+1.2E-02) | (+3.5E-02) (+3.2E-02) \n",
" | \n",
"mushy +8.6E-02 | mushy +1.5E+00 -1.5E+00 +0.0E+00 \n",
" (+1.9E-01) | (+6.5E-01) (+1.1E+00) \n",
"=====================================================================================================\n",
"\n",
"Beta Estimates (Robust SEs in Parentheses):\n",
"==========\n",
" prices \n",
"----------\n",
" -3.2E+01 \n",
"(+2.3E+00)\n",
"=========="
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"restricted_formulation = pyblp.Formulation('0 + income + age + child')\n",
"deleted_pi = np.delete(initial_pi, 1, axis=1)\n",
"restricted_problem = pyblp.Problem(\n",
" product_formulations,\n",
" product_data,\n",
" restricted_formulation,\n",
" agent_data\n",
")\n",
"restricted_problem.solve(\n",
" initial_sigma,\n",
" deleted_pi,\n",
" optimization=bfgs,\n",
" method='1s'\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The parameter estimates and standard errors are identical for both approaches. Based on the number of fixed point iterations, there is some evidence that the solver took a slightly different path for each problem, but both restricted problems arrived at identical answers. The $\\hat{\\Pi}_{12}$ interaction term is still insignificant."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.6"
}
},
"nbformat": 4,
"nbformat_minor": 2
}market_ids | city_ids | quarter | weights | nodes0 | nodes1 | nodes2 | nodes3 | income | income_squared | age | child | |
---|---|---|---|---|---|---|---|---|---|---|---|---|

0 | C01Q1 | 1 | 1 | 0.05 | 0.434101 | -1.500838 | -1.151079 | 0.161017 | 0.495123 | 8.331304 | -0.230109 | -0.230851 |

1 | C01Q1 | 1 | 1 | 0.05 | -0.726649 | 0.133182 | -0.500750 | 0.129732 | 0.378762 | 6.121865 | -2.532694 | 0.769149 |

2 | C01Q1 | 1 | 1 | 0.05 | -0.623061 | -0.138241 | 0.797441 | -0.795549 | 0.105015 | 1.030803 | -0.006965 | -0.230851 |

3 | C01Q1 | 1 | 1 | 0.05 | -0.041317 | 1.257136 | -0.683054 | 0.259044 | -1.485481 | -25.583605 | -0.827946 | 0.769149 |

4 | C01Q1 | 1 | 1 | 0.05 | -0.466691 | 0.226968 | 1.044424 | 0.092019 | -0.316597 | -6.517009 | -0.230109 | -0.230851 |