Coding: Rust (1987)

Setting

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

• 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

Actions

• 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} $$

Solving the Model

• 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

Code

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);

Static Utility

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;

Value Function

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;

Solving the Model

We can now solve for the value function.

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

DGP

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;

Generate the DAta

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 19207 investment decisions in 100000 observations

The Data

What does the data look like?

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

St	A	St1
3	FALSE	4
9	TRUE	1
8	TRUE	1
3	FALSE	4
9	FALSE	10
10	FALSE	10

Estimation - Lambda

• 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.8206570869594549 (true = 0.82)

Estimation - Theta

• 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

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 45937.866092460084

Estimating Theta

# 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.12063838656559037, -0.003220197034620527, 3.0865668144650487] (true = [0.13, -0.004, 3.1])

Starting Values

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

Recap

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) $$

CCPs

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.952419  0.0475814
##  0.923046  0.0769535
##  0.894443  0.105557
##  0.853306  0.146694
##  0.819293  0.180707
##  0.788935  0.211065
##  0.747248  0.252752
##  0.717915  0.282085
##  0.6947    0.3053
##  0.678452  0.321548

Transition Probabilities

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;

T

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.179343  0.820657  0.0       0.0       …  0.0       0.0       0.0
##  0.0       0.179343  0.820657  0.0          0.0       0.0       0.0
##  0.0       0.0       0.179343  0.820657     0.0       0.0       0.0
##  0.0       0.0       0.0       0.179343     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.820657  0.0       0.0
##  0.0       0.0       0.0       0.0          0.179343  0.820657  0.0
##  0.0       0.0       0.0       0.0          0.0       0.179343  0.820657
##  0.0       0.0       0.0       0.0          0.0       0.0       1.0

T (2)

Instead, the transitions conditional on investing are

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

## 10×10 Array{Float64,2}:
##  0.179343  0.820657  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
##  0.179343  0.820657  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
##  0.179343  0.820657  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
##  0.179343  0.820657  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
##  0.179343  0.820657  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
##  0.179343  0.820657  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
##  0.179343  0.820657  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
##  0.179343  0.820657  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
##  0.179343  0.820657  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
##  0.179343  0.820657  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0

Hotz & Miller Inversion

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;

From EV to EP

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;

Likelihood

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;

Estimation

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.12064911403839335, -0.003220614484856523, 3.086621855583483] (true = [0.13, -0.004, 3.1])

Aguirregabiria, Mira (2002)

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?

AM Likelihood Function

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;

Estimation

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.12063890836521114, -0.0032202282942220464, 3.086571461772538] (true = [0.13, -0.004, 3.1])

Not much changes in our case.

Speed

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.5477378368377686

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

## 0.3244161605834961

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

## 0.35499119758605957

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

Appendix

References [references]