20  Conditional Choice Probability Estimators

The nested fixed point algorithm from the DDC framework chapter requires us to solve the model every time we evaluate the likelihood. In our model of female labor force participation this takes a fraction of a second, but in many applications the state space is enormous, and solving the model even once can take minutes or hours. Realistic life-cycle models of labor supply, for example, might track experience in several occupations, a spouse’s employment, the ages of several children, and assets (Keane and Wolpin 1997; Keane, Todd, and Wolpin 2011). In dynamic games, the problem is compounded: solving the model means finding an equilibrium, which is harder than solving a single-agent problem and may not have a unique solution.

This chapter introduces a class of estimators that sidestep these problems by exploiting a simple fact: when all of the relevant state variables are observed, the conditional choice probabilities (CCPs) can be estimated directly from the data. The seminal insight of Hotz and Miller (1993) is that, armed with estimates of the CCPs, we can construct the continuation values that agents use to make decisions without solving the model. We’ll see how this insight leads to:

  1. Two-step estimators that require no model solution at all (Hotz and Miller 1993; Hotz et al. 1994).
  2. Recursive estimators that bridge the gap between two-step estimators and maximum likelihood (Aguirregabiria and Mira 2002).
  3. Estimators that exploit finite dependence to avoid computing continuation values over the entire state space (Altuğ and Miller 1998; Arcidiacono and Miller 2011).
  4. Tractable estimators for dynamic games (Aguirregabiria and Mira 2007; Pesendorfer and Schmidt-Dengler 2008; Bajari, Benkard, and Levin 2007).

20.1 Setup

We will use the model of female labor force participation as our running example. Recall that in each period \(t=1,\ldots,T\) a woman with experience \(h\), husband’s income \(y\), and young child indicator \(k\) chooses whether to work, with flow payoffs:

\[ u_{1}(h,y,k) = \log(y + w(h)),\qquad u_{0}(h,y,k) = \log(y) + \alpha_{0} + \alpha_{1}k \]

and extreme value shocks with scale \(\sigma\). The code for this model is collected in a script that we’ll load here. We’ll simulate the same dataset as in Example 18.2, estimate the wage equation, and calculate the NFXP estimates of \(\theta=(\alpha_{0},\alpha_{1},\sigma)\) as a benchmark. As before, we estimate \(\log(\sigma)\) and treat \(\beta\) as known.

include("../scripts/flfp_model.jl")
using Optim, ForwardDiff, Statistics

env = flfp_environment()
β = 0.95
θ_true = (;α0 = 0.3, α1 = 0.3, σ = 0.25)
x_true = [θ_true.α0, θ_true.α1, log(θ_true.σ)]
sol = solve_flfp(θ_true,env,β)
data = simulate_flfp(sol.logp1,env,2000; seed = 1)

# first step: the wage equation
env_hat = flfp_environment(γ = estimate_wages(data))

# benchmark: NFXP
x0 = [0.,0.,log(0.5)]
x_nfxp = optimize(x->-log_likelihood(x,data,env_hat,β),x0,Newton(),autodiff=:forward).minimizer
with_σ(x) = [x[1],x[2],exp(x[3])] #<- report σ rather than log(σ)
with_σ(x_nfxp)
3-element Vector{Float64}:
 0.2980938742851119
 0.31369294119805363
 0.2548586635371121

20.2 The Hotz-Miller Representation

Recall from our discussion of identification that, with extreme value shocks (and \(\sigma=1\)), the conditional value functions can be recovered from choice probabilities:

\[ v_{j}(x) - v_{0}(x) = \log p_{j}(x) - \log p_{0}(x). \]

This inversion holds much more generally. Hotz and Miller (1993) show that for any distribution of shocks with full support, there is a one-to-one mapping between choice probabilities and differences in conditional values. Moreover, for any choice \(j\), the integrated value function can be written as the conditional value of that choice plus a correction term that depends only on choice probabilities (Arcidiacono and Miller 2011):

\[ V(x) = v_{j}(x) + \psi_{j}(\mathbf{p}(x)). \]

In the logit case, \(\psi_{j}(\mathbf{p}) = \gamma - \log p_{j}\). Intuitively, \(\psi_{j}\) is the expected “option value” of being able to choose something other than \(j\) when the shocks are realized. If \(j\) is chosen with probability close to one, there is little option value and \(\psi_{j}\) is close to the mean of \(\varepsilon_{j}\) (\(\gamma\)).

An equally useful representation comes from averaging over choices. Since the agent chooses \(j\) with probability \(p_{j}(x)\), and the expected value of the shock conditional on choosing \(j\) is \(e_{j}(x) = \gamma - \log p_{j}(x)\) in the logit case, we can write:

\[ V(x) = \sum_{j}p_{j}(x)\left[u_{j}(x) + e_{j}(x) + \beta\sum_{x'}V(x')f_{j}(x'|x)\right]. \]

In a stationary, infinite-horizon model, we can write this in vector form. Letting \(\mathbf{F}^{\mathbf{p}} = \sum_{j}\text{diag}(\mathbf{p}_{j})\mathbf{F}_{j}\) as before:

\[ \mathbf{V} = \sum_{j}\mathbf{p}_{j}\circ\left(\mathbf{u}_{j} + \mathbf{e}_{j}\right) + \beta\mathbf{F}^{\mathbf{p}}\mathbf{V}\quad\Rightarrow\quad\mathbf{V}(\theta,\mathbf{p}) = \left(\mathbf{I}-\beta\mathbf{F}^{\mathbf{p}}\right)^{-1}\sum_{j}\mathbf{p}_{j}\circ\left(\mathbf{u}_{j}(\theta)+\mathbf{e}_{j}(\mathbf{p})\right). \tag{20.1}\]

This is just the value of following the (possibly suboptimal) policy encoded in \(\mathbf{p}\) forever. In a finite-horizon model like ours, choice probabilities depend on \(t\) and the same idea gives a backward recursion rather than a matrix inverse. With scale parameter \(\sigma\) and \(\mathbf{V}_{T+1}=\mathbf{0}\):

\[ \mathbf{V}_{t}(\theta,\mathbf{p}) = \sum_{j}\mathbf{p}_{jt}\circ\left(\mathbf{u}_{j}(\theta)+\sigma\mathbf{e}_{jt}(\mathbf{p})\right) + \beta\mathbf{F}^{\mathbf{p}}_{t}\mathbf{V}_{t+1}(\theta,\mathbf{p}),\qquad t=T,T-1,\ldots,1. \tag{20.2}\]

Notice that there is no maximization here: we are evaluating a given policy, not solving for the optimal one. Plugging \(\mathbf{V}\) into the conditional value functions and the logit formula gives us a mapping from payoff parameters and choice probabilities into choice probabilities:

\[ \Psi_{jt}(x;\theta,\mathbf{p}) = \frac{\exp\left(\tilde{v}_{jt}(x;\theta,\mathbf{p})/\sigma\right)}{\sum_{k}\exp\left(\tilde{v}_{kt}(x;\theta,\mathbf{p})/\sigma\right)},\qquad \tilde{\mathbf{v}}_{jt}(\theta,\mathbf{p}) = \mathbf{u}_{j}(\theta) + \beta\mathbf{F}_{j}\mathbf{V}_{t+1}(\theta,\mathbf{p}).\]

You can think of \(\Psi\) as a policy improvement step: given that the agent follows \(\mathbf{p}\) from tomorrow onwards, \(\Psi\) gives the optimal choice probabilities today. The key property is that the true choice probabilities are a fixed point of this mapping at the true parameters:

\[ \mathbf{p}_{0} = \Psi(\theta_{0},\mathbf{p}_{0}).\]

20.3 Two-Step Estimators

This fixed point property suggests a simple two-step strategy:

  1. Estimate the CCPs \(\hat{\mathbf{p}}\) (and the transitions \(\hat{\mathbf{F}}\)) directly from the data.
  2. Estimate \(\theta\) by maximizing the pseudo-likelihood: \[ \hat{\theta} = \arg\max_{\theta}\sum_{n=1}^{N}\sum_{t=1}^{T}\log\Psi_{D_{nt},t}(X_{nt};\theta,\hat{\mathbf{p}}).\]

No fixed point needs to be solved in the second step. In an infinite-horizon model, each trial value of \(\theta\) requires a single linear solve (and the matrix \(\mathbf{I}-\beta\mathbf{F}^{\hat{\mathbf{p}}}\) does not even depend on \(\theta\), so it can be factorized once). Better still, if flow payoffs are linear in parameters, \(\mathbf{u}_{j}(\theta)=\mathbf{Z}_{j}\theta\), then Equation 20.1 shows that \(\tilde{\mathbf{v}}_{j}\) is also linear in \(\theta\):

\[ \tilde{\mathbf{v}}_{j}(\theta,\hat{\mathbf{p}}) = \tilde{\mathbf{Z}}_{j}(\hat{\mathbf{p}})\theta + \tilde{\mathbf{e}}_{j}(\hat{\mathbf{p}}) \]

where \(\tilde{\mathbf{Z}}_{j}\) and \(\tilde{\mathbf{e}}_{j}\) depend only on the first stage. The same is true of the backward recursion in Equation 20.2. Once these objects are computed, the model never needs to be revisited: the second stage is just a logit with “offset” terms.

20.3.1 The First Stage

The two-step estimator can only be as good as its first stage. When the state space is small relative to the sample, we can use a frequency estimator: the fraction of observations in each state that make each choice. In our labor force participation model, the state in period \(t\) is \((h,y,k)\) with \(h\in\{0,\ldots,t-1\}\), so there are \(\sum_{t=1}^{30}10t=4,650\) feasible states in total. With 2,000 women observed for 30 years we have 60,000 observations, but they are spread very unevenly across these states: many states are never visited, and many more are visited only a handful of times. A frequency estimate in a cell with a few observations is very noisy, and an estimate of exactly zero or one would break the estimator, since the representation involves \(\log\hat{p}_{j}\).

The solution is to smooth the first stage. A common approach is a flexible logit: a logit model for the choice with a rich polynomial in the state variables and their interactions. This borrows strength across neighboring states, and always delivers probabilities strictly between zero and one. Other options include kernel and sieve estimators. The cost is that the first stage is now an approximation. If the approximation is poor, the second-stage estimates will be biased.

Example 20.1 First, let’s see how thinly spread the data are across the state space.

(;T,Ky) = env
counts = zeros(Int,T,Ky,2,T)
for n in axes(data.D,2), t in 1:T
    counts[data.H[t,n]+1,data.Y[t,n],data.K[t,n]+1,t] += 1
end
println("Feasible states: $(sum(10t for t in 1:T))")
println("States visited at least once: $(sum(counts .> 0))")
println("States visited at least 10 times: $(sum(counts .>= 10))")
Feasible states: 4650
States visited at least once: 3323
States visited at least 10 times: 1725

Fewer than half of the feasible states are visited at least ten times. So we’ll estimate the CCPs with a flexible logit. The function ccp_basis constructs a polynomial in \((t,h,\log y,k)\) with interactions. Since the logit log-likelihood is globally concave with a simple gradient and Hessian, we’ll write our own Newton iterations rather than calling Optim.

function ccp_basis(h,iy,k,t,env)
    ly = log(env.ygrid[iy]) - log(5.)
    τ = t/10; e = h/10 #<- rescale to keep the polynomial terms well-conditioned
    return [1., τ, τ^2, τ^3, e, e^2, e^3, τ*e, τ^2*e, τ*e^2, ly, ly*e, ly*τ, k, k*τ, k*e, k*ly]
end

# logit by Newton's method: gradient Z'(d-p), Hessian -Z'diag(p(1-p))Z
function logit_newton(Z,d; tol = 1e-10, max_iter = 50)
    b = zeros(size(Z,2))
    for iter in 1:max_iter
        p = 1 ./ (1 .+ exp.(-Z*b))
        g = Z' * (d .- p)
        H = Z' * (p .* (1 .- p) .* Z)
        b .+= H \ g
        if maximum(abs.(g)) < tol
            break
        end
    end
    return b
end

function flexible_logit_ccp(data,env)
    (;H,Y,K,D) = data
    (;T,Ky) = env
    Z = reduce(vcat,[ccp_basis(H[t,n],Y[t,n],K[t,n],t,env)' for n in axes(D,2) for t in 1:T])
    d = [D[t,n] for n in axes(D,2) for t in 1:T]
    b̂ = logit_newton(Z,d)
    logp1 = zeros(T,Ky,2,T)
    logp0 = zeros(T,Ky,2,T)
    for t in 1:T, k in 0:1, iy in 1:Ky, h in 0:t-1
        z = ccp_basis(h,iy,k,t,env)'*b̂
        logp1[h+1,iy,k+1,t] = -log(1+exp(-z))
        logp0[h+1,iy,k+1,t] = -log(1+exp(z))
    end
    return logp1,logp0
end

logp̂1,logp̂0 = flexible_logit_ccp(data,env)
# average absolute error in the first stage, weighting states by how often they are visited
sum(counts .* abs.(exp.(logp̂1) .- exp.(sol.logp1))) / sum(counts)
0.009796139843673295

On average, the estimated probabilities of working are within about one percentage point of the truth at the states we observe.

Now for the second stage. Notice that \(\theta=(\alpha_{0},\alpha_{1},\sigma)\) enters the flow payoffs and the expected shocks linearly, so by Equation 20.2 the value of following \(\hat{\mathbf{p}}\) is also linear in \(\theta\):

\[ \mathbf{V}_{t}(\theta,\hat{\mathbf{p}}) = \mathbf{V}^{y}_{t} + \alpha_{0}\mathbf{V}^{\alpha_{0}}_{t} + \alpha_{1}\mathbf{V}^{\alpha_{1}}_{t} + \sigma\mathbf{V}^{e}_{t} \]

where each component is computed by the backward recursion using only the corresponding piece of the flow payoff: \(\log\) income, the indicator for not working, the indicator for not working with a young child, and the expected shock \(e_{jt}\). The function hm_components below computes these components once, and stores the implied contribution of each one to the difference in conditional values \(\tilde{v}_{1t}-\tilde{v}_{0t}\). After that, evaluating \(\Psi\) at any \(\theta\) is just a weighted sum of four arrays.

# Value of following the policy p, decomposed as V = Vy + α0*Vα0 + α1*Vα1 + σ*Ve.
# Returns the contribution of each component to the difference in conditional values v1 - v0
function hm_components(logp1,logp0,env,β)
    (;T,Ky,ygrid,Fy,γ) = env
    V = zeros(4,T+1,Ky,2,T+1) #<- V[c,h+1,iy,k+1,t] for components c = (y,α0,α1,e)
    Δ = zeros(4,T,Ky,2,T)
    for t in T:-1:1
        Fk = child_transition(t,env)
        for k in 0:1, iy in 1:Ky, h in 0:t-1
            EV0 = zeros(4); EV1 = zeros(4)
            for k′ in 0:1, iy′ in 1:Ky
                f = Fy[iy,iy′]*Fk[k+1,k′+1]
                EV0 .+= f .* V[:,h+1,iy′,k′+1,t+1] #<- not working: experience stays at h
                EV1 .+= f .* V[:,h+2,iy′,k′+1,t+1] #<- working: experience rises to h+1
            end
            y = ygrid[iy]
            v1 = [log(y + (1+env.subsidy(t))*wage(h,γ)), 0., 0., 0.] .+ β .* EV1
            v0 = [log(y), 1., k, 0.] .+ β .* EV0
            l1 = logp1[h+1,iy,k+1,t]; l0 = logp0[h+1,iy,k+1,t]
            p1 = exp(l1); p0 = exp(l0)
            e = [0., 0., 0., p1*(γₑ - l1) + p0*(γₑ - l0)] #<- expected shock (per unit of σ)
            V[:,h+1,iy,k+1,t] .= p1 .* v1 .+ p0 .* v0 .+ e
            Δ[:,h+1,iy,k+1,t] .= v1 .- v0
        end
    end
    return Δ
end

# the mapping Ψ: log choice probabilities given θ, when p is followed in the future
function ccp_mapping(x,Δ)
    α0,α1,σ = x[1],x[2],exp(x[3])
    Δv = @views @. Δ[1,:,:,:,:] + α0*Δ[2,:,:,:,:] + α1*Δ[3,:,:,:,:] + σ*Δ[4,:,:,:,:]
    logp1 = @. -log(1 + exp(-Δv/σ))
    logp0 = @. -log(1 + exp(Δv/σ))
    return logp1,logp0
end
ccp_mapping (generic function with 1 method)

Before we estimate, let’s check that the true CCPs are a fixed point of \(\Psi\) at the true parameters. We only compare feasible states (those with \(h<t\)).

feasible = [h < t for h in 0:T-1, iy in 1:Ky, k in 0:1, t in 1:T]
Δ_true = hm_components(sol.logp1,sol.logp0,env,β)
logΨ1,_ = ccp_mapping(x_true,Δ_true)
maximum(abs.(logΨ1 .- sol.logp1)[feasible])
1.8489654252107357e-12

Great. The pseudo-likelihood is now just the likelihood of choices evaluated at the CCPs returned by \(\Psi\).

pseudo_loglikelihood(x,Δ,data) = choice_loglikelihood(ccp_mapping(x,Δ)...,data)

Δ̂ = hm_components(logp̂1,logp̂0,env_hat,β)
obj_2step(x) = -pseudo_loglikelihood(x,Δ̂,data)
x_2step = optimize(obj_2step,x0,Newton(),autodiff=:forward).minimizer
[with_σ(x_2step) with_σ(x_nfxp) with_σ(x_true)]
3×3 Matrix{Float64}:
 0.298492  0.298094  0.3
 0.313197  0.313693  0.3
 0.254493  0.254859  0.25

The two-step estimates (first column) are very close to the NFXP estimates (second column). Let’s compare the time it takes to run each estimator.

obj_nfxp(x) = -log_likelihood(x,data,env_hat,β)
optimize(obj_nfxp,x0,Newton(),autodiff=:forward); #<- run once to compile
optimize(obj_2step,x0,Newton(),autodiff=:forward);
println("NFXP:")
@time optimize(obj_nfxp,x0,Newton(),autodiff=:forward);
println("Two-step, first stage:")
@time begin
    logp̂1,logp̂0 = flexible_logit_ccp(data,env)
    Δ̂ = hm_components(logp̂1,logp̂0,env_hat,β)
end;
println("Two-step, second stage:")
@time optimize(obj_2step,x0,Newton(),autodiff=:forward);
NFXP:
  0.023886 seconds (7.60 k allocations: 48.270 MiB, 6.64% gc time)
Two-step, first stage:
  0.026233 seconds (380.84 k allocations: 135.073 MiB, 20.79% gc time)
Two-step, second stage:
  0.019772 seconds (574 allocations: 42.049 MiB, 20.62% gc time)

Interesting! The second stage is faster than NFXP, but not dramatically so, and once we include the first stage there is no gain at all. The reason is that, in a finite-horizon model with a few thousand states, solving the model is cheap: backward induction requires exactly one pass through the state space, which is the same work as evaluating Equation 20.2. The gains from CCP estimators become large when the model is expensive to solve: when the state space is very large, when a fixed point is hard to find (as in games), or when we can exploit finite dependence (see below) to avoid computing values over the whole state space.

20.3.2 Properties of the Two-Step Estimator

The two-step estimator is consistent and asymptotically normal as long as the first stage is consistent. In general, the asymptotic variance of a two-step estimator must account for estimation error in the first stage (see the theory of two-step estimators). However, Aguirregabiria and Mira (2002) point out a remarkable property of single-agent models: at the true parameters,

\[ \nabla_{\mathbf{p}}\Psi(\theta_{0},\mathbf{p}_{0}) = \mathbf{0}. \]

This zero Jacobian property is a consequence of optimality. Since \(\mathbf{p}_{0}\) is the optimal policy, small changes in the continuation policy have only second-order effects on continuation values (this is the envelope theorem at work). A consequence is that first-stage estimation error does not affect the asymptotic distribution of \(\hat{\theta}\): the two-step pseudo-likelihood estimator is asymptotically equivalent to the partial maximum likelihood estimator.

That’s good news, but it is an asymptotic result. In finite samples, the two-step estimator can perform poorly for two reasons:

  1. Noisy first stages. Frequency estimators are only feasible when the state space is small relative to the sample. In larger state spaces, the first stage must be smoothed (for example with a flexible logit or kernel estimator), and the quality of the second stage depends heavily on the quality of the first.
  2. Nonlinearity. The mapping \(\Psi\) is nonlinear in \(\mathbf{p}\) (it involves \(\log\mathbf{p}\)), so noise in \(\hat{\mathbf{p}}\) translates into finite-sample bias in \(\hat{\theta}\).

20.3.3 Nested Pseudo-Likelihood

Aguirregabiria and Mira (2002) proposed a way to reduce this finite-sample bias. The idea is to use the second-stage estimates to update the CCPs by applying \(\Psi\), and then re-estimate. Starting from \(\hat{\mathbf{p}}^{0}\), the nested pseudo-likelihood (NPL) algorithm iterates on:

  1. \(\hat{\theta}^{k} = \arg\max_{\theta}\sum_{n,t}\log\Psi_{D_{nt},t}(X_{nt};\theta,\hat{\mathbf{p}}^{k-1})\)
  2. \(\hat{\mathbf{p}}^{k} = \Psi(\hat{\theta}^{k},\hat{\mathbf{p}}^{k-1})\)

Each iteration replaces the CCPs from the first stage with CCPs that are consistent with the structure of the model. Aguirregabiria and Mira (2002) show that in single-agent models, if this algorithm converges, it converges to the maximum likelihood estimator. Notice that updating the CCPs with \(\Psi\) is one step of policy iteration, so the NPL algorithm “swaps” the order of the nested fixed point: instead of solving the model fully for each guess of \(\theta\), we take one solution step for each estimate of \(\theta\).

Example 20.2  

function npl(logp1,logp0,data,env,β; max_iter = 20, tol = 1e-8, verbose = true)
    (;T,Ky) = env
    feasible = [h < t for h in 0:T-1, iy in 1:Ky, k in 0:1, t in 1:T]
    x = [0.,0.,log(0.5)]
    for k in 1:max_iter
        Δ = hm_components(logp1,logp0,env,β)
        x = optimize(x->-pseudo_loglikelihood(x,Δ,data),x,Newton(),autodiff=:forward).minimizer
        l1,l0 = ccp_mapping(x,Δ)
        err = maximum(abs.(exp.(l1) .- exp.(logp1))[feasible])
        logp1,logp0 = l1,l0
        if verbose
            println("Iteration $k: θ = $(round.(with_σ(x),digits=4)), change in CCPs = $(round(err,sigdigits=3))")
        end
        if err<tol
            break
        end
    end
    return x
end

x_npl = npl(logp̂1,logp̂0,data,env_hat,β)
[with_σ(x_npl) with_σ(x_nfxp)]
Iteration 1: θ = [0.2985, 0.3132, 0.2545], change in CCPs = 0.121
Iteration 2: θ = [0.2981, 0.3137, 0.2549], change in CCPs = 0.00512
Iteration 3: θ = [0.2981, 0.3137, 0.2549], change in CCPs = 1.08e-5
Iteration 4: θ = [0.2981, 0.3137, 0.2549], change in CCPs = 2.31e-9
3×2 Matrix{Float64}:
 0.298094  0.298094
 0.313693  0.313693
 0.254859  0.254859

The algorithm converges in a few iterations, and the NPL estimates coincide with the NFXP estimates, just as Aguirregabiria and Mira (2002) predict. Notice that this is true even though our first stage was a (slightly misspecified) flexible logit: the limit of the NPL algorithm does not depend on where it starts.

20.4 Finite Dependence

The representations in Equation 20.1 and Equation 20.2 still require us to calculate values over the entire state space. In large models this can be just as burdensome as solving the model. Arcidiacono and Miller (2011) showed that in many models we can avoid this by exploiting finite dependence.

The idea is as follows. Recall that for any choice \(k\) we can write \(V(x') = v_{k}(x') + \psi_{k}(\mathbf{p}(x'))\). The difference in conditional values between choice \(j\) and choice \(0\) is:

\[ v_{j}(x) - v_{0}(x) = u_{j}(x) - u_{0}(x) + \beta\sum_{x'}\left[v_{k}(x') + \psi_{k}(\mathbf{p}(x'))\right]\left[f_{j}(x'|x)-f_{0}(x'|x)\right]. \]

Now suppose that there is a choice \(k\) (or more generally, a sequence of choices) such that, after choosing \(j\) today and \(k\) tomorrow, the distribution of the state the day after tomorrow is the same as after choosing \(0\) today and \(k\) tomorrow. Then the continuation values two periods ahead cancel out of the difference, and \(v_{j}(x)-v_{0}(x)\) depends only on flow payoffs and CCPs over a short horizon. When this happens after \(\rho\) periods, we say the model exhibits \(\rho\)-period finite dependence (Altuğ and Miller 1998; Arcidiacono and Miller 2011).

The simplest case is when there is a renewal or terminal action whose consequences do not depend on the past. Examples include replacing a bus engine (Rust 1987), retiring, or exiting a market. In the entry model, for example, choosing to be inactive resets the incumbency state regardless of whether the firm was active before, and this leads to one-period finite dependence (you will show this in Exercise 20.3).

Renewal actions are not the only source of finite dependence. In many labor models, the key state variable is accumulated experience, and no choice “resets” it. Nevertheless, Altuğ and Miller (1998) noticed that two different sequences of choices can lead to the same state. This is the case in our model of female labor force participation.

Example 20.3 Consider a woman in period \(t\) with experience \(h\). Compare two sequences of choices:

  1. Work in period \(t\) and don’t work in \(t+1\). Experience at \(t+2\) is \(h+1\).
  2. Don’t work in period \(t\) and work in \(t+1\). Experience at \(t+2\) is \(h+1\).

Since husband’s income and children evolve exogenously, both sequences lead to the same distribution of states in period \(t+2\). Hence the model exhibits two-period finite dependence. To use this, write the integrated value function in period \(t+1\) in terms of the conditional value of the choice prescribed by each sequence:

\[ V_{t+1}(h+1,y',k') = v_{0,t+1}(h+1,y',k') + \sigma\gamma - \sigma\log p_{0,t+1}(h+1,y',k') \] \[ V_{t+1}(h,y',k') = v_{1,t+1}(h,y',k') + \sigma\gamma - \sigma\log p_{1,t+1}(h,y',k'). \]

Both \(v_{0,t+1}(h+1,y',k')\) and \(v_{1,t+1}(h,y',k')\) lead to experience of \(h+1\) in period \(t+2\), so the continuation values in these terms are identical and cancel when we take the difference. This leaves:

\[ V_{t+1}(h+1,y',k') - V_{t+1}(h,y',k') = u_{0}(h+1,y',k') - u_{1}(h,y',k') + \sigma\left[\log p_{1,t+1}(h,y',k') - \log p_{0,t+1}(h+1,y',k')\right] \]

and hence, for \(t<T\):

\[ v_{1t}(h,y,k) - v_{0t}(h,y,k) = u_{1}(h,y,k) - u_{0}(h,y,k) + \beta\sum_{y',k'}f_{y}(y'|y)f_{kt}(k'|k)\left[V_{t+1}(h+1,y',k') - V_{t+1}(h,y',k')\right]. \]

In the final period, \(v_{1T}-v_{0T}=u_{1}-u_{0}\). Notice what this buys us: to evaluate the likelihood we only need flow payoffs and the CCPs one period ahead. We never need to solve the model.

The catch is that we need first-stage estimates of the CCPs one period ahead of every observed state, including states that are rarely (or never) visited. We’ll use the flexible logit estimates from Example 20.1. Now we can write the finite dependence log-likelihood:

function fd_loglikelihood(x,logp1,logp0,data,env,β)
    θ = (;α0 = x[1], α1 = x[2], σ = exp(x[3]))
    (;σ) = θ
    (;T,Ky,ygrid,Fy) = env
    (;H,Y,K,D) = data
    ll = zero(eltype(x))
    for n in axes(D,2), t in 1:T
        h = H[t,n]; iy = Y[t,n]; k = K[t,n]; y = ygrid[iy]
        Δv = u1(h,y,k,t,θ,env) - u0(h,y,k,t,θ,env)
        if t<T
            Fk = child_transition(t,env)
            for k′ in 0:1, iy′ in 1:Ky
                f = Fy[iy,iy′]*Fk[k+1,k′+1]
                y′ = ygrid[iy′]
                ΔV = u0(h+1,y′,k′,t+1,θ,env) - u1(h,y′,k′,t+1,θ,env) +
                     σ*(logp1[h+1,iy′,k′+1,t+1] - logp0[h+2,iy′,k′+1,t+1])
                Δv += β*f*ΔV
            end
        end
        ll += D[t,n]==1 ? -log(1+exp(-Δv/σ)) : -log(1+exp(Δv/σ))
    end
    return ll
end
fd_loglikelihood (generic function with 1 method)

As a check, if we plug in the true CCPs and the true parameters, the finite dependence likelihood should be identical to the full likelihood:

fd_loglikelihood(x_true,sol.logp1,sol.logp0,data,env,β) - log_likelihood(x_true,data,env,β)
5.820766091346741e-11

Finally, let’s estimate the model and compare to the two-step and NFXP estimates:

obj_fd(x) = -fd_loglikelihood(x,logp̂1,logp̂0,data,env_hat,β)
x_fd = optimize(obj_fd,x0,Newton(),autodiff=:forward).minimizer
[with_σ(x_fd) with_σ(x_2step) with_σ(x_nfxp) with_σ(x_true)]
3×4 Matrix{Float64}:
 0.284441  0.298492  0.298094  0.3
 0.30102   0.313197  0.313693  0.3
 0.249203  0.254493  0.254859  0.25

The finite dependence estimates (first column) are close to the NFXP estimates, but somewhat further away than the two-step estimates. Although it is cheap to compute, the finite dependence estimator uses the first-stage CCPs in a different (and less efficient) way: it relies on the estimated CCPs one period ahead directly, rather than on the value of following them for the rest of the life-cycle. Any error in the flexible logit therefore carries over more directly into the second stage. Let’s see how long it takes:

optimize(obj_fd,x0,Newton(),autodiff=:forward); #<- run once to compile
@time optimize(obj_fd,x0,Newton(),autodiff=:forward);
  0.554067 seconds (12.18 M allocations: 371.723 MiB, 3.28% gc time)

This is actually slower than the two-step estimator! The reason is that our implementation loops over every observation (and every possible state next period) each time the likelihood is evaluated, while the two-step estimator did the heavy lifting once, when it computed hm_components. We could speed it up by exploiting linearity in the same way. But the comparison misses the point. The cost of the finite dependence estimator depends only on the number of observations: it never needs to visit a state that is not observed in the data or one period ahead of an observed state. The cost of the two-step estimator (and of NFXP) grows with the size of the state space. In a model that tracks, say, the ages of several children and the employment history of a spouse, the state space can easily have billions of points and only finite dependence remains feasible.

Discussion

Finite dependence is a property of the model, not the data. Consider the following models and discuss whether they exhibit finite dependence:

  1. The labor supply model in Exercise 16.2, where human capital depreciates stochastically when not working.
  2. A version of the labor force participation model in which experience depreciates by one unit in every period that a woman does not work (with a floor of zero).
  3. A version of the labor force participation model in which a woman who did not work last period must pay a re-entry cost to work today.
  4. The model of welfare participation with time limits.

20.5 Dynamic Games

Recall the entry-exit model, where two firms make simultaneous entry decisions and each firm’s payoff depends on its rival’s decision. Estimating this model by NFXP faces two problems:

  1. Computation. For each trial value of \(\theta\), we must solve for a Markov perfect equilibrium, which involves finding a fixed point in both values and choice probabilities. Unlike the single-agent case, there is no guarantee that simple iteration on this mapping converges.
  2. Multiplicity. Dynamic games frequently have multiple equilibria. The likelihood is then a correspondence rather than a function of \(\theta\), and the econometrician must take a stand on which equilibrium is played in the data.

CCP estimators address both problems. Let \(p(x,a,a')\) be the equilibrium probability that a firm with incumbency status \(a\), facing a rival with status \(a'\), in a market with characteristic \(x\), enters. If the same equilibrium is played in all markets with the same observable characteristics, then \(p\) can be estimated directly from the data. Given an estimate \(\hat{p}\), each firm’s problem is a single-agent dynamic discrete choice problem: the rival’s behavior is just part of the (now known) transition process for the state. Thus, we can apply exactly the same logic as before to construct a best response mapping \(\Psi(\theta,\hat{p})\) and estimate \(\theta\) by pseudo-likelihood (Aguirregabiria and Mira 2007). The data effectively select the equilibrium for us.

Pesendorfer and Schmidt-Dengler (2008) showed that this idea can be cast as an asymptotic least squares problem:

\[ \hat{\theta} = \arg\min_{\theta}\left(\hat{\mathbf{p}} - \Psi(\theta,\hat{\mathbf{p}})\right)'\mathbf{W}\left(\hat{\mathbf{p}} - \Psi(\theta,\hat{\mathbf{p}})\right) \]

which is a special case of the minimum distance estimator. Different choices of \(\mathbf{W}\) yield different estimators in the literature, and there is an optimal choice. Bajari, Benkard, and Levin (2007) propose a related approach that uses forward simulation to compute continuation values, which is useful when states are continuous.

Example 20.4 We’ll load the code we wrote to solve the entry-exit model, which is collected in this script. Since that script uses some of the same function names as our labor force participation model (u0 and u1, for example), we load it inside a module and refer to its functions with the prefix EntryExit.. Following our discussion of identification in Example 17.4, we normalize the scrap value to zero.

module EntryExit
include("../models/entry_exit_model.jl")
end

ϕ_game = [1., 0.5, 1., 2., 0.] #<- (ϕ0,ϕ1,ϕ2,ϕ3,ϕ4)
xvals = collect(-1.:0.5:1.)
# equilibrium choice probabilities: p_eq[ix][a+1,a′+1]
p_eq = [EntryExit.calc_p(EntryExit.solve_model_newton(x,ϕ_game,β)) for x in xvals]
p_eq[3]
2×2 Matrix{Float64}:
 0.472387  0.343208
 0.886686  0.810081

Next we’ll simulate data from \(M=1,000\) markets over \(T=10\) periods. Each market is assigned one of the five values of \(x\).

function simulate_game(p_eq,M,T; seed = nothing)
    if !isnothing(seed)
        Random.seed!(seed)
    end
    Kx = length(p_eq)
    X = zeros(Int,M)
    A = zeros(Int,2,T,M) #<- A[f,t,m] is the incumbency status of firm f
    D = zeros(Int,2,T,M) #<- D[f,t,m] is the entry decision of firm f
    for m in 1:M
        ix = rand(1:Kx)
        X[m] = ix
        a = [0,0]
        for t in 1:T
            p = p_eq[ix]
            d = [rand() < p[a[1]+1,a[2]+1], rand() < p[a[2]+1,a[1]+1]]
            A[:,t,m] = a
            D[:,t,m] = d
            a = d
        end
    end
    return (;X,A,D)
end

data_game = simulate_game(p_eq,1000,10; seed = 2);

The first stage is a frequency estimator of \(p(x,a,a')\), pooling across both firms since the equilibrium is symmetric.

function frequency_ccp_game(data,Kx)
    (;X,A,D) = data
    num = zeros(Kx,2,2)
    den = zeros(Kx,2,2)
    for m in eachindex(X), t in axes(D,2), f in 1:2
        a = A[f,t,m]
        a′ = A[3-f,t,m]
        num[X[m],a+1,a′+1] += D[f,t,m]
        den[X[m],a+1,a′+1] += 1
    end
    return num ./ den
end

p̂_game = clamp.(frequency_ccp_game(data_game,length(xvals)),1e-3,1-1e-3)
[p̂_game[3,:,:] p_eq[3]]
2×4 Matrix{Float64}:
 0.484104  0.376812  0.472387  0.343208
 0.890646  0.8       0.886686  0.810081

Now the best response mapping. Fixing \(x\), the state for each firm is \((a,a')\). Given \(\hat{p}\), the rival enters with probability \(q = \hat{p}(x,a',a)\), which determines both the expected flow payoff from entry and the transition of the rival’s incumbency status. Otherwise, the code follows the same logic as Equation 20.1: since the game is stationary, we invert a matrix rather than iterating backwards in time.

function best_response(θ,p,x,β)
    ϕ0,ϕ1,ϕ2,ϕ3 = θ
    sidx = LinearIndices((2,2)) #<- s = sidx[a+1,a′+1]
    F0 = zeros(4,4)
    F1 = zeros(4,4)
    u1 = zeros(eltype(θ),4)
    u0 = zeros(4) #<- normalization
    for a in 0:1, a′ in 0:1
        s = sidx[a+1,a′+1]
        q = p[a′+1,a+1] #<- the rival's entry probability
        u1[s] = ϕ0 + ϕ1*x - ϕ2*q - ϕ3*(1-a)
        F0[s,sidx[1,2]] = q
        F0[s,sidx[1,1]] = 1-q
        F1[s,sidx[2,2]] = q
        F1[s,sidx[2,1]] = 1-q
    end
    P = p[:] #<- own choice probabilities, stacked in the same order
    FP = (1 .- P) .* F0 .+ P .* F1
    V = (I - β*FP) \ ((1 .- P) .* (u0 .+ γₑ .- log.(1 .- P)) .+ P .* (u1 .+ γₑ .- log.(P)))
    v0 = u0 .+ β * F0 * V
    v1 = u1 .+ β * F1 * V
    return reshape(1 ./ (1 .+ exp.(v0 .- v1)),2,2)
end

# check: the equilibrium is a fixed point of the best response at the true parameters
maximum(abs.(best_response(ϕ_game[1:4],p_eq[2],xvals[2],β) .- p_eq[2]))
6.661338147750939e-16

Finally, the pseudo-likelihood:

function pseudo_loglikelihood_game(θ,p̂,data,xvals,β)
    Ψ = [best_response(θ,p̂[ix,:,:],xvals[ix],β) for ix in eachindex(xvals)]
    (;X,A,D) = data
    ll = zero(eltype(θ))
    for m in eachindex(X), t in axes(D,2), f in 1:2
        P = Ψ[X[m]][A[f,t,m]+1,A[3-f,t,m]+1]
        ll += D[f,t,m]==1 ? log(P) : log(1-P)
    end
    return ll
end

res = optimize(θ->-pseudo_loglikelihood_game(θ,p̂_game,data_game,xvals,β),zeros(4),Newton(),autodiff=:forward)
[res.minimizer ϕ_game[1:4]]
4×2 Matrix{Float64}:
 0.974553  1.0
 0.49217   0.5
 0.965237  1.0
 1.9569    2.0

We recover the parameters of the game without ever solving for an equilibrium.

NPL in Games

It is tempting to apply the NPL algorithm to games to reduce finite-sample bias. Unfortunately, the zero Jacobian property does not hold in games (a rival’s choice probabilities are not chosen optimally by the firm), and Pesendorfer and Schmidt-Dengler (2010) show that NPL can converge to inconsistent estimates when the equilibrium in the data is not stable under best-response iteration. See Kasahara and Shimotsu (2012) for conditions under which NPL converges and Dearing and Blevins (2024) for an efficient and convergent alternative.

20.6 What Do CCP Estimators Assume?

The essential requirement for CCP estimators is that we can consistently estimate the choice probabilities that agents use when forming expectations about the future. This requires that:

  1. All persistent state variables are observed by the econometrician. If there is a persistent unobservable (say, a woman’s permanent taste for work), the CCPs that we estimate from the data are averages over the distribution of that unobservable among women in each observed state. These are not the CCPs that women use to form expectations, and the two-step estimator will be inconsistent.
  2. In games, the same equilibrium is played across markets (and over time) that share the same observable state.
Discussion

Suppose that women differ in a permanent, unobserved taste for work, so that \(\alpha_{0}\) varies across women, and that we ignore this when estimating the labor force participation model.

  1. Women who have worked a lot in the past (high \(h\)) will tend to be those with a low \(\alpha_{0}\). How will this affect the first-stage estimates of the CCPs, and the relationship between participation and experience that they display?
  2. Which parameters of the model do you expect to be biased, and in which direction? How does this problem relate to the classic distinction between state dependence and heterogeneity (Heckman 1981)?
  3. Could we fix the problem by including individual fixed effects in the first stage? What would you need for this to work?

These questions motivate the next part of the course, where we will see how to combine CCP estimators with the EM algorithm to handle unobserved heterogeneity (Arcidiacono and Miller 2011).

20.7 Exercises

Exercise

Exercise 20.1 Return to the labor supply model from Exercise 16.2 and assume that \(\beta\) and \(\pi\) are known.

  1. Write code to estimate the CCPs \(p(w,j_{t-1})\) by frequency estimation and use Equation 20.1 to implement the two-step pseudo-likelihood estimator of \((\kappa,\sigma)\). Hint: with scale parameter \(\sigma\), the expected shock conditional on choice \(j\) is \(\sigma(\gamma - \log p_{j})\) and the CCPs depend on \(\tilde{v}_{j}/\sigma\).
  2. Implement the NPL estimator.
  3. Run a Monte Carlo simulation to compare the two-step, NPL, and NFXP estimators for samples of \(N=100\) and \(N=1,000\) individuals observed for \(T=10\) periods. Comment on the finite sample performance of each estimator.
Exercise

Exercise 20.2 In this exercise you will implement a choice probability estimator for the model of labor supply and program participation with time limits. Use the parameters in default_pars() from the chapter (also available in the script), and the function solve to calculate the true choice probabilities.

  1. Use the function simulate from the chapter to simulate a panel of \(N=1000\) individuals over their life-cycle.
  2. For this sample of size \(N=1000\), write a function to estimate the array of choice probabilities using a frequency estimator. Note that in period \(t\), only experience levels up to \(t-1\) can be observed so you do not need to estimate this array in its entirety.
  3. Show that this model exhibits two-period finite dependence. Hint: compare the sequence “choose \(j\) today, then choose neither tomorrow” to the sequence “choose neither today, then choose \(j\) tomorrow”. Be careful with states in which the time limit binds.
  4. Write an estimation routine for the structural parameters \((\phi,\alpha,\varphi,\sigma,\gamma)\) that uses the first-stage CCP estimates and your finite dependence representation. Hint: you can estimate \(\phi\) directly using OLS.
  5. Can \(\beta\) be estimated using your approach? Try it.
  6. Write a Monte Carlo simulation in which you repeatedly draw a sample of size \(N=1000\) and produce an estimate. Use the results to check that your estimator is working properly, and compare its precision and speed to the NFXP estimator in Example 19.3.
Exercise

Exercise 20.3 Consider the single-agent entry model with the normalization \(\phi_{4}=0\), so that \(u_{0}(x,a)=0\).

  1. Show that choosing to be inactive (\(d=0\)) is a renewal action, and use this to show that: \[ v_{1}(x,a) - v_{0}(x,a) = \phi_{0} + \phi_{1}x - \phi_{3}(1-a) + \beta\sum_{x'}\mathbf{F}_{x}[x,x']\left[\log p_{0}(x',0) - \log p_{0}(x',1)\right]. \] Give an economic interpretation of the final term.
  2. Using the data from Example 17.2, estimate the CCPs by frequency and implement the finite dependence estimator, which is now just a binary logit with a known offset. Why might you need to trim the frequency estimates?
  3. Compare your estimates to the NFXP estimates from Example 17.3. Implement the two-step pseudo-likelihood estimator using Equation 20.1 and compare its estimates as well.
  4. What would go wrong with this argument if the scrap value \(\phi_{4}\) were not normalized to zero?
Exercise

Exercise 20.4 The club store data from Dearing and Blevins (2024) record whether each of three warehouse club chains is active in each of 1,610 county markets over 12 years, along with a categorical measure of market population (pop, taking values 1 through 5). In this exercise you will estimate a dynamic entry game using these data.

Let \(f\in\{1,2,3\}\) index chains and let \(N_{-f,t}\) be the number of rival chains active in period \(t\). Consider the flow payoffs: \[ u_{1f}(x,a_{f},N_{-f}) = \phi_{0f} + \phi_{1}\texttt{pop} - \phi_{2}N_{-f} - \phi_{3}(1-a_{f}),\qquad u_{0f} = 0 \] with type 1 extreme value shocks that are private information, and \(\beta=0.95\).

  1. What is the state space for each firm? How many states are there?
  2. Estimate the transition matrix for pop and the CCPs of each firm using frequency estimators. How many states are never (or rarely) visited? Propose and implement a way of smoothing the first-stage estimates.
  3. Implement a two-step pseudo-likelihood estimator of \((\phi_{01},\phi_{02},\phi_{03},\phi_{1},\phi_{2},\phi_{3})\).
  4. Calculate standard errors using a bootstrap that resamples markets.
  5. Use your estimates to compute the long-run (steady state) distribution of market structures when the entry cost \(\phi_{3}\) is reduced by half. Why must you solve for an equilibrium to answer this question, even though you did not need to solve one to estimate the model? What problems might you run into?