Lecture 5: Advanced Econometrics BLP and Rust ¶
Antoine Chapel (Sciences Po & PSE) ¶
Alfred Galichon's math+econ+code prerequisite class on numerical optimization and econometrics, in Python ¶
Class content by Antoine Chapel. Past and present support from Alfred Galichon's ERC grant CoG-866274 is acknowledged, as well as inputs from contributors listed here. If you reuse material from this class, please cite as:
Antoine Chapel, 'math+econ+code' prerequisite class on numerical optimization and econometrics, January 2023
In this last session, we will study two advanced econometrics topics that will be studied in the M+E+C masterclass. The masterclass is assuming knowledge of the principles for these two techniques, and will re-visit them using Optimal Transport. Today, we will build the basics for these two tools and introduce two fields of economics in the process:
- B(erry)-L(evinsohn)-P(akes): a structural model in industrial organization
- Rust: a structural model for dynamic discrete choice
These two objects are kept for the end because you now have all the intermediate tools that are necessary to understand them, and the mainstream algorithm that can be used for both of them has a similar structure called: nested fixed point (NFXP)
References¶
- "Automobile Prices in Market Equilibrium", (Berry, Levinsohn and Pakes), Econometrica Vol. 63, No. 4 (Jul., 1995), pp. 841-890
- "Best practices for differentiated products demand estimation with PyBLP", (Conlon and Gortmaker), The RAND Journal of Economics, volume 51, No. 4, 2020
- Lecture notes on Empirical IO (Marleen Marra)
- "Optimal Replacement of GMC Bus Engines: An Empirical Model of Harold Zurcher", (John Rust), Econometrica, Vol. 55, No. 5 (Sep., 1987), pp. 999-1033
- Nested Fixed Point Algorithm Documentation Manual (Rust), version 6, 2000
The data used here and the code for Rust was taken and adapted from an excellent Rust replication by Natasha Watkins: we invite you to take a look at the original notebook for a more complete replication of the Rust paper and a good tutorial on how to handle the data from the original Rust paper: https://notes.quantecon.org/submission/6234fe0f96e1ce001b61fad8
BLP: Random Coefficients Logit Model¶
We will study the simplest possible form of BLP and ignore many of the subtleties of the topic. Great papers (BLP 1995, Nevo 2000) have been written to explain the random coefficients logit model and we refer you to that for the details. The model also has a lot of variations, so what we present here is focused on the basic of idea of what BLP is and how the algorithm works.
Utility function is specified as follows, where $i$ indexes individuals, and $j$ indexes alternatives:
\begin{align} u_{ij} = x_j'\beta_i - \alpha p_j + \xi_j + \epsilon_{ij} \end{align}
- $x_j$ is a vector of characteristics specific to alternative $j$
- $p_j$ is the price of alternative $j$
- $\xi_j$ is the unobserved heterogeneity specific to characteristic $j$
- $\epsilon_{ij}$ is the econometric error of the model, which we assumed to be Gumbel distributed.
How do these differ from the models you are used to ? The coefficient $\beta$ is indexed per individual. Denoting $\overline{\beta}$ the mean of this "random coefficient", and considering that the vector $x_j$ contains $k$ characteristics:
\begin{align} x_j' \beta_i = x_j \overline{\beta} + \sum_{k=1}^K \sigma_k x_{jk} \nu_{ik} \end{align}
- $\nu_{ik} \sim \mathcal{N}(0, 1)$ is an individual-characteristic idiosyncratic shock. This is what drives the randomness in $\beta_i$
- $\overline{\beta}_k$ and $\sigma_k$ are the parameters to estimate
\begin{align*} u_{ij} &= x_j \overline{\beta} - \alpha p_j + \xi_j + \sum_{k=1}^K \sigma_k x_{jk} \nu_{ik} + \epsilon_{ij} \\ &= \delta_j + \sum_{k=1}^K \sigma_k x_{jk} \nu_{ik} + \epsilon_{ij}\\ &= \delta_j + \mu_{ij} \end{align*}
Where $\delta_j$ denotes the "mean utility" that alternative $j$ provides to individuals, and $\mu_{ij}$ denotes the deviation from that mean utility at the individual level.
For a given individual $i$, we can use the formula above to rewrite the choice probability of that individual using the multinomial logit form you know already:
\begin{align*} P_{ij|\nu_i} &= \frac{\exp(\delta_j + \sum_{k=1}^K \sigma_k x_{jk} \nu_{ik})}{1 + \sum_{h=1}^J \exp(\delta_h + \sum_{k=1}^K \sigma_k x_{hk} \nu_{ik})} \end{align*}
The only uncertain component in this expression is $\nu_{i}$, which varies across individuals. Therefore, the "market share" of alternative $j$, when integrated across the population, can be written as follows:
\begin{align} s_j = \int P_{ij|\nu_i} f_\nu(\nu_i) d\nu_i \end{align}
This integral is not necessarily easy to compute exactly: in their paper, BLP rely on simulation to compute it. That is, they will draw $N$ $v_i$ assuming they follow a normal distribution, and approximate the integral thereby, numerically.
One last element that's required for BLP: instruments. Indeed, prices are endogenous in this model, which is why we need to rely on instrumental variables $Z$. If you do not know anything about IV, check it in "Mostly Harmless Econometrics" or "Microeconometrics". As a reminder, instruments need to have two properties: they need to be strongly correlated with the endogenous variable (here, the price), and be uncorrelated with the variable we attempt to explain (here, $u_{ij}$).
The algorithm is composed of an outer loop and an inner loop. The outer loop is an IV GMM procedure, by which we try to find the parameters such that $E[\xi Z] = 0$ (moment condition for IV). Translated as a GMM problem, it seeks to determine the parameter $\theta$ that minimizes the (squared, weighted) empirical moment $\frac{1}{JT} \sum_{j, t} \xi_{jt}(\theta)'Z_{jt}$. The minimization process is the outer loop. We build the evaluation of the empirical moment in an inner loop, which follows the following steps.
Given a candidate $\theta_0 = (\beta_0, \sigma_0)$¶
- Set $\hat{\theta} = \theta_0$
- Using $\hat{\theta}$, determine the $\hat{\delta}_{jt}$ such that simulated market shares $\tilde{s}_j$ match observed market shares $s_j$. This is done through a contraction mapping that can be summarized as $\delta^{k+1} = \delta^k + \log(s_j) - \log(\tilde{s}_j)$
- Through linear regression, recover the unobserved heterogeneity that corresponds to the current guess of $\theta$: $\hat{\xi}(\hat{\theta}) = \hat{\delta}_{jt} + \alpha p_{jt} - x_{jt} \hat{\overline{\beta}}$
- Evaluate $\hat{E}[\hat{\xi}(\hat{\theta})'Z] = \frac{\sum_{j, t} z'_{jt} \hat{\xi}_{jt}(\theta)}{JT}$
- Return GMM objective $\big(\frac{\sum_{j, t} z'_{jt} \hat{\xi}_{jt}(\theta)}{JT}\big)'W\big(\frac{\sum_{j, t} z'_{jt} \hat{\xi}_{jt}(\theta)}{JT}\big)$
This is all fairly technical already, and we have ignored several complexities, in particular the contraction mapping that is used in BLP to do step 1. Fortunately, pyblp can do all of it on our behalf. The following crash tutorial is taken from the simplest version of BLP that is presented in the pyblp documentation, available here
#!pip install pyblp
#To run the pyblp part, you may need to revert your numpy installation to a former version (1.18.1 works).
import pandas as pd
import pyblp
import numpy as np
import matplotlib.pyplot as plt
#this (fake) cereal data is provided directly by the pyblp library to try out their code
product_data = pd.read_csv(pyblp.data.NEVO_PRODUCTS_LOCATION)
product_data.columns
Index(['market_ids', 'city_ids', 'quarter', 'product_ids', 'firm_ids', 'brand_ids', 'shares', 'prices', 'sugar', 'mushy', 'demand_instruments0', 'demand_instruments1', 'demand_instruments2', 'demand_instruments3', 'demand_instruments4', 'demand_instruments5', 'demand_instruments6', 'demand_instruments7', 'demand_instruments8', 'demand_instruments9', 'demand_instruments10', 'demand_instruments11', 'demand_instruments12', 'demand_instruments13', 'demand_instruments14', 'demand_instruments15', 'demand_instruments16', 'demand_instruments17', 'demand_instruments18', 'demand_instruments19'], dtype='object')
A short tutorial in pyblp¶
pyblp is a good tool to master if you are interested in industrial organization (IO) or if you want to RA for a professor who does IO. The official tutorial, from which the following was taken, is a great resource and we invite you to take a look at it if you want to go further than what is presented here.
- Load the data as a Pandas object (done already)
- Write a $X_1$ formulation which contains all variables for which you do not want random coefficients, i.e price in our example (in the example below, there is no constant $(0)$, but we add product fixed effect
- Write a $X_2$ formulation, which includes every variable for which you want random coefficients. $(1)$ is for the Constant
- Select a method for simulating the market shares
- write a pyblp.Problem object which aggregates the regression formulation, data and integration method.
#X1: Linear variables, no random coefficients
#X2: Nonlinear variables, random coefficients
X1_formulation = pyblp.Formulation('0 + prices', absorb='C(product_ids)')
X2_formulation = pyblp.Formulation('1 + prices + sugar + mushy')
product_formulations = (X1_formulation, X2_formulation)
#Numerical integration of the probability P_{ij}
mc_integration = pyblp.Integration('monte_carlo', size=50, specification_options={'seed': 0})
mc_integration
mc_problem = pyblp.Problem(product_formulations, product_data, integration=mc_integration)
mc_problem
Initializing the problem ... Absorbing demand-side fixed effects ... Initialized the problem after 00:00:00. Dimensions: ============================================ T N F I K1 K2 MD ED --- ---- --- ---- ---- ---- ---- ---- 94 2256 5 4700 1 4 20 1 ============================================ Formulations: =========================================================== Column Indices: 0 1 2 3 ----------------------------- ------ ------ ----- ----- X1: Linear Characteristics prices X2: Nonlinear Characteristics 1 prices sugar mushy ===========================================================
Dimensions: ============================================ T N F I K1 K2 MD ED --- ---- --- ---- ---- ---- ---- ---- 94 2256 5 4700 1 4 20 1 ============================================ Formulations: =========================================================== Column Indices: 0 1 2 3 ----------------------------- ------ ------ ----- ----- X1: Linear Characteristics prices X2: Nonlinear Characteristics 1 prices sugar mushy ===========================================================
pyblp.options.verbose = False
results = mc_problem.solve(sigma=np.ones((4, 4)), optimization=pyblp.Optimization('bfgs', {'gtol': 1e-4}))
results
Problem Results Summary: ================================================================================================================ GMM Objective Gradient Hessian Hessian Clipped Weighting Matrix Covariance Matrix Step Value Norm Min Eigenvalue Max Eigenvalue Shares Condition Number Condition Number ---- ------------- ------------- -------------- -------------- ------- ---------------- ----------------- 2 +1.483665E+02 +8.703747E-05 +8.515620E-02 +6.535576E+03 0 +5.150953E+07 +8.252073E+05 ================================================================================================================ Cumulative Statistics: =========================================================================== Computation Optimizer Optimization Objective Fixed Point Contraction Time Converged Iterations Evaluations Iterations Evaluations ----------- --------- ------------ ----------- ----------- ----------- 00:01:00 Yes 58 75 88300 271026 =========================================================================== Nonlinear Coefficient Estimates (Robust SEs in Parentheses): ================================================================================================================================================================= Sigma: 1 prices sugar mushy | Sigma Squared: 1 prices sugar mushy ------ --------------- --------------- --------------- --------------- | -------------- --------------- --------------- --------------- --------------- 1 +1.207566E+00 | 1 +1.458216E+00 -1.382374E+01 +7.314899E-02 -7.099420E-01 (+2.961805E+00) | (+7.153152E+00) (+5.188636E+01) (+2.222540E-01) (+2.272004E+00) | prices -1.144760E+01 +8.423600E+00 | prices -1.382374E+01 +2.020046E+02 -1.463091E+00 +1.492144E+00 (+1.774956E+01) (+1.157303E+01) | (+5.188636E+01) (+3.052988E+02) (+1.203963E+00) (+1.508886E+01) | sugar +6.057555E-02 -9.136790E-02 +3.783374E-02 | sugar +7.314899E-02 -1.463091E+00 +1.344888E-02 +2.034639E-02 (+2.485504E-01) (+2.282689E-01) (+8.298206E-02) | (+2.222540E-01) (+1.203963E+00) (+2.772954E-02) (+2.731843E-01) | mushy -5.879114E-01 -6.218281E-01 -2.261704E-02 +4.800048E-01 | mushy -7.099420E-01 +1.492144E+00 +2.034639E-02 +9.632262E-01 (+2.132815E+00) (+1.532679E+00) (+2.530561E+00) (+1.330829E+00) | (+2.272004E+00) (+1.508886E+01) (+2.731843E-01) (+3.964838E+00) ================================================================================================================================================================= Beta Estimates (Robust SEs in Parentheses): =============== prices --------------- -3.137350E+01 (+6.006485E+00) ===============
Rust 1987¶
In 1987, a paper by Rust introduced in a concise and understandable manner the estimation of dynamic discrete choice processes. So far, we have done static Discrete Choice, where an individual's choice is completely described by the alternatives available. Here, every choice has consequences on the future. Therefore, modelling of the individual's behaviour will combine dynamic programming and discrete choice models. Fortunately, knowing both tools, you are perfectly able to estimate the structural parameters in a dynamic discrete choice problem.
Dynamic Discrete Choice¶
In a quite general form, the dynamic discrete choice problem can be written this way. We focus here on an infinite horizon case, so the policy will be stationary. We denote present state $x$ and next period state $x'$.
\begin{align} V_\theta(x,\epsilon, d) &= \max_{d \in D(x)} \Big[ u(x, d, \theta) + \epsilon(d) + \beta\int_{x'} \int_{\epsilon'} V_\theta(x, \epsilon) \cdot π(x',\epsilon'|x, \epsilon, \theta) dx' d\epsilon' \Big]\\ V_\theta(x,\epsilon, d) &= \max_{d \in D(x)} \Big[ u(x, d, \theta) + \epsilon(d) + \beta E[V_\theta(x', \epsilon') ] \Big] \end{align}
Let's decompose this expression: when the agent is at state $x$, he faces a set of possible actions: $D(x)$. Given his decision, he will get a flow of immediate utility $u(x, d, \theta)$, some random shock $\epsilon(d)$, and a discounted value for future state $\beta E[V_\theta(x, \epsilon, d) ]$.
As you see there is some randomness in the value of future state. We integrate two times for the two sources of randomness: taking action $d$ does not guarantee us to get some state $x'$, and we ignore the value of the future random shock $\epsilon'$. So far this is a lot of uncertainty, and we can hardly go further. Thus we make two assumptions:
- $\epsilon \sim_{iid} \text{Gumbel}(0, 1)$
- $\pi(x'\epsilon' | x, e, \theta) = p(x'|x, d, \theta) \cdot q(\epsilon'|x, \theta)$
The first one is straightforward and implies as usual that we can rewrite conditional choice probabilities in closed form, but what does the second one mean ? We can separate the state transition from the random shock.
The Emax operator¶
Small parenthesis that will be useful anytime you tackle this type of problem. Let there be $j$ alternatives that provide systematic utility $U_j$ and a random shock $\epsilon_j$:
\begin{align} E\big[ \max \{U_j + \epsilon_j\} \big] = \gamma + \ln \sum_j e^{U_j} \end{align}
Where $\gamma = E[\epsilon] \approx 0.57$, the Euler-Mascheroni constant.
So:
\begin{align} &E\big[V_\theta (x)\big]\\ = &E\Big[ \max_{d∈D(x)} u(x, d, \theta) + \epsilon + \beta E[V_\theta(x')]\Big] \\ = &\int_{x} \Big[ \gamma + \log \big[ \sum_{d \in D(x)} \exp\{ u(x, d, \theta) + \beta E[V_\theta(x')\} \big]\Big] p(x|d, \theta) dx \end{align}
So, we integrate the $\epsilon$ uncertainty and rewrite the value function:
\begin{align} V_\theta(x) &= \max_{d∈D(x)} \Big[ u(x, d, \theta) + \epsilon + \beta E[V_\theta(x')\Big]\\ &= \\ &\max_{d∈D(x)} \Big[ u(x, d, \theta) + \epsilon + \beta\int_{x'} \Big[ \gamma +\log \big[ \sum_{d \in D(x')} \exp\{ u(x', d, \theta) + \beta E[V_\theta(x'', d)] \big] \Big] p(x'|d, \theta) dx' \end{align}
Remember our session on multinomial logit: we can use McFadden's formula: if alternative $j$ provides utility $U_{ij} = V_{ij} + \epsilon_{ij}$ to consumer $i$, the probability that he chooses alternative $j$ can be written:
\begin{align} P_{ij} = P(j = \arg \max U_{ij}) = \frac{e^{V_{ij}}}{\sum_k e^{V_{ik}}} \end{align}
Here, the equivalent of $U_{ij}$ in a dynamic setting is $V_\theta(x, d)$. Therefore, the Conditional Choice Probabilities are written as follows:
\begin{align} P(d|x, \theta) = \frac{\exp\{{u(x, d, \theta) + \beta E[V_\theta(x', d)]}\}}{\sum_{k \in D(x)} \exp\{{u(x', k, \theta) + \beta E[ V_\theta(x', k)] }\}} \end{align}
Let us now take an econometrics perspective: we can observe $(x, d)_{\{i, t\}}$, but not $\theta$. The estimation of $\hat{\theta}$ can be done through maximum likelihood:
\begin{align} LL(\theta) = \sum_t \sum_i \sum_j d_{ijt} \log \big( P(d_{ijt}|x_j, \theta) \big) \end{align}
Where $d_{ijt}$ is an indicator equal to $1$ if agent $i$ chose $j$ at time $t$. We can then find the ML estimate $\hat{\theta} = \arg \max_\theta LL(\theta)$
The Rust model (optimal bus engine replacement)¶
An agent is handling a fleet of buses. Given one bus, at every period, the agent can either replace the engine and incur a cost $RC$, or perform day-to-day reparations $c(x, \theta_1)$, where $x$ is the mileage of the bus. When the mileage of a bus increase, the cost of performing day-to-day reparations increases as well, so at some point, it is cost-efficient to just replace the engine, and reinitialize the mileage.
We denote the agent's utility as:
\begin{align} u(x, d, \theta) = \begin{cases} -RC(\theta) - c(0, \theta) + \epsilon & \text{if $d=1$ (replace)}\\ -c(x, \theta) + \epsilon & \text{if $d=0$ (not replace)} \end{cases} \end{align}
Mileage $x$ of a bus is reinitialized if $d=1$, and increases if $d=0$. Thus, we can write the transition probabilities as:
\begin{align} p(x'|x, d, \theta) = \begin{cases} g(x'- 0|\theta) &\text{if $d=1$}\\ g(x' - x|\theta) &\text{if $d=0$} \end{cases} \end{align}
from scipy.optimize import minimize
P_df = pd.read_csv(r"https://raw.githubusercontent.com/AntoineChapel/pedagogical_contents/main/rust/transition_matrix.csv").iloc[0:, 1:]
rust_data = pd.read_csv(r"https://raw.githubusercontent.com/AntoineChapel/pedagogical_contents/main/rust/rust_data.csv").iloc[:, 2:]
#There are ways to estimate it from the data, but to keep things accessible
#I am giving it here:
#How to interpret: every row/column index corresponds to a mileage category. 0: [0 5000]
# 1: [5000 10000], 2: [10000 15000]...
#If you are at mileage category 1, you have 35% chances to stay in that category, 64% chances to move above by one category,
# and 1% chances to move up by 2 categories. You can guess how these can be estimated from the data.
P_df
0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | ... | 80 | 81 | 82 | 83 | 84 | 85 | 86 | 87 | 88 | 89 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 0.35103 | 0.637445 | 0.011525 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.0 | 0.0 | 0.0 | ... | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.00000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 |
1 | 0.00000 | 0.351030 | 0.637445 | 0.011525 | 0.000000 | 0.000000 | 0.000000 | 0.0 | 0.0 | 0.0 | ... | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.00000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 |
2 | 0.00000 | 0.000000 | 0.351030 | 0.637445 | 0.011525 | 0.000000 | 0.000000 | 0.0 | 0.0 | 0.0 | ... | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.00000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 |
3 | 0.00000 | 0.000000 | 0.000000 | 0.351030 | 0.637445 | 0.011525 | 0.000000 | 0.0 | 0.0 | 0.0 | ... | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.00000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 |
4 | 0.00000 | 0.000000 | 0.000000 | 0.000000 | 0.351030 | 0.637445 | 0.011525 | 0.0 | 0.0 | 0.0 | ... | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.00000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 |
... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
85 | 0.00000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.0 | 0.0 | 0.0 | ... | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.35103 | 0.637445 | 0.011525 | 0.000000 | 0.000000 |
86 | 0.00000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.0 | 0.0 | 0.0 | ... | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.00000 | 0.351030 | 0.637445 | 0.011525 | 0.000000 |
87 | 0.00000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.0 | 0.0 | 0.0 | ... | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.00000 | 0.000000 | 0.351030 | 0.637445 | 0.011525 |
88 | 0.00000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.0 | 0.0 | 0.0 | ... | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.00000 | 0.000000 | 0.000000 | 0.351030 | 0.648970 |
89 | 0.00000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.0 | 0.0 | 0.0 | ... | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.00000 | 0.000000 | 0.000000 | 0.000000 | 1.000000 |
90 rows × 90 columns
P = P_df.to_numpy()
rust_data
state | decision | |
---|---|---|
0 | 0.0 | 0.0 |
1 | 0.0 | 0.0 |
2 | 1.0 | 0.0 |
3 | 2.0 | 0.0 |
4 | 3.0 | 0.0 |
... | ... | ... |
8255 | 68.0 | 0.0 |
8256 | 68.0 | 0.0 |
8257 | 69.0 | 0.0 |
8258 | 69.0 | 0.0 |
8259 | 69.0 | 0.0 |
8260 rows × 2 columns
#Parameters used in the original paper
T = 90
β = 0.9999
scale=1e-3
γ = np.euler_gamma
x = np.arange(T)
data = rust_data.to_numpy()
#θ is a vector of two parameters that we wish to estimate. θ_0 denotes the cost of reparations per mileage
#and θ_1 denotes the replacement cost.
def u(x, d, θ):
if d==1: #replace
uval = -θ[1] - (scale*x*θ[0])
elif d==0: #not replace
uval = -(scale*x*θ[0])
return uval
Computational trick to avoid overflow issues:
\begin{align} EV_{new} &= \log(e^{rep} + e^{wait}) \\ EV_{new} &= \log(e^{rep} + e^{wait}) - \log(e^{EV}) + EV \\ EV_{new} &= \log(e^{rep - EV} + e^{wait - EV}) + EV \end{align}
def return_EV(θ, tol=1e-3, maxiter=300000, verbose=False):
#fixed-point algorithm
n_iter=0
EV = np.zeros((T, 1))
error = 1e5
while error > tol or n_iter > maxiter:
EV_new = np.empty((T, 1))
u_not_replace = u(x, 0, θ).reshape(T, 1) + β*EV
u_replace = u(x[0], 1, θ) + β*EV[0]
EV_new = P@(np.log(np.exp(u_replace - EV) + np.exp(u_not_replace - EV)) + EV)
error = np.max(np.abs(EV - EV_new))
n_iter +=1
if verbose==True:
print(f'Iteration {n_iter}, error: {error}')
EV = EV_new.copy()
return EV
def return_ccp(θ):
EV = return_EV(θ)
state_ccp_map = np.empty((T, 2))
for state in x:
u_not_replace = u(state, 0, θ) + β*(P[state, :] @ EV)
u_replace = u(0, 1, θ) + β*EV[0]
proba_not_replace = (1/(1 + np.exp(u_replace - u_not_replace)))[0]
proba_replace = (1/(1 + np.exp(u_not_replace - u_replace)))[0]
state_ccp_map[state, :] = np.array([proba_not_replace, proba_replace])
return state_ccp_map
def ll(θ):
CCP = return_ccp(θ)
logL = 0
for s, d in data:
if int(d)==0:
logL += np.log(CCP[int(s)][0])
elif int(d)==1:
logL += np.log(CCP[int(s)][1])
return -logL
#we start close to the optimal solution to speed things up, but feel free to try alternate starting values.
#It will already be very long, so I'd advise not running it.
theta_star = minimize(ll, x0=np.array([2, 10])).x
theta_star
array([ 2.61806004, 10.03903027])