21  Finite Mixtures and the EM Algorithm

In this chapter we introduce the workhorse approach to unobserved heterogeneity in structural models: assume that individuals belong to one of a finite number of types, and that the econometrician does not know which. This approach has a long history in econometrics. J. Heckman and Singer (1984) showed that the nonparametric maximum likelihood estimator of a mixing distribution in duration models is discrete, which provided early justification for approximating heterogeneity with a small number of types. We also saw in the discrete choice chapter that latent class models can approximate rich patterns of heterogeneity in tastes. And in practice, discrete types are computationally convenient since the model must only be solved once for each type.

We’ll begin by showing what can go wrong when we ignore heterogeneity, then introduce the EM algorithm, show how it combines with CCP estimators, and finally extend it to unobservables that change over time. Throughout, our running example is the model of female labor force participation.

21.1 Why Unobserved Heterogeneity Matters

One of the most robust facts about women’s labor supply is that it is persistent: a woman who worked last year is much more likely to work this year than one who did not. There are two very different explanations for this fact. The first is state dependence: working itself changes the payoffs of future work, for example because it builds human capital that raises future wages, or because re-entering the labor force is costly. The second is heterogeneity: some women simply have a stronger taste for work (or face lower costs of working) than others, and these differences persist over time. A woman who worked last year is more likely to be the kind of woman who likes to work.

This is the question that motivated J. J. Heckman (1981a), who studied exactly this problem using the participation decisions of married women, and it has been revisited many times since (e.g. Hyslop 1999). The distinction matters enormously for policy. If persistence reflects state dependence, then a policy that draws women into work today (say, a temporary wage subsidy) will have lasting effects on employment. If persistence reflects heterogeneity, the effects will disappear as soon as the policy ends.

Our model of female labor force participation has a source of state dependence built in: experience raises wages. Let’s suppose now that women also come in two types, \(\ell\in\{1,2\}\), that differ in their value of time at home, \(\alpha_{0\ell}\):

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

Women know their own type, but the econometrician does not. Let \(\pi_{\ell}\) be the fraction of women of type \(\ell\). (We use \(\ell\) to index types since \(k\) is already taken by the young child indicator.)

Example 21.1 We’ll load our code for the participation model and write a function to simulate women with types. We’ll assume that half of women have no particular distaste for work (\(\alpha_{0}=0\)) while the other half place a high value on time at home (\(\alpha_{0}=0.6\)). The other parameters are the same as in Example 18.1.

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

env = flfp_environment()
β = 0.95

# two types of women who differ in their value of time at home
θ_types = [(;α0 = 0.0, α1 = 0.3, σ = 0.25),
           (;α0 = 0.6, α1 = 0.3, σ = 0.25)]
π_true = [0.5,0.5]
sol_types = [solve_flfp(θ,env,β) for θ in θ_types]

function simulate_types(sol_types,πk,env,N; seed = nothing)
    if !isnothing(seed)
        Random.seed!(seed)
    end
    type = rand(Categorical(πk),N)
    T = env.T
    H = zeros(Int,T,N); Y = zeros(Int,T,N); K = zeros(Int,T,N); D = zeros(Int,T,N)
    logW = fill(NaN,T,N)
    # simulate each type separately and then stitch the panels together
    for ℓ in eachindex(πk)
        I = type.==ℓ
        dℓ = simulate_flfp(sol_types[ℓ].logp1,env,sum(I))
        H[:,I] .= dℓ.H; Y[:,I] .= dℓ.Y; K[:,I] .= dℓ.K; D[:,I] .= dℓ.D; logW[:,I] .= dℓ.logW
    end
    return (;H,Y,K,D,logW,type)
end

data = simulate_types(sol_types,π_true,env,2000; seed = 3);
for ℓ in 1:2
    println("Type $ℓ: participation rate = $(round(mean(data.D[:,data.type.==ℓ]),digits=3))")
end
Type 1: participation rate = 0.848
Type 2: participation rate = 0.26

The two types behave very differently: the first works most of the time, the second only occasionally. Now let’s estimate the model ignoring heterogeneity. Notice that the wage equation can still be estimated by OLS on workers. Since types only affect preferences, and the measurement error in wages does not affect choices, there is no selection problem in the wage equation.

γ̂ = estimate_wages(data)
env_hat = flfp_environment(γ = γ̂)
x_naive = optimize(x->-log_likelihood(x,data,env_hat,β),[0.,0.,log(0.5)],Newton(),autodiff=:forward).minimizer
println("Estimated (α0, α1, σ): $(round.([x_naive[1:2]; exp(x_naive[3])],digits=3))")
println("True α1 = 0.3, true σ = 0.25")
Estimated (α0, α1, σ): [0.425, 0.114, 0.157]
True α1 = 0.3, true σ = 0.25

The estimated scale of the taste shocks, \(\sigma\), is much too small, and the effect of a young child is badly attenuated. To see why, let’s compare the persistence of participation in the data to the persistence implied by the naive model.

# probability of working given work status last period
function persistence(D)
    p11 = mean(D[2:end,:][D[1:end-1,:].==1])
    p01 = mean(D[2:end,:][D[1:end-1,:].==0])
    return round.((p11,p01),digits=3)
end
sol_naive = solve_flfp((;α0 = x_naive[1], α1 = x_naive[2], σ = exp(x_naive[3])),env_hat,β)
data_naive = simulate_flfp(sol_naive.logp1,env_hat,20_000; seed = 5)
println("Data: P(work | worked) and P(work | didn't work) = $(persistence(data.D))")
println("Naive model: P(work | worked) and P(work | didn't work) = $(persistence(data_naive.D))")
Data: P(work | worked) and P(work | didn't work) = (0.742, 0.329)
Naive model: P(work | worked) and P(work | didn't work) = (0.638, 0.457)

The naive model’s only source of persistence is experience. To explain the strong persistence in the data, it must make choices very sensitive to the wage gains from experience, which it does by shrinking \(\sigma\). Even so, it cannot generate as much persistence as there is in the data. In other words, the naive model attributes heterogeneity to state dependence. The effect of children is attenuated for a related reason: much of the variation in participation is driven by the omitted types, and the model can only fit the average response.

Does this matter for policy? Let’s calculate the effect of a permanent 10% wage subsidy on the average participation rate over the life-cycle, using the true model and the naive estimates.

# average participation over the life-cycle for a mixture of types
function participation(θ_types,πk,env,β)
    return sum(πk[ℓ] * mean(life_cycle_profile(solve_flfp(θ_types[ℓ],env,β).logp1,env).part) for ℓ in eachindex(πk))
end
subsidy_effect(θ_types,πk,γ) = participation(θ_types,πk,flfp_environment(γ = γ, subsidy = t->0.1),β) - participation(θ_types,πk,flfp_environment(γ = γ),β)

θ_naive = [(;α0 = x_naive[1], α1 = x_naive[2], σ = exp(x_naive[3]))]
println("Effect of subsidy, true model: $(round(subsidy_effect(θ_types,π_true,env.γ),digits=3))")
println("Effect of subsidy, naive model: $(round(subsidy_effect(θ_naive,[1.],γ̂),digits=3))")
Effect of subsidy, true model: 0.028
Effect of subsidy, naive model: 0.078

The naive model overstates the response to the subsidy by a factor of almost three. Since it attributes all of the persistence in the data to human capital, it believes that participation is far more sensitive to wages than it really is.

21.2 Finite Mixture Models

Let’s set things up in some generality. Let \(\mathbf{w}_{n}=(X_{nt},D_{nt})_{t=1}^{T}\) collect the data for person \(n\), and let \(\theta_{\ell}\) be the parameters that govern behavior for type \(\ell\) (some elements of \(\theta_{\ell}\) may be common across types). The likelihood contribution of person \(n\), conditional on being type \(\ell\), is:

\[ \mathcal{L}_{n\ell}(\theta_{\ell}) = \prod_{t=1}^{T}p_{D_{nt}}(X_{nt};\theta_{\ell}).\]

Since the type is unobserved, the likelihood of the data for person \(n\) is a mixture:

\[ \mathcal{L}_{n}(\theta,\pi) = \sum_{\ell=1}^{L}\pi_{\ell}\mathcal{L}_{n\ell}(\theta_{\ell}). \]

If transitions of the observed states depended on type, we would include them in \(\mathcal{L}_{n\ell}\) too. For example, if types also differed in their wages, the likelihood of observed wages would enter \(\mathcal{L}_{n\ell}\), and wages would then help to reveal each woman’s type.

Initial Conditions

In our example, every woman is observed from the very beginning of her working life, at age 25 with no experience. Since the initial state is the same for everyone, it is independent of the type. Real data rarely look like this. In a panel such as the PSID, we typically start observing women part-way through their careers, and their experience at that point is the outcome of past participation decisions that depended on their type. A woman with a lot of experience at age 35 is more likely to be a type with a low value of time at home. Treating initial experience as exogenous then confounds the effect of experience with the effect of type. This is the initial conditions problem (J. J. Heckman 1981b). Common solutions are to let the type probabilities depend on the initial state, \(\pi_{\ell}(X_{n1})\), or to model the distribution of the initial state as an outcome of the model (for example by simulating each type’s history from the start of working life).

We could simply maximize \(\sum_{n}\log\mathcal{L}_{n}(\theta,\pi)\) directly. This is feasible in small problems, but the log of a sum is not separable across types, the likelihood surface tends to have many local maxima, and every evaluation requires solving the model for every type. The EM algorithm provides a more convenient approach.

21.3 The EM Algorithm

The Expectation-Maximization algorithm (Dempster, Laird, and Rubin 1977) is based on a simple inequality. For any probability distribution \(q_{n}=(q_{n1},\ldots,q_{nL})\) over types, Jensen’s inequality gives:

\[ \log\mathcal{L}_{n}(\theta,\pi) = \log\sum_{\ell}q_{n\ell}\frac{\pi_{\ell}\mathcal{L}_{n\ell}(\theta_{\ell})}{q_{n\ell}} \geq \sum_{\ell}q_{n\ell}\log\left(\frac{\pi_{\ell}\mathcal{L}_{n\ell}(\theta_{\ell})}{q_{n\ell}}\right). \]

The inequality holds with equality when \(q_{n\ell}\) is the posterior probability that person \(n\) is of type \(\ell\):

\[ q_{n\ell} = \frac{\pi_{\ell}\mathcal{L}_{n\ell}(\theta_{\ell})}{\sum_{\ell'}\pi_{\ell'}\mathcal{L}_{n\ell'}(\theta_{\ell'})}. \]

So the right-hand side defines a lower bound on the log-likelihood that touches it at the current parameter guess. The EM algorithm alternates between making the bound tight and maximizing the bound. Starting from a guess \((\theta^{0},\pi^{0})\), iterate:

  1. E-step: Calculate the posterior type probabilities \(q^{s}_{n\ell}\) using the formula above, evaluated at \((\theta^{s},\pi^{s})\).
  2. M-step: Maximize the expected complete-data log-likelihood: \[ (\theta^{s+1},\pi^{s+1}) = \arg\max_{\theta,\pi}\sum_{n}\sum_{\ell}q^{s}_{n\ell}\left[\log\pi_{\ell} + \log\mathcal{L}_{n\ell}(\theta_{\ell})\right]. \]

The M-step is much simpler than the original problem. The solution for the type probabilities is just the average posterior, \(\pi^{s+1}_{\ell} = \frac{1}{N}\sum_{n}q^{s}_{n\ell}\), and the problem for \(\theta\) is a weighted maximum likelihood problem, where each person contributes \(L\) observations, one for each type, weighted by \(q^{s}_{n\ell}\). If the parameters are all type-specific, this problem separates into \(L\) independent problems.

Theorem 21.1 (Monotonicity of EM (Dempster, Laird, and Rubin 1977)) Let \(\Lambda(\theta,\pi)=\sum_{n}\log\mathcal{L}_{n}(\theta,\pi)\). Each iteration of the EM algorithm weakly increases the log-likelihood: \(\Lambda(\theta^{s+1},\pi^{s+1})\geq\Lambda(\theta^{s},\pi^{s})\).

The proof follows directly from the inequality above. Writing \(Q(\theta,\pi;q)\) for the lower bound,

\[ \Lambda(\theta^{s+1},\pi^{s+1}) \geq Q(\theta^{s+1},\pi^{s+1};q^{s}) \geq Q(\theta^{s},\pi^{s};q^{s}) = \Lambda(\theta^{s},\pi^{s}).\]

The first inequality is Jensen’s, the second holds because the M-step maximizes \(Q\), and the equality holds because \(q^{s}\) are the posteriors at \((\theta^{s},\pi^{s})\).

A few practical points are worth emphasizing:

  1. Local maxima. EM is guaranteed to increase the likelihood, but not to find the global maximum. It is good practice to run the algorithm from several starting values.
  2. Slow convergence. EM typically converges linearly, and can be very slow when types are hard to distinguish. A common strategy is to use EM to get close to the solution and then switch to a gradient-based method on the full likelihood.
  3. Standard errors. Since EM converges to a maximum of the likelihood, standard errors should be computed from the full log-likelihood \(\Lambda(\theta,\pi)\) (for example using automatic differentiation), not from the weighted M-step problem, which treats the posteriors as known.
  4. Label switching. The likelihood is invariant to permuting the labels of the types, so types are only identified up to relabeling. This is harmless, but you must be careful when comparing estimates across samples (for example in a Monte Carlo or bootstrap).

Example 21.2 Let \(\theta = (\alpha_{01},\alpha_{02},\alpha_{1},\log\sigma)\). We’ll start with the E-step. Notice that we calculate the posteriors in logs and subtract the maximum before exponentiating. With \(T=30\) observations per woman, the likelihoods \(\mathcal{L}_{n\ell}\) can become very small, and working in levels risks numerical underflow.

# log-likelihood of the choices of woman n given log choice probabilities
function person_loglikelihood(logp1,logp0,data,n)
    (;H,Y,K,D) = data
    ll = zero(eltype(logp1))
    for t in axes(D,1)
        i = (H[t,n]+1,Y[t,n],K[t,n]+1,t)
        ll += D[t,n]==1 ? logp1[i...] : logp0[i...]
    end
    return ll
end

# ccps is a vector (one for each type) of named tuples with fields logp1 and logp0
# returns an L x N matrix of posterior type probabilities, and the log-likelihood
function e_step(ccps,πk,data)
    N = size(data.D,2)
    L = length(πk)
    q = zeros(L,N)
    ll = 0.
    for n in 1:N
        logL = [log(πk[ℓ]) + person_loglikelihood(ccps[ℓ].logp1,ccps[ℓ].logp0,data,n) for ℓ in 1:L]
        lmax = maximum(logL)
        q[:,n] .= exp.(logL .- lmax)
        total = sum(q[:,n])
        q[:,n] ./= total
        ll += lmax + log(total)
    end
    return q,ll
end
e_step (generic function with 1 method)

For the M-step, we need the weighted log-likelihood, which requires solving the model for each type.

unpack_types(x,L) = [(;α0 = x[ℓ], α1 = x[L+1], σ = exp(x[L+2])) for ℓ in 1:L]

function weighted_loglikelihood(x,q,data,env,β)
    L = size(q,1)
    ll = zero(eltype(x))
    for (ℓ,θ) in enumerate(unpack_types(x,L))
        (;logp1,logp0) = solve_flfp(θ,env,β)
        for n in axes(data.D,2)
            ll += q[ℓ,n] * person_loglikelihood(logp1,logp0,data,n)
        end
    end
    return ll
end

function em_full(x,πk,data,env,β; max_iter = 500, tol = 1e-6, verbose = false)
    L = length(πk)
    ll_old = -Inf
    for iter in 1:max_iter
        sols = [solve_flfp(θ,env,β) for θ in unpack_types(x,L)]
        q,ll = e_step(sols,πk,data) #<- E-step
        πk = mean(q,dims=2)[:] #<- M-step: type probabilities
        x = optimize(z->-weighted_loglikelihood(z,q,data,env,β),x,Newton(),autodiff=:forward).minimizer #<- M-step: preferences
        if verbose
            println("Iteration $iter, log-likelihood = $(round(ll,digits=3))")
        end
        if ll - ll_old < tol
            return (;x,πk,ll,iter)
        end
        ll_old = ll
    end
    return (;x,πk,ll = ll_old,iter = max_iter)
end
em_full (generic function with 1 method)

Let’s start from the naive estimates, perturbing the value of home time to break the symmetry between types.

x0 = [x_naive[1] - 0.3, x_naive[1] + 0.3, x_naive[2], x_naive[3]]
em = em_full(x0,[0.5,0.5],data,env_hat,β; verbose = true)
println("Estimates (α01, α02, α1, σ): $(round.([em.x[1:3]; exp(em.x[4])],digits=3))")
println("Type probabilities: $(round.(em.πk,digits=3))")
Iteration 1, log-likelihood = -33639.231
Iteration 2, log-likelihood = -28705.026
Iteration 3, log-likelihood = -28704.337
Iteration 4, log-likelihood = -28704.336
Iteration 5, log-likelihood = -28704.336
Estimates (α01, α02, α1, σ): [0.012, 0.596, 0.296, 0.246]
Type probabilities: [0.504, 0.496]

The EM algorithm recovers the true parameters in just a few iterations. Convergence is fast here because the two types behave so differently that the posteriors are close to zero or one for most women after a single E-step. What does the estimated model say about the wage subsidy?

θ_em = unpack_types(em.x,2)
println("Effect of subsidy, EM estimates: $(round(subsidy_effect(θ_em,em.πk,γ̂),digits=3))")
Effect of subsidy, EM estimates: 0.028

Once we account for heterogeneity, we get the policy effect right.

21.4 EM and CCP Estimators

The M-step in the example above requires solving the model for each type many times, which brings us back to the computational problem that CCP estimators were designed to solve. Arcidiacono and Miller (2011) showed how to combine the two approaches. The key insight is that we can also treat the type-specific CCPs, \(p_{j}(x,\ell)\), as parameters to be updated in the M-step. Since the M-step is a weighted likelihood problem, the nonparametric maximum likelihood estimate of the CCPs is a weighted frequency estimator:

\[ \hat{p}_{1}(x,\ell) = \frac{\sum_{n,t}q_{n\ell}D_{nt}\mathbf{1}\{X_{nt}=x\}}{\sum_{n,t}q_{n\ell}\mathbf{1}\{X_{nt}=x\}}. \]

When the state space is large, many cells will be thinly populated, and we can instead use any flexible model of the CCPs (such as a flexible logit) estimated by weighted maximum likelihood. This suggests a two-stage approach:

  1. Run the EM algorithm using flexible, type-specific CCPs. This delivers estimates of the CCPs \(\hat{p}(x,\ell)\), the type probabilities, and the posteriors \(q_{n\ell}\), all without imposing any structure from the economic model.
  2. Estimate the structural parameters using any of the CCP estimators, treating the type as known but weighting each person-type combination by \(q_{n\ell}\).

Arcidiacono and Miller (2011) also propose a version that imposes the structure of the model in the first stage, by updating the CCPs in each iteration with the mapping \(\Psi\) instead of the flexible estimator. This is more efficient and is useful when the state space is too large for flexible CCPs to be precisely estimated.

Example 21.3 The participation model has 4,650 states, too many for a frequency estimator of the CCPs. Instead, we’ll use a flexible logit for each type: a logit for \(d_{t}\) with a polynomial in age, experience, log husband’s income, and the child indicator. The E-step is unchanged. For the M-step, we need a weighted logit. Since the logit log-likelihood is globally concave with simple derivatives, we’ll write our own Newton routine:

\[ \nabla = \sum_{i}w_{i}(d_{i}-p_{i})\mathbf{z}_{i},\qquad \nabla^{2} = -\sum_{i}w_{i}p_{i}(1-p_{i})\mathbf{z}_{i}\mathbf{z}_{i}'. \]

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

# stack the data into a design matrix with one row per person-period (t varies fastest)
function design(data,env)
    (;H,Y,K,D) = data
    T,N = size(D)
    Z = reduce(vcat,[ccp_basis(H[t,n],Y[t,n],K[t,n],t,env)' for n in 1:N for t in 1:T])
    d = [D[t,n] for n in 1:N for t in 1:T]
    return Z,d
end

# weighted logit by Newton's method, starting from b
function weighted_logit(Z,d,w,b; tol = 1e-8, max_iter = 50)
    for iter in 1:max_iter
        p = 1 ./ (1 .+ exp.(-Z*b))
        g = Z' * (w .* (d .- p))
        Hs = Z' * (Z .* (w .* p .* (1 .- p)))
        step = Hs \ g
        b = b .+ step
        if maximum(abs.(step)) < tol
            break
        end
    end
    return b
end

# log choice probabilities over the whole state space implied by logit coefficients b
function ccp_arrays(b,env)
    (;T,Ky) = env
    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] = -log1p(exp(-z))
        logp0[h+1,iy,k+1,t] = -log1p(exp(z))
    end
    return (;logp1,logp0)
end

function em_ccp(b_types,πk,data,env,Z,d; max_iter = 1000, tol = 1e-6)
    T = size(data.D,1)
    L = length(πk)
    ll_old = -Inf
    for iter in 1:max_iter
        ccps = [ccp_arrays(b,env) for b in b_types]
        q,ll = e_step(ccps,πk,data)
        πk = mean(q,dims=2)[:]
        # each person's posterior weight applies to all T of their observations
        b_types = [weighted_logit(Z,d,repeat(q[ℓ,:],inner=T),b_types[ℓ]) for ℓ in 1:L]
        if ll - ll_old < tol
            return (;b_types,πk,q,ll,iter)
        end
        ll_old = ll
    end
    ccps = [ccp_arrays(b,env) for b in b_types]
    q,ll = e_step(ccps,πk,data)
    return (;b_types,πk,q,ll,iter = max_iter)
end
em_ccp (generic function with 1 method)

We’ll initialize by estimating a pooled logit and shifting the intercept up for the first type and down for the second.

Z,d = design(data,env)
b_pool = weighted_logit(Z,d,ones(length(d)),zeros(size(Z,2)))
shift(b,c) = b .+ [c; zeros(length(b)-1)]
em1 = em_ccp([shift(b_pool,0.5),shift(b_pool,-0.5)],[0.5,0.5],data,env,Z,d)
println("Converged in $(em1.iter) iterations. Type probabilities: $(round.(em1.πk,digits=3))")
ccps_hat = [ccp_arrays(b,env) for b in em1.b_types]
for ℓ in 1:2
    est = life_cycle_profile(ccps_hat[ℓ].logp1,env).part
    truth = life_cycle_profile(sol_types[ℓ].logp1,env).part
    println("Type $ℓ participation at ages 25, 35, 45, 54. Estimated: $(round.(est[[1,11,21,30]],digits=2)). True: $(round.(truth[[1,11,21,30]],digits=2))")
end
Converged in 18 iterations. Type probabilities: [0.505, 0.495]
Type 1 participation at ages 25, 35, 45, 54. Estimated: [0.81, 0.82, 0.88, 0.92]. True: [0.8, 0.83, 0.88, 0.91]
Type 2 participation at ages 25, 35, 45, 54. Estimated: [0.21, 0.22, 0.29, 0.33]. True: [0.22, 0.22, 0.28, 0.33]

The type-specific CCPs are estimated well, without imposing any structure from the model. For the second stage, we’ll use the two-period finite dependence representation from Example 20.3: a woman who works today and not tomorrow arrives at the same experience level as a woman who doesn’t work today and works tomorrow. We now allow for type-specific values of home time, use type-specific CCPs, and weight each person-type combination by its posterior probability.

function fd_loglikelihood_types(x,ccps,q,data,env,β)
    L = length(ccps)
    (;T,Ky,ygrid,Fy) = env
    (;H,Y,K,D) = data
    σ = exp(x[L+2])
    ll = zero(eltype(x))
    for ℓ in 1:L
        θ = (;α0 = x[ℓ], α1 = x[L+1], σ)
        (;logp1,logp0) = ccps[ℓ]
        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
                    y′ = ygrid[iy′]
                    # V(h+1,⋅) - V(h,⋅) next period, using the "don't work" and "work" choices respectively
                    Δ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 += β * Fy[iy,iy′] * Fk[k+1,k′+1] * ΔV
                end
            end
            ll += q[ℓ,n] * (D[t,n]==1 ? -log1p(exp(-Δv/σ)) : -log1p(exp(Δv/σ)))
        end
    end
    return ll
end

x_ccp = optimize(x->-fd_loglikelihood_types(x,ccps_hat,em1.q,data,env_hat,β),[0.,0.5,0.,log(0.5)],Newton(),autodiff=:forward).minimizer
println("Estimates (α01, α02, α1, σ)")
println("  CCP-EM: $(round.([x_ccp[1:3]; exp(x_ccp[4])],digits=3))")
println("  Full EM: $(round.([em.x[1:3]; exp(em.x[4])],digits=3))")
println("  Truth: $([0.,0.6,0.3,0.25])")
Estimates (α01, α02, α1, σ)
  CCP-EM: [-0.014, 0.566, 0.271, 0.228]
  Full EM: [0.012, 0.596, 0.296, 0.246]
  Truth: [0.0, 0.6, 0.3, 0.25]

The two-stage CCP estimates are close to the full EM estimates and to the truth, although the estimate of \(\sigma\) is somewhat too low. As we saw in the CCP chapter, the price of not solving the model is that the second stage inherits the noise (and any misspecification) of the flexible first stage. The advantage is that we never solved the model: the second stage is just a weighted logit with an offset.

21.4.1 Choosing the Number of Types

How many types should we use? Since adding types can only increase the likelihood, we need a criterion that penalizes complexity. A common choice is the Bayesian Information Criterion:

\[ \text{BIC}(L) = -2\Lambda(\hat{\theta}_{L},\hat{\pi}_{L}) + n_{L}\log(N) \]

where \(n_{L}\) is the number of free parameters and \(N\) is the number of independent units (here, women). We choose the \(L\) that minimizes BIC.

Example 21.4 Let’s use the flexible first stage from Example 21.3, where each type has its own set of 17 logit coefficients.

N = size(data.D,2)
for L in 1:3
    b_init = [shift(b_pool,0.6*(ℓ - (L+1)/2)) for ℓ in 1:L]
    res = em_ccp(b_init,fill(1/L,L),data,env,Z,d)
    nL = L*size(Z,2) + L - 1
    println("L = $L: log-likelihood = $(round(res.ll,digits=1)), BIC = $(round(-2res.ll + nL*log(N),digits=1)), type probabilities = $(round.(res.πk,digits=3))")
end
L = 1: log-likelihood = -29534.6, BIC = 59198.4, type probabilities = [1.0]
L = 2: log-likelihood = -28699.2, BIC = 57664.4, type probabilities = [0.495, 0.505]
L = 3: log-likelihood = -28681.8, BIC = 57766.4, type probabilities = [0.489, 0.008, 0.504]

The BIC correctly selects two types. With three types, the likelihood improves only slightly, not enough to justify 18 additional parameters. Notice also that the third type is given a tiny probability: the algorithm is essentially telling us that it has nothing to do.

Discussion

The BIC penalizes model complexity using the number of parameters in the first stage in this example. In applications with structural models, the number of types often matters most through how it changes the counterfactual predictions of the model.

  1. Suppose two specifications (with \(L=2\) and \(L=3\)) have similar BIC but different predictions for the effect of a wage subsidy. How would you decide between them?
  2. Is it reasonable to think that heterogeneity in women’s taste for work is truly discrete? If not, how would you interpret the types?
  3. In our example, types differ only in preferences. What additional data would help to distinguish types if they also differed in their earnings potential?

21.5 Time-Varying Unobserved States

So far we have assumed that the unobserved type is fixed over time. In many applications the unobservables evolve: a woman’s health may deteriorate or improve, demands on her time at home may rise when an elderly parent needs care, a worker’s human capital may appreciate or depreciate. We can model this by letting the unobserved state \(\ell_{t}\) follow a Markov chain with transition matrix \(\Pi\), where \(\Pi[\ell,\ell'] = \mathbb{P}(\ell_{t+1}=\ell'|\ell_{t}=\ell)\), and initial distribution \(\mu\). This is a hidden Markov model (HMM). Notice that permanent types are the special case where \(\Pi=\mathbf{I}\).

The difficulty is that the likelihood now requires summing over all possible paths of the unobserved state:

\[ \mathcal{L}_{n} = \sum_{\ell_{1},\ldots,\ell_{T}}\mu(\ell_{1})g_{n1}(\ell_{1})\prod_{t=2}^{T}\Pi[\ell_{t-1},\ell_{t}]g_{nt}(\ell_{t}) \]

where \(g_{nt}(\ell) = p_{D_{nt}}(X_{nt},\ell)\) is the probability of the observed choice given the unobserved state. There are \(L^{T}\) such paths, so direct calculation quickly becomes infeasible. Fortunately, the Markov structure allows us to calculate this sum recursively.

21.5.1 The Forward-Backward Algorithm

Suppressing the individual subscript, define the forward probabilities as the joint probability of the unobserved state today and the observed choices up to today:

\[ \alpha_{t}(\ell) = \mathbb{P}(\ell_{t}=\ell, D_{1},\ldots,D_{t}) \]

and the backward probabilities as the probability of future choices given the state today:

\[ \beta_{t}(\ell) = \mathbb{P}(D_{t+1},\ldots,D_{T}|\ell_{t}=\ell) \]

where all probabilities are implicitly conditioned on the observed states. These satisfy simple recursions:

\[ \alpha_{1}(\ell) = \mu(\ell)g_{1}(\ell),\qquad \alpha_{t}(\ell) = \left[\sum_{\ell'}\alpha_{t-1}(\ell')\Pi[\ell',\ell]\right]g_{t}(\ell) \] \[ \beta_{T}(\ell) = 1,\qquad \beta_{t}(\ell) = \sum_{\ell'}\Pi[\ell,\ell']g_{t+1}(\ell')\beta_{t+1}(\ell'). \]

The likelihood is \(\mathcal{L}=\sum_{\ell}\alpha_{T}(\ell)\), and for the E-step we need two sets of posterior probabilities:

\[ q_{t}(\ell) = \mathbb{P}(\ell_{t}=\ell|D_{1},\ldots,D_{T}) = \frac{\alpha_{t}(\ell)\beta_{t}(\ell)}{\mathcal{L}} \] \[ \xi_{t}(\ell,\ell') = \mathbb{P}(\ell_{t}=\ell,\ell_{t+1}=\ell'|D_{1},\ldots,D_{T}) = \frac{\alpha_{t}(\ell)\Pi[\ell,\ell']g_{t+1}(\ell')\beta_{t+1}(\ell')}{\mathcal{L}}. \]

The M-step then updates each piece of the model using weighted frequency estimators (Rabiner 1989):

\[ \hat{\Pi}[\ell,\ell'] = \frac{\sum_{n}\sum_{t<T}\xi_{nt}(\ell,\ell')}{\sum_{n}\sum_{t<T}q_{nt}(\ell)},\qquad \hat{\mu}(\ell) = \frac{1}{N}\sum_{n}q_{n1}(\ell) \]

and the CCPs are updated as before, using weights \(q_{nt}(\ell)\) that now vary over time.

Numerical Stability

Since \(\alpha_{t}\) is a joint probability of \(t\) events, it shrinks geometrically with \(t\) and will underflow for long panels. The standard fix is to normalize \(\alpha_{t}\) to sum to one in each period, keeping track of the normalizing constants \(c_{t}\). The log-likelihood is then \(\sum_{t}\log c_{t}\). The backward probabilities can be normalized in the same way since any rescaling of \(\beta_{t}\) cancels out of the posteriors.

Example 21.5 Suppose now that, instead of being a permanent type, a woman’s value of time at home switches between the two values from Example 21.1. Think of the high value (\(\alpha_{0}=0.6\)) as periods in which demands at home are high, for example because of a health problem or the need to care for a family member. The regime follows a Markov chain with transition matrix:

\[ \Pi = \left[\begin{array}{cc}0.9 & 0.1 \\ 0.2 & 0.8\end{array}\right] \]

and each regime is equally likely at age 25. We’ll assume, for now, that women do not anticipate changes in the regime, so that the regime-specific CCPs are the same as the type-specific CCPs from Example 21.1. (In Exercise 21.3 you will relax this.)

Π_true = [0.9 0.1; 0.2 0.8]
μ_true = [0.5,0.5]

function simulate_regimes(sol_types,Π,μ,env,N; σξ = 0.3, seed = nothing)
    if !isnothing(seed)
        Random.seed!(seed)
    end
    (;T,Fy,γ) = env
    π∞ = stationary(Fy)
    H = zeros(Int,T,N); Y = zeros(Int,T,N); K = zeros(Int,T,N); D = zeros(Int,T,N)
    R = zeros(Int,T,N) #<- the unobserved regime
    logW = fill(NaN,T,N)
    for n in 1:N
        h = 0
        iy = rand(Categorical(π∞))
        k = rand() < 0.3
        r = rand(Categorical(μ))
        for t in 1:T
            d = rand() < exp(sol_types[r].logp1[h+1,iy,k+1,t])
            H[t,n] = h; Y[t,n] = iy; K[t,n] = k; D[t,n] = d; R[t,n] = r
            if d==1
                logW[t,n] = log(wage(h,γ)) + σξ*randn()
            end
            h += d
            iy = rand(Categorical(Fy[iy,:]))
            k = rand() < child_transition(t,env)[k+1,2]
            r = rand(Categorical(Π[r,:]))
        end
    end
    return (;H,Y,K,D,logW,R)
end

data_hmm = simulate_regimes(sol_types,Π_true,μ_true,env,2000; seed = 4);

Here is the forward-backward algorithm for a single woman. It takes an \(L\times T\) matrix \(g\) of choice probabilities in each unobserved state, fills in the posteriors \(q\), adds the pairwise posteriors to \(\xi\), and returns the log-likelihood. We write out the loops explicitly and use the @views macro (so that slices like α[:,t] do not allocate copies) since this function will be called many times.

@views function forward_back!(q,ξ,α,βb,g,Π,μ)
    L,T = size(g)
    loglik = 0.
    # forward pass
    for t in 1:T
        c = 0.
        for ℓ in 1:L
            α[ℓ,t] = t==1 ? μ[ℓ]*g[ℓ,1] : sum(α[ℓ′,t-1]*Π[ℓ′,ℓ] for ℓ′ in 1:L) * g[ℓ,t]
            c += α[ℓ,t]
        end
        α[:,t] ./= c #<- normalize to avoid underflow
        loglik += log(c)
    end
    # backward pass
    βb[:,T] .= 1.
    for t in reverse(1:T-1)
        for ℓ in 1:L
            βb[ℓ,t] = sum(Π[ℓ,ℓ′]*g[ℓ′,t+1]*βb[ℓ′,t+1] for ℓ′ in 1:L)
        end
        βb[:,t] ./= sum(βb[:,t])
    end
    # posteriors
    for t in 1:T
        for ℓ in 1:L
            q[ℓ,t] = α[ℓ,t]*βb[ℓ,t]
        end
        q[:,t] ./= sum(q[:,t])
    end
    # pairwise posteriors
    for t in 1:T-1
        c = 0.
        for ℓ in 1:L, ℓ′ in 1:L
            c += α[ℓ,t]*Π[ℓ,ℓ′]*g[ℓ′,t+1]*βb[ℓ′,t+1]
        end
        for ℓ in 1:L, ℓ′ in 1:L
            ξ[ℓ,ℓ′] += α[ℓ,t]*Π[ℓ,ℓ′]*g[ℓ′,t+1]*βb[ℓ′,t+1] / c
        end
    end
    return loglik
end
forward_back! (generic function with 1 method)

The E-step loops over women, and the M-step updates the flexible logit CCPs (using the weighted logit from Example 21.3 with weights that now vary by period), the transition matrix, and the initial distribution.

function e_step_hmm!(Q,ξ,workspace,ccps,Π,μ,data)
    (;H,Y,K,D) = data
    (;α,βb,g) = workspace
    T,N = size(D)
    L = length(ccps)
    fill!(ξ,0.)
    ll = 0.
    for n in 1:N
        for t in 1:T, ℓ in 1:L
            i = (H[t,n]+1,Y[t,n],K[t,n]+1,t)
            g[ℓ,t] = D[t,n]==1 ? exp(ccps[ℓ].logp1[i...]) : exp(ccps[ℓ].logp0[i...])
        end
        ll += @views forward_back!(Q[:,:,n],ξ,α,βb,g,Π,μ)
    end
    return ll
end

function em_hmm(b_types,Π,μ,data,env,Z,d; max_iter = 2000, tol = 1e-6)
    T,N = size(data.D)
    L = length(b_types)
    Q = zeros(L,T,N)
    ξ = zeros(L,L)
    workspace = (;α = zeros(L,T), βb = zeros(L,T), g = zeros(L,T))
    ll_old = -Inf
    for iter in 1:max_iter
        ccps = [ccp_arrays(b,env) for b in b_types]
        ll = e_step_hmm!(Q,ξ,workspace,ccps,Π,μ,data) #<- E-step
        # M-step. Q[ℓ,:,:][:] stacks the weights with t varying fastest, matching the rows of Z
        b_types = [weighted_logit(Z,d,Q[ℓ,:,:][:],b_types[ℓ]) for ℓ in 1:L]
        Π = ξ ./ sum(ξ,dims=2)
        μ = mean(Q[:,1,:],dims=2)[:]
        if ll - ll_old < tol
            return (;b_types,Π,μ,Q,ll,iter)
        end
        ll_old = ll
    end
    return (;b_types,Π,μ,Q,ll = ll_old,iter = max_iter)
end
em_hmm (generic function with 1 method)

Let’s run it.

Z_hmm,d_hmm = design(data_hmm,env)
b_pool_hmm = weighted_logit(Z_hmm,d_hmm,ones(length(d_hmm)),zeros(size(Z_hmm,2)))
b_init = [shift(b_pool_hmm,0.5),shift(b_pool_hmm,-0.5)]
hmm = em_hmm(b_init,[0.8 0.2; 0.2 0.8],[0.5,0.5],data_hmm,env,Z_hmm,d_hmm)
println("Converged in $(hmm.iter) iterations")
println("Estimated transition matrix: $(round.(hmm.Π,digits=3))")
println("Estimated initial distribution: $(round.(hmm.μ,digits=3))")
ccps_hmm = [ccp_arrays(b,env) for b in hmm.b_types]
for ℓ in 1:2
    est = life_cycle_profile(ccps_hmm[ℓ].logp1,env).part
    truth = life_cycle_profile(sol_types[ℓ].logp1,env).part
    println("Regime $ℓ participation at ages 25, 35, 45, 54. Estimated: $(round.(est[[1,11,21,30]],digits=2)). True: $(round.(truth[[1,11,21,30]],digits=2))")
end
Converged in 241 iterations
Estimated transition matrix: [0.9 0.1; 0.209 0.791]
Estimated initial distribution: [0.572, 0.428]
Regime 1 participation at ages 25, 35, 45, 54. Estimated: [0.78, 0.81, 0.89, 0.92]. True: [0.8, 0.83, 0.88, 0.91]
Regime 2 participation at ages 25, 35, 45, 54. Estimated: [0.16, 0.24, 0.29, 0.31]. True: [0.22, 0.22, 0.28, 0.33]

To summarize the regime-specific CCPs, we’ve calculated the participation profile that would result if a woman stayed in each regime for her whole life. We recover the transition matrix well, and the regime-specific CCPs reasonably well, although the estimated initial distribution is less precise. Notice that it took a couple of hundred iterations for the algorithm to converge, compared to a handful for the model with permanent types. This illustrates the slow convergence of EM when the unobserved states are only weakly revealed by the data: a single choice says much less about a woman’s current regime than thirty choices say about her permanent type. How well can we tell which regime a woman is in at a point in time?

accuracy = mean((hmm.Q[2,:,:] .> 0.5) .== (data_hmm.R .== 2))
println("Fraction of person-periods correctly classified: $(round(accuracy,digits=3))")
Fraction of person-periods correctly classified: 0.834

So the posterior classifies the regime correctly a bit more than 80% of the time. This is a useful reminder that the EM algorithm does not need to classify each observation correctly to deliver consistent estimates. It averages over the uncertainty using the posterior weights.

21.5.2 Sparse Implementations

In the example above there were only two unobserved states. In richer applications the unobserved state may have many points of support. For example, if we model unobserved human capital or a persistent wage shock on a grid, the unobserved state may have hundreds of values. The forward step \(\alpha_{t} = (\Pi'\alpha_{t-1})\circ\ell_{t}\) then requires a matrix-vector product that costs \(O(K^{2})\) operations if \(\Pi\) is stored as a dense matrix.

Fortunately, transition matrices in economic models are often very sparse: from any given state, only a handful of states are reachable next period. Storing \(\Pi\) as a sparse matrix reduces the cost of each step to the number of non-zero transitions.

Example 21.6 Consider a hidden state with \(K=1,000\) points where the state can only move to adjacent grid points. Here is a version of the forward pass written with matrix operations, which we can call with either a dense or a sparse transition matrix.

using SparseArrays

function forward_filter!(α,ℓ,Π,μ)
    α[:,1] .= μ .* ℓ[:,1]
    c = sum(α[:,1])
    α[:,1] ./= c
    ll = log(c)
    for t in 2:size(ℓ,2)
        @views mul!(α[:,t],Π',α[:,t-1]) #<- in-place matrix-vector product
        @views α[:,t] .*= ℓ[:,t]
        c = sum(@view α[:,t])
        @views α[:,t] ./= c
        ll += log(c)
    end
    return ll
end

K = 1000
T = 50
Π_dense = zeros(K,K)
for k in 1:K, k′ in max(1,k-1):min(K,k+1)
    Π_dense[k,k′] = 1.
end
Π_dense ./= sum(Π_dense,dims=2)
Π_sparse = sparse(Π_dense)

ℓ = rand(K,T)
α = zeros(K,T)
μ = fill(1/K,K)
forward_filter!(α,ℓ,Π_dense,μ); forward_filter!(α,ℓ,Π_sparse,μ) #<- compile
@time ll_dense = forward_filter!(α,ℓ,Π_dense,μ)
@time ll_sparse = forward_filter!(α,ℓ,Π_sparse,μ)
ll_dense - ll_sparse
  0.000861 seconds (9 allocations: 23.797 KiB)
  0.000131 seconds (9 allocations: 23.797 KiB)
0.0

The answers are identical, but the sparse version is much faster. The gains grow with \(K\): the dense product scales with \(K^{2}\) while the sparse product scales with \(K\). In a full EM implementation this calculation is repeated for every unit in every iteration.

21.6 Exercises

Exercise

Exercise 21.1 Take the labor supply model from Exercise 16.2 and introduce some time-invariant unobserved heterogeneity by allowing the cost parameter \(\kappa\) to depend on an individual’s type \(k\in\{1,2,...,K\}\).

  1. Write code to simulate panel data \(\{d_{n,t},w_{n,t}\}_{n=1,t=1}^{N,T}\) from a model with \(K=2\), where \(d_{n,t}\) indicates whether agent \(n\) worked in period \(t\). You can choose your own parameters for \((\kappa_{1},\kappa_{2})\) and assume \(\pi_{1}=\pi_{2}=0.5\). To avoid the initial conditions problem, assume that everyone in the data begins with \(d_{n,0}=0\).
  2. Estimate the model ignoring heterogeneity. Which parameters are biased, and in what direction?
  3. Now consider the problem of estimating the type probabilities and the type-specific CCPs: \[ \mathbf{P} = \{\pi_{k}, P[d=1|d_{t-1},k,w] \}_{k\in\{1,2\},\ w\in\{w_{L},w_{H}\}, d_{t-1}\in\{0,1\}}. \] Write code to estimate these probabilities using the EM algorithm, assuming that \(\pi\) and \((w_{L},w_{H})\) are known. Show that the M-step for the CCPs is a weighted frequency estimator.
  4. Use your first-stage estimates to estimate \((\kappa_{1},\kappa_{2},\sigma)\) using a CCP estimator from Exercise 20.1.
  5. Demonstrate the performance of your estimation algorithm by running a Monte Carlo simulation. Be careful about label switching!
Exercise

Exercise 21.2 In this exercise you will estimate a simple hidden Markov model with missing data. There is a discrete state variable \(k\in\{1,2,3,..,K\}\) and a binary outcome \(j\in\{0,1\}\) that:

  1. Is determined probabilistically by the state; and
  2. Moves the state up one grid point if \(j=1\), i.e. \(k_{t+1} = \min\{K,k_{t} + j\}\).

Let \(p\) be a \(K\)-dimensional vector where \(p_{k}\) holds the probability that \(j=1\) given that the model is in state \(k\). The state \(k\) is never observed and the outcome \(j\) is only observed half the time. Thus, define \(j^*\) to be: \[ j^* = \left\{ \begin{array}{cl} j & \text{with probability}\ 0.5 \\ -1 & \text{with probability}\ 0.5 \end{array} \right. \] Your task is to estimate the vector \(p\) non-parametrically using the EM algorithm, with the forward-backward algorithm to calculate the E-step. The code below simulates panel data \((j^*_{nt})_{t=1,n=1}^{T,N}\), assuming that \(k_{n1} = 1\) for all \(n\).

using Random, Distributions
K = 5
p_true = 1 ./ (1 .+ exp.(LinRange(-1,1,K))) #<- probability that j=1 in each state
k_next(k,j,K) = min(K,k+j)

function simulate(N,T,p)
    J = zeros(Int,T,N)
    K = length(p)
    for n in axes(J,2)
        k = 1
        for t in axes(J,1)
            j = rand() < p[k]
            J[t,n] = rand() < 0.5 ? -1 : j #<- j is missing half of the time
            k = k_next(k,j,K)
        end
    end
    return J
end

J_data = simulate(1000,10,p_true)

Since the outcomes are periodically unobserved, there are two ways to set up the problem: (1) define a composite hidden state \(s = (k,j)\) and run the forward-backward algorithm over this composite state; or (2) run the forward-backward algorithm over \(k\) only, and sum over potential realizations of \(j\) when it is missing. Since \(j\) is binary here, we’ll use the second approach.

  1. Write a function that (given a guess of \(p\)) takes a sequence \((j^*_{t})_{t=1}^{T}\) and runs the forward-backward algorithm. Your function should fill in: (a) a \(K\times T\) array of forward probabilities, \(\alpha_{t}(k)=\mathbb{P}[k_{t}=k,j^*_{1},\ldots,j^*_{t}]\); (b) a \(K\times T\) array of backward probabilities, \(\beta_{t}(k) = \mathbb{P}[j^*_{t+1},\ldots,j^*_{T} | k_{t}=k]\); and (c) a \(K\times T\) array of posterior probabilities \(q_{t}(k) = \mathbb{P}[k_{t}=k | j^*_{1},\ldots,j^*_{T}]\). Hint: notice that the transition of the hidden state depends on the outcome. It may help to construct three transition matrices, one for each value of \(j^*\in\{-1,0,1\}\).
  2. Write a function that iterates over all observations and fills in a \(K\times T \times N\) array of posterior weights.
  3. Given posterior weights \(q_{ntk}\), the expected log-likelihood for the M-step is: \[ \mathcal{L}(p) = \sum_{n}\sum_{t}\sum_{k}q_{ntk}\left(\mathbf{1}\{j^*_{nt}=1\}\log(p_{k}) + \mathbf{1}\{j^*_{nt}=0\}\log(1-p_{k})\right). \] Show that the maximizer is a frequency estimator: \[ \hat{p}_{k} = \frac{\sum_{n}\sum_{t}\mathbf{1}\{j^*_{nt}=1\} q_{ntk}}{\sum_{n}\sum_{t}\mathbf{1}\{j^*_{nt}\neq -1\} q_{ntk}} \] and write a function to calculate it.
  4. Run the EM routine by iterating on the E-step and M-step until \(\hat{p}\) converges. Does it look like you can recover the true parameters? Run a Monte Carlo simulation to check.
  5. (Extra) Adjust the model so that the state moves up with probability \(\lambda_{k}\) when \(j=1\) (and stays put otherwise). Extend your algorithm to estimate \(\lambda\) as part of the M-step.
  6. (Extra) Rewrite your forward-backward routine using sparse transition matrices and compare its speed to a dense implementation for \(K=200\).
Exercise

Exercise 21.3 In Example 21.5 we assumed that women do not anticipate changes in the demands on their time at home. Suppose instead that women know the transition matrix \(\Pi\) and take it into account when making decisions.

  1. Write down the conditional value functions for this model. What is the state space?
  2. Show that the model still exhibits two-period finite dependence, and derive an expression for \(v_{1t}(h,y,k,\ell)-v_{0t}(h,y,k,\ell)\) in terms of flow payoffs, the transitions of \((y,k)\), \(\Pi\), and the CCPs.
  3. Simulate data from this version of the model. Use the EM algorithm from Example 21.5 to estimate the regime-specific CCPs and \(\Pi\), and then estimate \((\alpha_{01},\alpha_{02},\alpha_{1},\sigma)\) in a second stage.
  4. Does the first stage of the EM algorithm need to change? Why or why not?
Exercise

Exercise 21.4 Recall our discussion of identification of the search model with unobserved heterogeneity, where the distribution of unemployment durations is a finite mixture of geometric distributions: \[ \mathbb{P}(t_U=t) = \sum_{k}\tilde{\pi}_{k}h_{k}(1-h_{k})^{t}. \]

  1. Using the CPS data from Example 8.1, select a sample of unemployed workers without a college degree and construct their (elapsed) unemployment durations using the variable DURUNEMP.
  2. Write an EM algorithm to estimate \((\tilde{\pi}_{k},h_{k})_{k=1}^{K}\). Derive the M-step for \(h_{k}\) in closed form.
  3. Estimate the model for \(K=1,2,3,4\) and use the BIC to choose \(K\).
  4. Compare your preferred estimates to those for workers with a college degree. What do the differences suggest about the sources of differences in unemployment durations across education groups?