Introduction

Gauss code and data released in the 1990s are still available (Rust 2000), and several papers that have replicated elements of the original analysis. In particular Larsen et al. (2012) report replications of some original results using Matlab code while testing the validity of assuming the extreme value distribution for choice errors. Augirreibiria and Mira (2002) use the same data and provide code for their pseudo MLE procedure but do not report replications of Rust's results. Su and Judd (2012) conduct Monte Carlo experiments on their method using Rust's point estimates but do not replicate the estimates nor use the original data. The software package ruspy (Blesch 2019) implements both the original and Su and Judd solution techniques. Müller and Reich (2021) use the original data and estimate the original model but use a different technique and do not report agreement with the original maximum likelihood estimation.

This paper's modest goal is to replicate MLE estimates in Rust (1987) in order to demonstrate the software package niqlow described in Ferrall (2022). The replication is based not on bespoke code but rather high-level statements in niqlow. In particular, a main statistical inference reported in the paper is a likelihood ratio test of the null hypothesis of myopic behavior in replacing bus engines, which corresponds to estimates setting the discount factor at 0, versus a forward-looking model with the discount factor fixed at 0.9999. That specification requires the nested fixed point solution developed in the model. This replication reproduces both the MLE estimates of the parameter values and the results of these hypothesis tests.

Despite a 35 year gap in computing hardware and software, and a fundamentally different approach to coding the model, the solutions of the DP model are identical for a given set of parameter values. The MLE estimates do differ with the original values because of the difference in the transformation of the original raw data mentioned above. Using the data files in Larsen et al. (2012) an error in the discretization of the data was resolved. As with their analysis, re-estimation with the original data encoding results in identical estimates and test statistics. Using the data produced by the corrected encoding procedure results in different estimates. However, the shifts in likelihood are similar in the myopic and forward-looking versions of the models so test statistics are essentially the same.

Model and Its Solution

Panel B of Table 1 summarizes the model using the niqlow framework. A new class named Zurcher, is derived from the pre-defined Rust class. The state variable $x$ is an object of the Renewal class, following the term Rust used to describe the process. This objected-oriented approach makes it possible to build specialized and extended versions of existing models replacing only the code or features that differ.

Choice Probability Replication

Before turning to the ML estimation procedure, consider output at one of the estimated models in the original paper. Figure 3 in the paper compares $P^\star(1;x),$ the replacement probability, for discount factors of $\delta=0$ and $\delta=0.9999$. This does not involve the data. It is a check on the recoding of the Newton-Kantorovich algorithm, described in Rust (1988), to solve the fixed point in the value function along with the logit specification for the CCPs. 1 displays the original figure and the replication. The additional information in the original figure is the empirical hazard rate. The replication figure is shown alone on the right and is overlayed on the original scanned figure demonstrating they coincide exactly given the scan's resolution.

---

Left: original printed Figure 3 digitized with the replicated replacement probability curves overlayed (after matching the axes scales). Right: the replication shown alone.

---

Data and Estimates Replication

Observed Mileage Transitions

The transition or innovation $j = x_{s+1}-x_s(1-d_s)$ is constructed from the path of the discretized data $x_s$. The transitions are IID. Both $x$ and $d$ are observed so the transition probabilities $\theta_3$ form a multinomial problem not involving optimal replacement decisions (Rust's first stage estimation). 2 and 3 report the estimates of $\theta_3$ (fractions of transitions $j$) for bus groups 1-3 and group 4 from the original paper for 90 and 175 bins, respectively. The estimated transitions are compared with the values from the replication data.

It is in these transitions that a small discrepancy occurs. Both the original Gauss code and the Matlab code used in Larsen et al. (2012) compute the discrete $x_s$ using the ceiling function when it should be the floor function. For months without a bus replacement this has no effect because when $d_s=0$ the transition $j= x_{s+1}-x_{s}$ is the same regardless if both values are increased by one. However, in replacement months this results in $ j = x_{s+1}$ and this had minimum value 1 instead of 0. (It is not surprising that in months when an engine is replaced the added mileage is within the lowest mileage category which never occurred in the incorrect coding). The first stage transition probabilities are therefore incorrect which affects the second state estimates through $\theta_3$.

These differences are reported in final column labeled $Δ$ in 2 and 3. In the 4 different samples the multinomial likelihood on the transition probabilities differs by roughly 1% between the two codings of the discrete data. Since transition probabilities enter into the problem solved by the agent when determining conditional choices to replace an engine these differences will shift all MLE estimates.

Table 3 starts with the same continuous data but uses 175 discrete bins instead of 90. Since the width of the bins is smaller, the maximum jump increases from 2 to 5 bins. Comparing the replication to the original, the differences in likelihoods and cell percentages are somewhat larger than the corresponding values in Table 2. Multiplying the percentage differences by the number observations yields roughly 23 bus-month values are misclassified in each panel.

Likelihood

A sequence of outcomes for a single bus over $\hat{T}+1$ months is denoted $\{Y\}_{s=0}^{\hat T}$, where $Y_s = \{\al_s,\th_s\}$ is simply the decision and state vectors combined. In niqlow these are in general vectors, but in this model they reduce to two scalars $(d_s,x_s).$ The likelihood of a single sequence has the form: $$L\l(\ \hat\psi \, ; \,\{Y\}\ \r)\ =\ {\prod}_{s=0}^{\hat{T}}\quad \l\{ P^{\,\star}(\al_s; \th_s ) \l[ P\l( \th_{s+1}; \al_s, \th_s\r)\r]^{I\{s\lt \hat{T}\}}\r\}.\label{Lpath}$$

That is, the contribution in observed month $s$ is the CCP for the observed decision times the probability of the observed transition to the state next period. The transition is not observed in the last observed month for a bus. This function is generated automatically by niqlow from the model and the data. The general form has been specialized to this model for clarity. As in all "structural estimation" the observed data are inserted into the theoretical probabilities produced by the fixed point in the value function. In this model the transition probability is simply the multinomial parameter $\theta_{3j}.$

Equation (14.5) in Rust (1987) differs in two ways from \eqref{Lpath}. Translating terms and counting from 0 instead of 1, that equation would be written $$L^{1987}\l(\ \hat\psi \, ; \,\{Y\}\ \r)\ = \ {\prod}_{s=1}^{\hat{T}}\quad \l\{ P^{\,\star}(\al_s; \th_s )\l[ P\l(\th_{s}; \al_{s-1}, \th_{s-1}\r)\r]\r\}.\label{R14-5}$$ That is, \eqref{R14-5} removes the first choice probability at $s=0.$ The effect is coincidentally negligible. The estimated probability of not replacing an engine in most first observed months is very close to 1.0, because engines are only replaced many months into service. Thus near the maximum the difference between the log-likelihoods \eqref{Lpath} and \eqref{R14-5} is $\approx log(1) = 0.0$.

The second difference is also inconsequential for this case. \eqref{Lpath} multiplies the choice probability at $s$ by the transition to the current state $x_s$ not the transition from current state to the next state. This incorrect pairing is irrelevant since they are interchangeable in \eqref{Lpath}. However, for models with permanent or transitory unobserved states and actions they are not interchangeable. The likelihood should add across unobserved states as well as multiply across time so the conditional choice must be multiplied by the correct transition before summing as in \eqref{Lpath}, which niqlow does automatically for such models.

Estimation

The ML estimates were reported in original Table IX ($N=90$) and Table X ($N=175$). These models share the linear cost of bus maintenance already shown, and three different bus type samples (Groups 1-3, Group 4, and Groups 1-4 combined). The original estimates and the replications for the first two groups are reported in Tables 3 and 4. The values reported are "third stage" estimates in which all parameters are estimated including the first stage transition probabilities. Standard errors are also reported. The original procedure used BHHH and semi-closed forms for the derivatives of the value function with respect to structural parameters (derived in Rust 1988.) The replication uses BHHH and numerical central step derivatives. For $\delta=0.9999$ the convergence in the stage 3 estimates is weak: the norm of gradients remain large but progress ends. This combines large differences in the value function (which enter exponents) and precision in the convergence in the Newton-Kantorovich.

The likelihood values are within 1% of the original values. This magnitude is consistent with the exact match with choice probabilities for the same parameter values in Figure 1 and the coding discrepancy in the values of $j.$ Also note that the differences in log-likelihoods is approximately the same for $\delta=0$ as for $\delta=0.9999.$ In the former case convergence to the fix point is immediate, so any differences in tolerances or other aspects of the Newton-Kantorovich implementation have no effect on that case. Standard errors are all similar in magnitude between the original and replicated values.

Inference: Replicating the Myopia Test

One way to aggregate the overall effect of the data discrepancies and their effect on the estimates is to recompute the likelihood ratio test of the myopia restriction (Panel B versus Panel A in each table, one degree of freedom). The results are reported in Panel C of the tables. The p-values are essentially the same between the original and replicated results. The the myopic and forward-looking versions share the same transition data ($j$) and the effects on the distribution of $j$ is almost a constant shift in log-likelihood resulting in nearly identical LR statistics. This is likely because mileage $x$ and subsequent costs make modest changes month-to-month relative to the large reset when replacement occurs. Shifting the distribution of $j$ has little effect on the replacement dynamics.

Conclusion

Perhaps more important is the fact that the original model and estimation procedure have been independently replicated within a modern open-source platform. This reduces the cost of replication existing empirical DP results, and it makes modifying existing models for both teaching and research purposes substantially easier.

Acknowledgements

References

Supplemental Material

• Results were produced by niqlow version 5, available at
  https://github.com/ferrall/niqlow/releases/latest
• The code and data files specific to the results are contained in examples/replications. The relevant code files are:

RustEmet1987.h      RustEmet1987.ox         documentation
RustEmet1987mle.h   RustEmet1987mle.ox      documentation

• The data files are produced by RustEmet1987readdata2022.ox from the original ASCII files. Data sets used in the replication are:

  RustEmet1987_type1_col1_ALL.dta  RustEmet1987_type1_col2_ALL.dta RustEmet1987_type1_col3_ALL.dta
  

  Data identical to the original encoding using ceil() are contained in corresponding "type2" files.
• Raw output files for all three estimation stages are located in examples/output in files of the form R_{N}_{C}_{D}.txt where $N$ is either 0 (for 90) or 1 (for 175), C is column 0, 1 or 2, and D is either 0 ($\delta=0.9999$) or 1 ($\delta=0$).
• All estimation can be re-done using the interactive menu explained here:
https://ferrall.github.io/niqlow/examples/main.ox.html