© 2018-2026 by Alfred Galichon. Past and present support from NSF grant DMS-1716489, ERC grant CoG-866274 are acknowledged, as well as inputs from contributors listed here.
If you reuse material from this masterclass, please cite as:
Alfred Galichon, 'math+econ+code' masterclass. https://github.com/math-econ-code/mec_notebooks
Rust, John (1987), Optimal replacement of GMC bus engines: an empirical model of Harold Zurcher. Econometrica.
Ferrall, C. (2023). Was Harold Zurcher myopic after all? Replicating Rust's engine replacement estimates. Journal of Applied Econometrics, 38(7), 1093–1100. https://doi.org/10.1002/jae.3001
Galichon, Alfred. Discrete Choice Models: Mathematical Methods, Data Science, and Econometrics. Princeton Univeristy Press.
Galichon, Alfred. "math-econ-code." GitHub. https://github.com/math-econ-code.
Lin, Ranie. "Rust-1987-Replication." GitHub. https://github.com/ranielin/Rust-1987-Replication.
Before detailing the model, it is worth emphasizing this paper's foundational role in the dynamic discrete choice literature. Despite its relative simplicity, the structural framework introduced by Rust (1987) remains widely utilized today.The economic agent in this model is Harold Zurcher, superintendent of maintenance at the municipal bus depot in Madison, Wisconsin.
Each period, he faces a discrete binary choice: replace a given bus engine or keep it running. This constitutes an optimal stopping problem, where the objective is to find a decision rule that perfectly balances the tradeoff between minimizing ongoing maintenance costs and preventing unexpected engine failures. The optimal solution takes the form of a threshold rule, characterized by a critical mileage cutoff, $x^*$.
Rust formalizes a regenerative optimal stopping problem.
The State Variable ($x_t$): The accumulated mileage on the bus engine at time $t$.
The Decision ($i_t$): Each month, Zurcher faces a discrete binary choice:
Every month $t$, Zurcher observes the state of the bus engine, primarily its accumulated mileage $x_t$, and a set of unobserved shocks $\epsilon_t(0)$ and $\epsilon_t(1)$ that affect his costs and then make a choice.
$$u(x_t, i_t, \theta_1) + \epsilon_t(i_t) = \begin{cases} -c(x_t, \theta_1) + \epsilon_t(0) & \text{if } i_t = 0 \text{ (Keep)} \\ -RC - c(0, \theta_1) + \epsilon_t(1) & \text{if } i_t = 1 \text{ (Replace)} \end{cases}$$Where:
Zurcher is forward-looking and wants to maximize his expected discounted utility over an infinite horizon.
$$V_\theta(x_t, \epsilon_t) = \max_{i \in \{0,1\}} \left[ u(x_t, i, \theta_1) + \epsilon_t(i) + \beta EV_\theta(x_t, \epsilon_t, i) \right]$$The parameter vector we are estimating is $\theta = (\beta, \theta_1, RC, \theta_3)$ via Maximum Likelihood Estimation:
$\theta_3$ (Transition Probabilities): How mileage evolves from month to month.
$\theta_1$ (Maintenance Cost Parameters): How quickly regular costs scale with mileage.
$RC$ (Replacement Cost Parameter): The cost of an engine swap.
$\beta$ (Discount Factor): How much Zurcher cares about the future. (Rust tests both a myopic model where $\beta = 0$ and a dynamic model where $\beta = 0.9999$ ).
We briefly describe the dataset. It comprises maintenance records for 162 buses spanning from December 1974 to May 1985. The data provides monthly odometer readings and a detailed maintenance diary for each bus. Of the various maintenance tasks recorded, Rust isolates major engine overhauls and replacements. Consequently, the core objective is to estimate a structural model capable of predicting the time and mileage at which an engine replacement becomes optimal."
We started with raw text files (g870.asc, rt50.asc, etc.) containing unlabelled, column-major matrices.
In Rust's data, Group 1 buses were tracked for exactly 36 months. Group 2 was tracked for 60 months, Group 3 for 81 months, and Group 4 for 128 months, we will use it to track bus by bus the data.
To make the infinite-horizon dynamic programming computationally tractable, continuous mileage cannot be used as the state variable. Instead, the continuous mileage since the last replacement is discretized into 5,000-mile bins (e.g., state 0 is 0-5,000 miles, state 1 is 5,000-10,000 miles, etc.). The following cleaning functions identify the exact months where engine replacements ($i_t=1$) occurred, reset the continuous mileage tracker, and bin the data into these discrete states.
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import scipy.optimize as opt
import scipy.stats as stats
!pip install mec
from mec.data import load_Rust_data
Requirement already satisfied: mec in c:\users\alfre\anaconda3\lib\site-packages (0.210) Requirement already satisfied: gurobipy in c:\users\alfre\anaconda3\lib\site-packages (from mec) (11.0.2)
[notice] A new release of pip is available: 25.1.1 -> 26.1.1 [notice] To update, run: C:\Users\alfre\anaconda3\python.exe -m pip install --upgrade pip
# We used NumPy to read the flat files and reshape them into 2D arrays
# using Fortran order (order='F').
base_url = 'https://raw.githubusercontent.com/math-econ-code/mec_datasets/main/dynamicchoice_Rust/datafiles/'
def load_and_reshape(filepath, num_buses):
url = base_url + filepath
flat_data = np.loadtxt(url)
nrow = len(flat_data) // num_buses
return flat_data.reshape((nrow, num_buses), order='F')
def clean_dat(df_raw, group_num):
od_rep_1 = df_raw[5, :]
od_rep_2 = df_raw[8, :]
df = df_raw[11:, :]
rep_1 = (df >= od_rep_1) & (od_rep_1 > 0) #The engine is replaced once.
rep_2 = (df >= od_rep_2) & (od_rep_2 > 0) #The engine is replaced twice.
mileage_continuous = df.copy()
mileage_continuous -= (rep_1 & ~rep_2) * od_rep_1
mileage_continuous -= rep_2 * od_rep_2
x = np.floor(mileage_continuous / 5000) # We discretize into 5,000-mile bins.
diff_rep_1 = np.diff(rep_1.astype(int), axis=0)
diff_rep_2 = np.diff(rep_2.astype(int), axis=0)
i_inner = ((diff_rep_1 > 0) | (diff_rep_2 > 0)).astype(int)
zeros_row = np.zeros((1, df.shape[1]), dtype=int)
i = np.vstack([i_inner, zeros_row])
num_months, num_buses = df.shape
bus_ids = np.tile(np.arange(1, num_buses + 1), (num_months, 1))
time_idx = np.tile(np.arange(1, num_months + 1).reshape(-1, 1), (1, num_buses))
df_out = pd.DataFrame({
'group': group_num,
'bus_id': bus_ids.flatten(order='F'),
'time_idx': time_idx.flatten(order='F'),
'raw_odom': df.flatten(order='F'),
'mileage': mileage_continuous.flatten(order='F'),
'state_x': x.flatten(order='F'),
'replace_i': i.flatten(order='F')
})
df_out = df_out[df_out['raw_odom'] > 0].copy()
return df_out
df1_raw = load_and_reshape("g870.asc", 15)
df2_raw = load_and_reshape("rt50.asc", 4)
df3_raw = load_and_reshape("t8h203.asc", 48)
df4_raw = load_and_reshape("a530875.asc", 37)
df1 = clean_dat(df1_raw, group_num=1)
df2 = clean_dat(df2_raw, group_num=2)
df3 = clean_dat(df3_raw, group_num=3)
df4 = clean_dat(df4_raw, group_num=4)
bus_dat = pd.concat([df1, df2, df3, df4], ignore_index=True)
Table IIa summarizes the replacement data for the subsample of buses that underwent at least one engine replacement.
# It increments every time an engine is replaced for a specific bus
bus_dat['replac_id'] = bus_dat.groupby(['group', 'bus_id'])['replace_i'].cumsum()
# Shift_id by 1 so the month of replacement is grouped with the engine being replaced
bus_dat['replac_id'] = bus_dat.groupby(['group', 'bus_id'])['replac_id'].shift().fillna(0)
completed_replac = bus_dat.groupby(['group', 'bus_id', 'replac_id']).filter(lambda x: x['replace_i'].iloc[-1] == 1)
spell_summary = completed_replac.groupby(['group', 'bus_id', 'replac_id']).agg(
mileage_at_replacement=('mileage', 'max'),
elapsed_months=('mileage', 'size')
).reset_index()
# We generate Table IIa
table_IIa = spell_summary.groupby('group').agg(
Mileage_Max=('mileage_at_replacement', 'max'),
Mileage_Min=('mileage_at_replacement', 'min'),
Mileage_Mean=('mileage_at_replacement', 'mean'),
Mileage_Std=('mileage_at_replacement', 'std'),
Months_Max=('elapsed_months', 'max'),
Months_Min=('elapsed_months', 'min'),
Months_Mean=('elapsed_months', 'mean'),
Months_Std=('elapsed_months', 'std'),
Observations=('mileage_at_replacement', 'count')
).reset_index()
# We re-index to ensure Groups 1 and 2 show up with 0s (since they get filtered out by having no replacements)
all_groups = pd.DataFrame({'group': [1, 2, 3, 4]})
table_IIa = all_groups.merge(table_IIa, on='group', how='left').fillna(0)
# @title
# We round the statistics for clean display
table_IIa = table_IIa.round({'Mileage_Mean': 0, 'Mileage_Std': 0, 'Months_Mean': 1, 'Months_Std': 1})
print("TABLE IIa: SUMMARY OF REPLACEMENT DATA (Groups 1-4)")
print(table_IIa.to_string(index=False))
TABLE IIa: SUMMARY OF REPLACEMENT DATA (Groups 1-4)
group Mileage_Max Mileage_Min Mileage_Mean Mileage_Std Months_Max Months_Min Months_Mean Months_Std Observations
1 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
2 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
3 273369.0 121326.0 198853.0 37864.0 69.0 30.0 52.9 11.0 27.0
4 387282.0 120709.0 256417.0 65317.0 112.0 27.0 71.2 22.9 33.0
Table IIb presents the right-censored data, representing the subsample of buses whose engines were not replaced during the observation period.
bus_dat['elapsed_months'] = bus_dat.groupby(['group', 'bus_id', 'replac_id']).cumcount() + 1
replaced_obs = bus_dat[bus_dat['replace_i'] == 1]
final_rows = bus_dat.groupby(['group', 'bus_id', 'replac_id']).tail(1)
kept_obs = final_rows[final_rows['replace_i'] == 0]
replacements_per_bus = bus_dat.groupby(['group', 'bus_id'])['replace_i'].sum()
never_replaced_buses = replacements_per_bus[replacements_per_bus == 0].reset_index()
# We filter our main dataset to only include these censored buses
censored_bus_dat = bus_dat.merge(never_replaced_buses[['group', 'bus_id']], on=['group', 'bus_id'])
# Since they were never replaced, their highest mileage and total elapsed time will be located in the very last row of their data sequence.
final_observations = censored_bus_dat.groupby(['group', 'bus_id']).tail(1)
# We generate Table IIb
table_IIb = final_observations.groupby('group').agg(
Mileage_Max=('mileage', 'max'),
Mileage_Min=('mileage', 'min'),
Mileage_Mean=('mileage', 'mean'),
Mileage_Std=('mileage', 'std'),
Months_Max=('elapsed_months', 'max'),
Months_Min=('elapsed_months', 'min'),
Months_Mean=('elapsed_months', 'mean'),
Months_Std=('elapsed_months', 'std'),
Observations=('bus_id', 'count')
).reset_index()
table_IIb = table_IIb.round({
'Mileage_Max': 0, 'Mileage_Min': 0, 'Mileage_Mean': 0, 'Mileage_Std': 0,
'Months_Max': 0, 'Months_Min': 0, 'Months_Mean': 1, 'Months_Std': 2
})
print("TABLE IIb: CENSORED DATA (Groups 1-4)")
print(table_IIb.to_string(index=False))
TABLE IIb: CENSORED DATA (Groups 1-4)
group Mileage_Max Mileage_Min Mileage_Mean Mileage_Std Months_Max Months_Min Months_Mean Months_Std Observations
1 120151.0 65643.0 100117.0 12929.0 25 25 25.0 0.0 15
2 161748.0 142009.0 151182.0 8530.0 49 49 49.0 0.0 4
3 280802.0 199626.0 250766.0 21325.0 70 70 70.0 0.0 21
4 352450.0 310910.0 337222.0 17802.0 117 117 117.0 0.0 5
Following Rust, we restrict our focus to bus groups 1 through 4. The buses in these groups were the newest in the fleet and were used on the most active routes. Rust isolates these specific groups for two main reasons: the availability of reliable replacement cost data, and the fact that the monthly mileage distributions for each bus are relatively homogeneous within each group. Because the estimation procedure allows for between-group heterogeneity but assumes within-group homogeneity, selecting the most uniform groups is theoretically justified.
# @title
plt.figure(figsize=(10, 8))
plt.scatter(
kept_obs['mileage'] / 1000,
kept_obs['elapsed_months'],
marker='.',
color='black',
label=f'keep ({len(kept_obs)} obs)'
)
plt.scatter(
replaced_obs['mileage'] / 1000,
replaced_obs['elapsed_months'],
marker='+',
color='black',
s=60,
label=f'replace ({len(replaced_obs)} obs)'
)
plt.title('Bus Replacement Data: Full Sample\n+ = replace, . = keep', fontsize=14)
plt.xlabel('Mileage since last replacement (Thousands)', fontsize=12)
plt.ylabel('Elapsed time (months)', fontsize=12)
plt.xlim(0, 400)
plt.ylim(0, 130)
plt.xticks(range(0, 401, 100))
plt.yticks(range(0, 131, 10))
plt.tight_layout()
plt.show()
Figure 1 illustrates the considerable variation in both the time and mileage at which bus engine replacements occur, ranging roughly from 125,000 to over 300,000 miles. This wide dispersion contradicts the naive hypothesis of a single, strictly deterministic stopping barrier. Consequently, this variance motivates the core premise of Rust's structural model: Zurcher's replacement decisions are not based solely on observed odometer readings, but are also driven by other state variables that are known to him but remain unobserved by the econometrician.
Since we discretized mileage into 5,000-mile bins (90 intervals of 5,000 miles, same as Rust)., a bus will typically "jump" by 0, 1, or 2 bins in a single month. We define $\theta_3$ as the set of transition probabilities:
def get_mileage_transitions(df):
df = df.sort_values(by=['bus_id', 'time_idx']).copy()
df['prev_state_x'] = df.groupby('bus_id')['state_x'].shift(1)
df['prev_replace_i'] = df.groupby('bus_id')['replace_i'].shift(1)
valid_transitions = df.dropna(subset=['prev_state_x']).copy()
jumps = np.where(
valid_transitions['prev_replace_i'] == 1,
valid_transitions['state_x'],
valid_transitions['state_x'] - valid_transitions['prev_state_x']
)
valid_jumps = jumps[jumps >= 0]
valid_jumps = np.clip(valid_jumps, 0, 2)
return valid_jumps.astype(int)
def estimate_transition_parameters(jumps):
n_obs = len(jumps)
counts = np.bincount(jumps, minlength=3)
theta = counts / n_obs
ll = np.sum(counts * np.log(np.where(theta > 0, theta, 1e-10)))
return theta, ll, n_obs
# @title
# We generate Table V: Within Group Estimates
print("TABLE V: WITHIN GROUP ESTIMATES OF MILEAGE PROCESS")
print("-" * 65)
print(f"{'Group':<10} | {'Theta 31':<12} | {'Theta 32':<12} | {'Theta 33 | Restricted Log-Like'}")
print("-" * 65)
group_results = {}
for g in [1, 2, 3, 4]:
jumps = get_mileage_transitions(bus_dat[bus_dat['group'] == g])
theta, ll, n_obs = estimate_transition_parameters(jumps)
group_results[g] = {'theta': theta, 'll': ll, 'n_obs': n_obs, 'jumps': jumps}
print(f"Group {g:<4} | {theta[0]:.4f} | {theta[1]:.4f} | {theta[2]:.4f} | {ll:.2f}")
TABLE V: WITHIN GROUP ESTIMATES OF MILEAGE PROCESS ----------------------------------------------------------------- Group | Theta 31 | Theta 32 | Theta 33 | Restricted Log-Like ----------------------------------------------------------------- Group 1 | 0.1972 | 0.7889 | 0.0139 | -203.99 Group 2 | 0.3906 | 0.5990 | 0.0104 | -138.57 Group 3 | 0.3149 | 0.6751 | 0.0100 | -2235.67 Group 4 | 0.3996 | 0.5876 | 0.0128 | -3153.83
Do individual buses within the same group accumulate mileage at the same rate?
Our tests fail to reject the null hypothesis, indicating that a restricted model (where transition probabilities are identical across all buses) fits the data just as well as an unrestricted model. This validates Rust's assumption of within-group homogeneity.
# We need a 3-state transition function (as we discretizes with 5,000 miles bins (no city bus will ever drive more than 15,000 miles in a given month)).
def get_mileage_transitions_3state(df):
df = df.sort_values(by=['group', 'bus_id', 'time_idx']).copy()
df['prev_state_x'] = df.groupby(['group', 'bus_id'])['state_x'].shift(1)
df['prev_replace_i'] = df.groupby(['group', 'bus_id'])['replace_i'].shift(1)
valid_transitions = df.dropna(subset=['prev_state_x']).copy()
jumps = np.where(
valid_transitions['prev_replace_i'] == 1,
valid_transitions['state_x'],
valid_transitions['state_x'] - valid_transitions['prev_state_x']
)
valid_jumps = jumps[jumps >= 0]
valid_jumps = np.clip(valid_jumps, 0, 2)
return valid_jumps.astype(int)
def estimate_transition_parameters_3state(jumps):
n_obs = len(jumps)
counts = np.bincount(jumps, minlength=3)
theta = counts / n_obs
ll = np.sum(counts * np.log(np.where(theta > 0, theta, 1e-10)))
return theta, ll
# We generate Table VI: Between Group Estimates
print("\nTABLE VI: BETWEEN GROUP ESTIMATES OF MILEAGE PROCESS")
print("-" * 85)
print(f"{'Groups Pooled':<15} | {'Theta 30':<8} | {'Theta 31':<8} | {'Theta 32':<8} | {'Restricted LL':<15} | {'LR Stat':<10}")
print("-" * 85)
pooled_definitions = {
'1, 2, 3, 4': [1, 2, 3, 4],
'1, 2, 3': [1, 2, 3]
}
for label, groups in pooled_definitions.items():
# Restricted Model
pooled_data = bus_dat[bus_dat['group'].isin(groups)]
combined_jumps = get_mileage_transitions_3state(pooled_data)
theta_restricted, ll_restricted = estimate_transition_parameters_3state(combined_jumps)
# Unrestricted Model (We estimate separately for every individual bus)
ll_unrestricted = 0
# We identify unique buses by their group and bus_id.
unique_buses = pooled_data[['group', 'bus_id']].drop_duplicates()
num_buses = len(unique_buses)
for _, row in unique_buses.iterrows():
# We isolate the exact bus using both identifiers
bus_jumps = get_mileage_transitions_3state(
pooled_data[(pooled_data['group'] == row['group']) & (pooled_data['bus_id'] == row['bus_id'])]
)
if len(bus_jumps) > 0:
_, ll_bus = estimate_transition_parameters_3state(bus_jumps)
ll_unrestricted += ll_bus
# LR Test Statistic
lr_stat = -2 * (ll_restricted - ll_unrestricted)
# Degrees of Freedom
df = (num_buses - 1) * 3
p_val = stats.chi2.sf(lr_stat, df)
print(f"{label:<15} | {theta_restricted[0]:.4f} | {theta_restricted[1]:.4f} | {theta_restricted[2]:.4f} | {ll_restricted:<15.2f} | {lr_stat:.2f}")
print(f" LR Test (df={df}): {p_val:.4f}\n")
TABLE VI: BETWEEN GROUP ESTIMATES OF MILEAGE PROCESS ------------------------------------------------------------------------------------- Groups Pooled | Theta 30 | Theta 31 | Theta 32 | Restricted LL | LR Stat ------------------------------------------------------------------------------------- 1, 2, 3, 4 | 0.3561 | 0.6323 | 0.0116 | -5785.82 | 335.76 LR Test (df=309): 0.1416 1, 2, 3 | 0.3077 | 0.6819 | 0.0104 | -2592.90 | 166.64 LR Test (df=198): 0.9489
Group 4 clearly accumulates mileage at a different rate. Therefore, to maximize our sample size without introducing bias, we will pool Groups 1, 2, and 3 for our main estimation and exclude Group 4.
Our maximum likelihood estimation of the dynamic discrete choice model requires a nested fixed-point algorithm.
The Inner Loop: Solving the Fixed Point
For Zurcher to make a decision, he needs to know the expected long-term costs of keeping versus replacing the engine. For a given set of cost parameters, we iteratively solve the Bellman equation to compute the expected value function (the fixed point). This quantifies the forward-looking expected discounted utility of the decision-maker. Once the fixed point is found, the inner loop can translate those expected values into Choice Probabilities.
The Outer Loop: Maximum Likelihood
We use the derived expected value functions and unobserved shock distributions to calculate the conditional probability of Zurcher choosing to replace or keep the engine. We construct the log-likelihood function and use an optimization routine to maximize it with respect to the cost parameters ($\theta_1$ and $RC$).
We run this estimation routine for different assumptions about how much Zurcher cares about the future. Specifically, test a myopic model (setting the discount factor $\beta = 0$) and a dynamic model (setting $\beta = 0.9999$).
# --
# INNER LOOP: Value Function Iteration (Contraction Mapping)
# --
def solve_dynamic_program(RC, theta_1, theta_3, beta, n_states=90, tol=1e-6):
"""
Computes the Expected Value function for a given set of parameters.
"""
# We use a Linear Specification: c(x) = theta_1 * x)
states = np.arange(n_states)
maint_cost = theta_1 * states * 0.001
# We build Transition Matrices based on Stage 1 estimates (theta_3)
P_keep = np.zeros((n_states, n_states))
P_replace = np.zeros((n_states, n_states))
for i in range(n_states):
P_keep[i, i] = theta_3[0]
if i + 1 < n_states:
P_keep[i, i+1] = theta_3[1]
else:
P_keep[i, i] += theta_3[1]
if i + 2 < n_states:
P_keep[i, i+2] = theta_3[2]
else:
P_keep[i, -1] += theta_3[2]
P_replace[i, 0] = theta_3[0]
P_replace[i, 1] = theta_3[1]
P_replace[i, 2] = theta_3[2]
# We run the Contraction Mapping
EV = np.zeros(n_states) # Initial guess for Expected Value
max_iter = 2000
for _ in range(max_iter):
# Choice-specific value functions: V(x, d) = u(x, d) + beta * E[EV(x')]
# If Keep: Pay maintenance cost, state transitions normally
v_keep = -maint_cost + beta * (P_keep @ EV)
# If Replace: Pay RC and maintenance for state 0, state resets
v_replace = -RC - maint_cost[0] + beta * (P_replace @ EV)
# Log-sum formula
# EV_new = ln(exp(v_keep) + exp(v_replace))
max_v = np.maximum(v_keep, v_replace)
EV_new = max_v + np.log(np.exp(v_keep - max_v) + np.exp(v_replace - max_v))
if np.max(np.abs(EV_new - EV)) < tol:
break
EV = EV_new
return EV, v_keep, v_replace
# --
# OUTER LOOP: Maximum Likelihood Estimation
# --
def negative_log_likelihood(params, df, theta_3, beta):
"""
Calculates the negative log-likelihood of the observed choices given parameters.
"""
RC, theta_1 = params
# We run the inner loop to find Harold Zurcher's optimal policy for these params
EV, v_keep, v_replace = solve_dynamic_program(RC, theta_1, theta_3, beta)
state = df['state_x'].values.astype(int)
choice = df['replace_i'].values.astype(int)
# We get the choice-specific values for the observed states in the data
v_k = v_keep[state]
v_r = v_replace[state]
# We calculate log choice probabilities
max_v = np.maximum(v_k, v_r)
log_sum = max_v + np.log(np.exp(v_k - max_v) + np.exp(v_r - max_v))
# P(d=0 | x) and P(d=1 | x)
log_prob_keep = v_k - log_sum
log_prob_replace = v_r - log_sum
# We sum the log-likelihood of the actual choices made in the data
ll_array = np.where(choice == 1, log_prob_replace, log_prob_keep)
ll = np.sum(ll_array)
return -ll
# -- We replicate Table IX for Beta = 0.9999 and Beta = 0.0 --
print("-" * 85)
print(f"{'TABLE IX: STRUCTURAL ESTIMATES (LINEAR COST SPECIFICATION)':^85}")
print("-" * 85)
groups_to_run = {
'Groups 1, 2, 3': [1, 2, 3],
'Group 4': [4],
'Groups 1, 2, 3, 4': [1, 2, 3, 4]
}
betas_to_test = [0.9999, 0.0]
results = {beta: {} for beta in betas_to_test}
for beta in betas_to_test:
for name, grps in groups_to_run.items():
df_pool = bus_dat[bus_dat['group'].isin(grps)]
jumps = get_mileage_transitions_3state(df_pool)
theta_3, ll_stage1 = estimate_transition_parameters_3state(jumps)
res = opt.minimize(
negative_log_likelihood,
[10.0, 2.0],
args=(df_pool, theta_3, beta),
method='L-BFGS-B',
bounds=[(0, None), (0, None)],
options={'ftol': 1e-5}
)
results[beta][name] = {
'RC': res.x[0],
'theta_1': res.x[1],
'theta_3': theta_3,
'll_stage2': -res.fun,
'll_stage1': ll_stage1,
'full_ll': ll_stage1 - res.fun
}
# We display Table IX
for beta in betas_to_test:
print(f"\n>> DISCOUNT FACTOR: Beta = {beta}")
print(f"{'Parameter':<25} | {'Groups 1, 2, 3':<17} | {'Group 4':<17} | {'Groups 1, 2, 3, 4'}")
print("-" * 85)
print(f"{'Replacement Cost (RC)':<25} | {results[beta]['Groups 1, 2, 3']['RC']:<17.4f} | {results[beta]['Group 4']['RC']:<17.4f} | {results[beta]['Groups 1, 2, 3, 4']['RC']:.4f}")
print(f"{'Maintenance Cost (Theta 1)':<25} | {results[beta]['Groups 1, 2, 3']['theta_1']:<17.4f} | {results[beta]['Group 4']['theta_1']:<17.4f} | {results[beta]['Groups 1, 2, 3, 4']['theta_1']:.4f}")
print("-" * 85)
if beta == 0.9999:
print(f"{'Theta 30 (Jump = 0)':<25} | {results[beta]['Groups 1, 2, 3']['theta_3'][0]:<17.4f} | {results[beta]['Group 4']['theta_3'][0]:<17.4f} | {results[beta]['Groups 1, 2, 3, 4']['theta_3'][0]:.4f}")
print(f"{'Theta 31 (Jump = 1)':<25} | {results[beta]['Groups 1, 2, 3']['theta_3'][1]:<17.4f} | {results[beta]['Group 4']['theta_3'][1]:<17.4f} | {results[beta]['Groups 1, 2, 3, 4']['theta_3'][1]:.4f}")
print("-" * 85)
print(f"{'Full Log-Likelihood':<25} | {results[beta]['Groups 1, 2, 3']['full_ll']:<17.2f} | {results[beta]['Group 4']['full_ll']:<17.2f} | {results[beta]['Groups 1, 2, 3, 4']['full_ll']:.2f}")
-------------------------------------------------------------------------------------
TABLE IX: STRUCTURAL ESTIMATES (LINEAR COST SPECIFICATION)
-------------------------------------------------------------------------------------
>> DISCOUNT FACTOR: Beta = 0.9999
Parameter | Groups 1, 2, 3 | Group 4 | Groups 1, 2, 3, 4
-------------------------------------------------------------------------------------
Replacement Cost (RC) | 11.7247 | 9.9281 | 9.7696
Maintenance Cost (Theta 1) | 4.7905 | 2.2201 | 2.6164
-------------------------------------------------------------------------------------
Theta 30 (Jump = 0) | 0.3077 | 0.3996 | 0.3561
Theta 31 (Jump = 1) | 0.6819 | 0.5876 | 0.6323
-------------------------------------------------------------------------------------
Full Log-Likelihood | -2725.28 | -3317.42 | -6086.07
>> DISCOUNT FACTOR: Beta = 0.0
Parameter | Groups 1, 2, 3 | Group 4 | Groups 1, 2, 3, 4
-------------------------------------------------------------------------------------
Replacement Cost (RC) | 8.3087 | 7.6417 | 7.3154
Maintenance Cost (Theta 1) | 110.1428 | 71.6201 | 70.4670
-------------------------------------------------------------------------------------
Full Log-Likelihood | -2727.66 | -3319.31 | -6092.54
It is worth noting that while our estimates show some discrepancies with Rust's original 1987 results, they are entirely consistent with modern replication efforts (e.g., Ferrall, 2023). To formalize this, Rust conducts a 'myopia test' by restricting the discount factor to $\beta = 0$. This essentially assumes Zurcher acts myopically by replacing an engine only when its immediate maintenance cost exceeds the replacement cost, completely ignoring future maintenance costs expectations. However, the forward-looking specification (with a high discount factor) provides a significantly better fit to the data, allowing us to reject the null hypothesis of a myopic decision-maker at the 5% significance level.