From Rust (1988)

• An agent owns a fleet to buses

• Buses get old over time

• The older the bus is, the most costly it is to maintain

• The agent can decide to replace the bus engine with a new one, at a cost

• Dynamic trade-off

  • What is the best moment to replace the engine?

  • You don’t want to replace an engine too early

    • doesn’t change much

  • You don’t want to replace an engine too late

    • avoid unnecessary maintenance costs

• State: mileage of the bus

  \[s_t \in \lbrace 1, ..., 10 \rbrace \]

• State transitions: with probability \(\lambda\) the mileage of the bus increases

  \[ s_{t+1} = \begin{cases} \min \lbrace s_t + 1,10 \rbrace & \text { with probability } \lambda \newline s_t & \text { with probability } 1 - \lambda \end{cases} \]

  Note that \(\lambda\) does not depend on the value of the state

• Action: replacement decision \[ a_t \in \lbrace 0, 1 \rbrace \]

• Payoffs

  • Per-period maintenance cost

  • Cost of replacement \[ u\left(s_{t}, a_{t}, \epsilon_{1 t}, \epsilon_{2 t} ; \theta\right)= \begin{cases} -\theta_{1} s_{t}-\theta_{2} s_{t}^{2}+\epsilon_{0 t}, & \text { if } a_{t}=0 \newline -\theta_{3} + \epsilon_{1t}, & \text { if } a_{t}=1 \end{cases} \]

• Start with an initial expected value function \(V(s_t)=0\)

• Compute the alternative-specific value function \[ \bar V(s_t) = \begin{cases} -\theta_1 s_t - \theta_2 s_t^2 + \beta \Big[(1-\lambda) V(s_t) + \lambda V(\min \lbrace s_t+1,10 \rbrace ) \Big] , & \text { if } a_t=0 \newline -\theta_3 + \beta \Big[(1-\lambda) V(0) + \lambda V(1) \Big] , & \text { if } a_t=1 \end{cases} \]

• Compute the new expected value function \[ V'(a_t) = \log \Big( e^{\bar V(a_t|s_t=0)} + e^{\bar V(a_t|s_t=1)} \Big) \]

• Repeat until convergence

First we set the parameter values.

# Set parameters
θ = [0.13; -0.004; 3.1];
λ = 0.82;
β = 0.95;

Then we set the state space.

# State space
k = 10;
s = Vector(1:k);

First, we can compute static utility. \[ u\left(s_{t}, a_{t}, \epsilon_{1 t}, \epsilon_{2 t} ; \theta\right)= \begin{cases} -\theta_{1} s_{t}-\theta_{2} s_{t}^{2}+\epsilon_{0 t}, & \text { if } a_{t}=0 \newline -\theta_{3} + \epsilon_{1 t}, & \text { if } a_{t}=1 \end{cases} \]

function compute_U(θ::Vector, s::Vector)::Matrix
    """Compute static utility"""
    u1 = - θ[1]*s - θ[2]*s.^2       # Utility of not investing
    u2 = - θ[3]*ones(size(s))       # Utility of investing
    U = [u1 u2]                     # Combine in a matrix
    return U
end;

We can now set up the value function iteration

function compute_Vbar(θ::Vector, λ::Number, β::Number, s::Vector)::Matrix
    """Compute value function by Bellman iteration"""
    k = length(s)                                 # Dimension of the state space
    U = compute_U(θ, s)                           # Static utility
    index_λ = Int[1:k [2:k; k]];                  # Mileage index
    index_A = Int[1:k ones(k,1)];                 # Investment index
    γ = Base.MathConstants.eulergamma             # Euler's gamma

    # Iterate the Bellman equation until convergence
    Vbar = zeros(k, 2);
    Vbar1 = Vbar;
    dist = 1;
    iter = 0;
    while dist>1e-8
        V = γ .+ log.(sum(exp.(Vbar), dims=2))     # Compute value
        expV = V[index_λ] * [1-λ; λ]               # Compute expected value
        Vbar1 =  U + β * expV[index_A]             # Compute v-specific
        dist = max(abs.(Vbar1 - Vbar)...);         # Check distance
        iter += 1;
        Vbar = Vbar1                               # Update value function
    end
    return Vbar
end;

We can now solve for the value function.

# Compute value function
V_bar = compute_Vbar(θ, λ, β, s);

Now that we know how to compute the equilibrium, we can simulate the data.

function generate_data(θ::Vector, λ::Number, β::Number, s::Vector, N::Int)::Tuple
    """Generate data from primitives"""
    Vbar = compute_Vbar(θ, λ, β, s)             # Solve model
    ε = rand(Gumbel(0,1), N, 2)                 # Draw shocks
    St = rand(s, N)                             # Draw states
    A = (((Vbar[St,:] + ε) * [-1;1]) .> 0)      # Compute investment decisions
    δ = (rand(Uniform(0,1), N) .< λ)            # Compute mileage shock
    St1 = min.(St .* (A.==0) + δ, max(s...))    # Compute neSr state
    df = DataFrame(St=St, A=A, St1=St1)         # Dataframe
    CSV.write("../data/rust.csv", df)
    return St, A, St1
end;

We can now generate the data

# Generate data
N = Int(1e5);
St, A, St1 = generate_data(θ, λ, β, s, N);

How many investment decisions do we observe?

print("we observe ", sum(A), " investment decisions in ", N, " observations")

## we observe 19390 investment decisions in 100000 observations

What does the data look like?

# Read data
df = fread("../data/rust.csv")
kable(df[1:6,], digits=4)

• First we can estimate the value of lambda as the probability of mileage increase

• Conditional on not investing

• And not being in the last state (mileage cannot increase any more)

  \[ \hat \lambda = \mathbb E_n \Big[ (s_{t+1}-s_t) \mid a_{t}=0 \wedge s_{t}<10 \Big] \]

# Estimate lambda
Δ = St1 - St;
λ_ = mean(Δ[(A.==0) .& (St.<10)]);

print("Estimated lambda: $λ_ (true = $λ)")

## Estimated lambda: 0.8195976902431751 (true = 0.82)

• Take a parameter guess \(\theta_0\)

• Compute the alternative-specific value function \(\bar V(s_t ; \hat \lambda, \theta_0)\) by iteration

• Compute the implied choice probabilities

• Compute the likelihood \[ \mathcal{L}(\theta) = \prod_{t=1}^{T}\left(\hat{\operatorname{Pr}}\left(a=1 \mid s_{t}, \theta\right) \mathbb{1}\left(a_{t}=1\right)+\left(1-\hat{\operatorname{Pr}}\left(a=0 \mid s_{t}, \theta\right)\right) \mathbb{1}\left(a_{t}=0\right)\right) \]

• Repeat the above to find a minimum of the likelihood function

function logL_Rust(θ0::Vector, λ::Number, β::Number, s::Vector, St::Vector, A::BitVector)::Number
    """Compute log-likelihood functionfor Rust problem"""
    # Compute value
    Vbar = compute_Vbar(θ0, λ_, β, s)

    # Expected choice probabilities
    EP = exp.(Vbar[:,2]) ./ (exp.(Vbar[:,1]) + exp.(Vbar[:,2]))

    # Likelihood
    logL = sum(log.(EP[St[A.==1]])) + sum(log.(1 .- EP[St[A.==0]]))
    return -logL
end;

We can check the likelihood at the true value:

# True likelihood value
logL_trueθ = logL_Rust(θ, λ, β, s, St, A);
print("The likelihood at the true parameter is $logL_trueθ")

## The likelihood at the true parameter is 46371.98307795484

# Select starting values
θ0 = Float64[0,0,0];

# Optimize
θ_R = optimize(x -> logL_Rust(x, λ, β, s, St, A), θ0).minimizer;
print("Estimated thetas: $θ_R (true = $θ)")

## Estimated thetas: [0.12430277505205653, -0.0037182583936097658, 3.0416099489521438] (true = [0.13, -0.004, 3.1])

Starting values are important!

# Not all initial values are equally good
θ0 = Float64[1,1,1];

# Optimize
θ_R2 = optimize(x -> logL_Rust(x, λ, β, s, St, A), θ0).minimizer;
print("Estimated thetas: $θ_R2 (true = $θ)")

## Estimated thetas: [1.0, 1.0, 1.0] (true = [0.13, -0.004, 3.1])

Hotz & Miller estimation procedure works as follows

1. Estimate the CCPs from the data

2. Hotz & Miller inversion \[ \hat V = \Big[I - \beta \ \sum_a P_a .* T_a \Big]^{-1} \ * \ \left( \sum_a P_a \ .* \ \bigg[ u_a + \mathbb E [\epsilon_a] \bigg] \right) \]

3. Compute EP from EV \[ \hat \Pr(a=1 ; \theta) = \frac{\exp (u_1 +\beta T_1 \hat V )}{\sum_{a} \exp (u_a +\beta T_a \hat V )} \]

4. Compute the objective function: the (log)likelihood \[ \mathcal{L}(\theta) = \prod_{t=1}^{T}\left(\hat{\operatorname{Pr}}\left(a=1 \mid s_{t}; \theta\right) \mathbb{1}\left(a_{t}=1\right)+\left(1-\hat{\operatorname{Pr}}\left(a=0 \mid s_{t}; \theta\right)\right) \mathbb{1}\left(a_{t}=0\right)\right) \]

First, we need to estimate the Conditional Choice Proabilities (CCP)

• can be done non-parametrically
• i.e. just look at the frequency of investment in each state

# Estimate CCP
P = [mean(A[St.==i]) for i=s];
CCP = [(1 .- P) P]

## 10×2 Array{Float64,2}:
##  0.949131  0.0508692
##  0.92266   0.0773402
##  0.889882  0.110118
##  0.848416  0.151584
##  0.815519  0.184481
##  0.782453  0.217547
##  0.750602  0.249398
##  0.726972  0.273028
##  0.687072  0.312928
##  0.686552  0.313448

NeSr, we need \(T\), the matrices of transition probabilities, conditional on the investment choice.

function compute_T(k::Int, λ_::Number)::Array
    """Compute transition matrix"""
    T = zeros(k, k, 2);

    # Conditional on not investing
    T[k,k,1] = 1;
    for i=1:k-1
        T[i,i,1] = 1-λ_
        T[i,i+1,1] = λ_
    end

    # Conditional on investing
    T[:,1,2] .= 1-λ_;
    T[:,2,2] .= λ_;

    return(T)
end;

What form does the transition matrix \(T\) take?

# Compute T
T = compute_T(k, λ_);

# Conditional on not investing
T[:,:,1]

## 10×10 Array{Float64,2}:
##  0.180402  0.819598  0.0       0.0       …  0.0       0.0       0.0
##  0.0       0.180402  0.819598  0.0          0.0       0.0       0.0
##  0.0       0.0       0.180402  0.819598     0.0       0.0       0.0
##  0.0       0.0       0.0       0.180402     0.0       0.0       0.0
##  0.0       0.0       0.0       0.0          0.0       0.0       0.0
##  0.0       0.0       0.0       0.0       …  0.0       0.0       0.0
##  0.0       0.0       0.0       0.0          0.819598  0.0       0.0
##  0.0       0.0       0.0       0.0          0.180402  0.819598  0.0
##  0.0       0.0       0.0       0.0          0.0       0.180402  0.819598
##  0.0       0.0       0.0       0.0          0.0       0.0       1.0

Instead, the transitions conditional on investing are

# T Conditional on investing
T[:,:,2]

## 10×10 Array{Float64,2}:
##  0.180402  0.819598  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
##  0.180402  0.819598  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
##  0.180402  0.819598  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
##  0.180402  0.819598  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
##  0.180402  0.819598  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
##  0.180402  0.819598  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
##  0.180402  0.819598  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
##  0.180402  0.819598  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
##  0.180402  0.819598  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
##  0.180402  0.819598  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0

We now have all the pieces to compute the expected value function \(V\) through the Hotz & Miller inversion. \[ \hat V = \left[I - \beta \ \sum_a P_a .* T_a \right]^{-1} \ * \ \left( \sum_a P_a \ .* \ \bigg[ u_a + \mathbb E [\epsilon_a] \bigg] \right) \]

function HM_inversion(CCP::Matrix, T::Array, U::Matrix, β::Number)::Vector
    """Perform HM inversion"""

    # Compute LHS (to be inverted)
    γ = Base.MathConstants.eulergamma
    LEFT = I - β .* (CCP[:,1] .* T[:,:,1] + CCP[:,2] .* T[:,:,2])

    # Compute LHS (not to be inverted)
    RIGHT = γ .+ sum(CCP .* (U .- log.(CCP)) , dims=2)

    # Compute V
    EV_ = inv(LEFT) * RIGHT
    return vec(EV_)
end;

We can now compute the expected policy function from the expected value function \[ \hat \Pr(a=1 ; \theta) = \frac{\exp (u_1 +\beta T_1 \hat V )}{\sum_{a} \exp (u_a +\beta T_a \hat V )} \]

function from_EV_to_EP(EV_::Vector, T::Array, U::Matrix, β::Number)::Vector
    """Compute expected policy from expected value"""
    E = exp.( U + β .* [(T[:,:,1] * EV_) (T[:,:,2] * EV_)] )
    EP_ = E[:,2] ./ sum(E, dims=2)
    return vec(EP_)
end;

We now have all the pieces to build the likelihood function \[ \mathcal{L}(\theta) = \prod_{t=1}^{T} \left(\hat \Pr \left(a=1 \mid s_{t}; \theta\right) \mathbb{1} \left(a_{t}=1\right) + \left(1-\hat \Pr \left(a=0 \mid s_{t}; \theta\right)\right) \mathbb{1} \left(a_{t}=0\right)\right) \]

function logL_HM(θ0::Vector, λ::Number, β::Number, s::Vector, St::Vector, A::BitVector, T::Array, CCP::Matrix)::Number
    """Compute log-likelihood function for HM problem"""
    # Compute static utility
    U = compute_U(θ0, s)

    # Espected value by inversion
    EV_ = HM_inversion(CCP, T, U, β)

    # Implies choice probabilities
    EP_ = from_EV_to_EP(EV_, T, U, β)

    # Likelihood
    logL = sum(log.(EP_[St[A.==1]])) + sum(log.(1 .- EP_[St[A.==0]]))
    return -logL
end;

We can now estimate the parameters

# Optimize
θ0 = Float64[0,0,0];
θ_HM = optimize(x -> logL_HM(x, λ, β, s, St, A, T, CCP), θ0).minimizer;
print("Estimated thetas: $θ_HM (true = $θ)")

## Estimated thetas: [0.12420961622694673, -0.003710699952181827, 3.0413892912800002] (true = [0.13, -0.004, 3.1])

With Hotz and Miller, we have generated a mapping of the form

\[ \bar P(\cdot ; \theta) = g(h(\hat P(\cdot) ; \theta); \theta) \]

Aguirregabiria and Mira (2002): why don’t we iterate it?

The likelihood function in Aguirregabiria and Mira (2002) is extremely similar to Hotz and Miller (1993)

function logL_AM(θ0::Vector, λ::Number, β::Number, s::Vector, St::Vector, A::BitVector, T::Array, CCP::Matrix, K::Int)::Number
    """Compute log-likelihood function for AM problem"""
    # Compute static utility
    U = compute_U(θ0, s)
    EP_ = CCP[:,2]

    # Iterate HM mapping
    for _=1:K
        EV_ = HM_inversion(CCP, T, U, β)    # Expected value by inversion
        EP_ = from_EV_to_EP(EV_, T, U, β)   # Implies choice probabilities
        CCP = [(1 .- EP_) EP_]
    end

    # Likelihood
    logL = sum(log.(EP_[St[A.==1]])) + sum(log.(1 .- EP_[St[A.==0]]))
    return -logL
end;

We can now estimate the parameters

# Set number of iterations
K = 2;

# Optimize
θ0 = Float64[0,0,0];
θ_AM = optimize(x -> logL_AM(x, λ, β, s, St, A, T, CCP, K), θ0).minimizer;
print("Estimated thetas: $θ_AM (true = $θ)")

## Estimated thetas: [0.12430338308104717, -0.0037183107997343385, 3.0416115047879764] (true = [0.13, -0.004, 3.1])

Not much changes in our case.

We can compare the methods in terms of speed.

# Compare speed
θ0 = Float64[0,0,0];
optimize(x -> logL_Rust(x, λ, β, s, St, A), θ0).time_run

## 0.5985829830169678

optimize(x -> logL_HM(x, λ, β, s, St, A, T, CCP), θ0).time_run

## 0.31352901458740234

optimize(x -> logL_AM(x, λ, β, s, St, A, T, CCP, K), θ0).time_run

## 0.3311779499053955

Even in this simple example with a very small state space, the difference is significant.