[ Search |  Up Level |  Project home |  Index |  Class hierarchy ]

 RustEmet1987.ox

RustEmet1987.ox replicates Table IX, Column 2, Rows 1 and 2 in Rust (Econometrica 1987), the celebrated bus engine replacement model.

Skip down to documentation of items defined in DP.ox

1. Overview
2. Model
3. Implementation
4. Remarks
5. Replication Output
6. Simulation

1. Overview

The program produces engine replacement probabilities that match those shown in Figure 3 on page 1025 of the paper. The code segments are not exactly the same as in the Ox file, which includes additional code beyond the basics shown here.

The paper models the decision to replace a bus engine as a function of mileage and unobserved IID shocks. Replacement is a one-time cost that resets the odometer to x=0. The cost of monthly maintenance increases linearly with x.

See RustEmet1987mle to see estimation using data as in the original paper.

2. The Model

An example of Rust class DDP:

$\zeta$ specification: ex ante extreme value shocks with $\rho = 1.0.$

Clock: Ergodic (infinite horizon with no terminal or absorbing states)

Action: $\alpha = (d)$, d.N = 2, the binary decision to replace the engine or not.

Endogenous State

$\theta = (x),$ x.N = 90.

x is the mileage reading on the bus at its monthly maintenance, grouped into 5000 mile categories.

Transition: mileage next month is 0, 1 or 2 mileage bins more.

    x' = (1-d)x    w/ prob. θ31
         (1-d)x +1 w/ prob. θ32
         (1-d)x +2 w/ prob. 1-θ31-θ32

If x = 88 or 89 then the increment probabilities changed accordingly.

Semi-Endogenous: none

Endogenous: none

Utility:

The cost of regular maintenance or engine replacement.

U = -[ dRC + (1-d)θ1mx ] + n + zd    (5.1) p. 1015

z_d is the standard IID extreme value shock. (In the paper it is denoted ε)

Parameters:

Parameter	Row 1	Row 2
δ	0.9999;	0
RC	10.07	7.6538
θ1	2.293	71.5133
θ3	< 0.3919,0.5953>	< 0.3919,.5953>

The article suggests x enters the code in terms of thousands of miles (0 to 445,000), but I was not able to replicate results based on that. I was able to replicate replacement probabilities only when x takes on values between 0 and 89.

RC is the bus engine replacement cost. θ₁ is the operating cost given mileage x. The parameter m=0.001 is a scaling factor, which appears to rescale odometer readings. However, results were replicated when x is the bin number (0 to 89) not the odometer category.

Normalization

The paper refers to a normalization when displaying V() (see Figure 2, p. 1025). The normalization is not explained, but my calculations seem to explain why it is necessary. The closed form extreme-value value function iteration involves taking exponents of negative values (operating costs). When δ=0.9999 and n=0 above the result is numerical overflow in exp() With n>0 set properly the arguments are kept closer to 0. In the replication the normalization n is set to the maintenance cost for x=89/2, the median mileage category.

3. Implementation in DDP
Declarations in the Header File
Header File

Executable Code in the Ox File
Ox File

4. Remarks
Zurcher is derived from the Rust class
code>Rust class embeds an infinite horizon clock and a single binary action (d). It defines constants and declares a static member x>. Zurcher then declares the required method Utility().
Utility must have that name and cannot be static.
The form of the expression returned as the value should look related to the mathematical expression U() above. The parameters are static members of the class (they are the shared by each state θ). The current value of x is accessed as CV(x). The value of x is set to be the correct value at θ before Utility() is called.
Utility must return a column vector of values corresponding to the feasible action vector A[Aind]. In this case the action variable is simply d, and engine replacement is always feasible.
In mathematical notation, we would write the choice set as $\alpha \in A(\theta)$. Rather than accessing that directly, the column of values for d is accessed using CV(d). This makes the code closer to the notation and more general since it will work even if other actions are added to the model.
The tasks required to set up and solve the model are in Zurcher::Run(). Like Run(), a user's program must do these tasks in order:
1. Call Initialize(new Zurcher()) for the parent DDP class. Since Zurcher is of type Rust, the program is calling Rust::Initialize(), which will set the clock and add the binary action variable to the model, d.
2. Create state variables (using the new operator) and add them to the model.
3. Call CreateSpaces();.
The state variable x is of type Renewal, which captures the process in Rust (1987) but is more general than in the empirical specification.
To properly capture the transition it must be coordinated with the choice variable d. To do this d is passed to new Renewal() as an argument, along with how many categories x takes on and the vector of increment probabilities. The length of the vector determines how different states are feasible next period. In the empirical work, 3 of 90 odometer readings are feasible next period.
Once CreateSpaces() has been called new actions and new state variables cannot be added to the model because the state spaces have be set up already.
To solve the model the following steps are carried out
1. Set the discount factor δ. (This can be done at any time after Initialize() before solving the model with that δ.)
2. Create a new instance of a solution Method, such as value function iteration.
3. Call Solve() for the solution method to compute $V(\theta)$ and $P^\star(\alpha|\theta)$.
4. Use the results, then if necessary change some settings or parameters and solve again.
The code inside the for() loop solves the model twice, for Row 1 and Row 2 parameters.
For each set of parameters the value function is solved through Bellman iteration by Emax->Solve(). It calls the Output() routine to capture results to be printed and compared with Figure 3 in the paper.
5. Output

       Output produced: niqlow/examples/output/RustEmet1987.txt

The values of EV do not match those shown in Figure 2 of the paper. This is due to Rust using a different normalization than the value used here. However, the choice probabilities appear to match very well those shown if Figure 3, as the image below suggests.

Comparison to Figure 3


The paper labels the x axis as thousands of miles since engine replacement. This means that x should take on values 0, 5, ..., 445. It also suggests that the factor on costs scales 300,000 to 300. However, using this scaling I was not able to replicate the choice probabilities. Only when x takes on values 0, ..., 89 and the scaling factor still equals 0.001 was I able to reproduce the choice probabilities.

6. Simulation
Another use of the model is to simulate data by creating a Panel object and simulating a sample of data from the process.
The call to Panel::Simulate creates a data set of 10 buses observed over 400 months each. The initial mileage for each bus will be drawn from the ergodic distribution.
The output file shows a simulated panel of buses, including the months and mileage at every replacement. The simulate routine returns a panel of observations.
But in the case of this model it is sufficient to show only months when an engine is replaced, the age of the bus when it occurs and the mileage reading. So the replication program deletes simulated months when no replacement occurred. It also appends the row number for the model.
For example, the first simulated bus data looks like this

Simulated data
Simulated data
  Fxed  path     t     n     T   Aj|   s11   s21     x     t    t'     r     f  |ai|     a
     0     0    92    50     0     0     0     0    50     0     0     0     0     1     1
     0     0   164    44     0     0     0     0    44     0     0     0     0     1     1
     0     0   248    53     0     0     0     0    53     0     0     0     0     1     1
     0     0   344    69     0     0     0     0    69     0     0     0     0     1     1

Bus 0 had four replacements, which occurred in months 92, 164, 248 and 344 of its lifetime. The mileage on the odometer at replacement was in the interval [5000x,5000(x+1)].
The column labeled Const is a placeholder for exogenous state variables.
If the model has no such variables then DDP inserts a constant there. This makes it possible to construct indices into the state space without requiring a check of zero dimensions.
Since the model is stationary there is no separate time variable. If this were a finite horizon or non-stationary environment then t would contain the age of the process.
DDP allows for random aging so t would not necessarily be equivalent to the months column.
The t' variable (t″ in the notes) is used during iteration on the value function. During simulation it is not used and would always have the value 0.

Author:

© 2011-2018 Christopher Ferrall

---

Documentation of Items Defined RustEmet1987.ox

 Zurcher

Public fields
bins	static const	Number of bins by table.
COL	static	column of Table
dfactor	static const	Discount factor by row.
DMETH	static	Discretization method: Month or Mileage
mfact	static const	scaling on cost
normalization	static	added to U() to avoid underflow
NX	static	# of bins, 90 or 175
parlist	static const	2x3x2x4 array of reported parameter estimates.
pars	static	parameter vector $\psi$.
rc	static	value of RC
ROW	static	row of table/column.
th1	static	value of $\theta_1$
x	static	mileage state, x
Public methods
Output	static
Run	static	Setup and solve the model for each discount factor and generate the figure.
SetSpec	static	Set the target of the replication.
Utility		Computes the one period linear cost as 2x1 vector $$U = dRC+(1-d)\theta_1 x + n.$$
Enumerations
RNpars		tags for table/column/row choices.
Zparams		tags for estimated parameters.

 Zurcher

 bins

Number of bins by table.

 COL

column of Table

 dfactor

Discount factor by row.

 DMETH

Discretization method: Month or Mileage

 mfact

scaling on cost

 normalization

added to U() to avoid underflow

 NX

# of bins, 90 or 175

 Output

 parlist

2x3x2x4 array of reported parameter estimates.

 pars

parameter vector $\psi$.

 rc

value of RC

 ROW

row of table/column.

 Run

Setup and solve the model for each discount factor and generate the figure.

 SetSpec

Set the target of the replication.

This allows the basic code to be re-used for different specifications.

Parameters:

targ	is array of integer values for the Table, COlumn and Row from the original paper.

 th1

value of $\theta_1$

 Utility

Computes the one period linear cost as 2x1 vector $$U = dRC+(1-d)\theta_1 x + n.$$

 x

mileage state, x