22 Discretizing Heterogeneity with Clustering

The EM algorithm treats each unit’s type as a latent variable and integrates over it using posterior probabilities. This is a powerful approach, but it has some drawbacks. It requires specifying the full likelihood, it can be slow to converge, it must be run from many starting values, and the number of parameters grows quickly with the number of types. An alternative strategy, which has become popular in recent years, is to classify units into groups using the data, and then estimate the model treating the group memberships as known. This chapter introduces two versions of this idea: the grouped fixed effects estimator of Bonhomme and Manresa (2015) and the two-step approach of Bonhomme, Lamadon, and Manresa (2022).

Both approaches rely on the panel dimension of the data to classify units. The basic logic is that with enough observations per unit, the data reveal which group each unit belongs to with high probability. The key questions are: how large does \(T\) need to be, and what should we use to classify?

22.1 Grouped Fixed Effects

Bonhomme and Manresa (2015) consider linear panel data models of the form:

\[ y_{nt} = \mathbf{x}_{nt}'\theta + \alpha_{g(n),t} + v_{nt} \]

where \(g(n)\in\{1,\ldots,G\}\) is the (unknown) group membership of unit \(n\) and \(\alpha_{gt}\) is a group-specific time effect. The grouped fixed effects (GFE) estimator solves:

\[ (\hat{\theta},\hat{\alpha},\hat{g}) = \arg\min_{\theta,\alpha,g}\sum_{n=1}^{N}\sum_{t=1}^{T}\left(y_{nt} - \mathbf{x}_{nt}'\theta - \alpha_{g(n),t}\right)^{2}. \]

Notice that the group memberships \(g\) are treated as parameters to be estimated, just like fixed effects. The difference from a standard fixed effects estimator is that the number of distinct values of the unobserved heterogeneity is restricted to \(G\), so the number of parameters does not grow with \(N\).

This problem can be solved by iterating between two steps, in the spirit of the k-means algorithm:

Estimation: given group memberships \(g^{n}\), estimate \((\theta^{n+1},\alpha^{n+1})\) by least squares with group-time dummies.
Assignment: given \((\theta^{n+1},\alpha^{n+1})\), assign each unit to the group that minimizes its sum of squared residuals: \[ g^{n+1}(n) = \arg\min_{g}\sum_{t}\left(y_{nt} - \mathbf{x}_{nt}'\theta^{n+1} - \alpha^{n+1}_{gt}\right)^{2}. \]

Each step weakly reduces the objective, so the algorithm converges, but (like EM) only to a local minimum. Bonhomme and Manresa (2015) recommend running it from many starting values.

How does the GFE estimator behave? The key result in Bonhomme and Manresa (2015) is that, as long as groups are well-separated, the probability of misclassifying any given unit vanishes exponentially fast as \(T\) grows. As a consequence, when \(N\) and \(T\) grow together (with \(N/T^{\nu}\rightarrow0\) for some \(\nu>0\)), \(\hat{\theta}\) is asymptotically equivalent to the infeasible estimator that knows the true group memberships. In other words, the classification error does not affect inference asymptotically. In finite samples with short panels, however, misclassification can be substantial.

Example: Classifying with Means versus Least Squares

Example 22.1 Consider a simple model with heterogeneous coefficients:

\[ y_{nt} = \beta_{0,k(n)} + \beta_{1,k(n)}x_{nt} + \epsilon_{nt} \]

where \(x_{nt}\) is an AR(1) process:

\[ x_{nt} = \rho x_{nt-1} + \xi_{nt}\qquad \xi_{nt}\sim\mathcal{N}(0,\sigma^2_\xi). \]

Here’s code to simulate data from the model with two types.

using Random, Distributions, Clustering, LinearAlgebra, Statistics

function simulate_panel(N,T,pars)
    (;ρ,β,σξ,σϵ) = pars
    K = size(β,2) #<- β[:,k] holds the intercept and slope for type k
    X = zeros(T,N)
    Y = zeros(T,N)
    type = rand(Categorical(fill(1/K,K)),N)
    σ_stat = sqrt(σξ^2/(1-ρ^2)) #<- stationary standard deviation of x
    for n in axes(X,2)
        k = type[n]
        x = rand(Normal(0,σ_stat))
        for t in axes(X,1)
            X[t,n] = x
            Y[t,n] = β[1,k] + β[2,k]*x + rand(Normal(0,σϵ))
            x = ρ*x + rand(Normal(0,σξ))
        end
    end
    return (;type,X,Y)
end

pars = (;β = [0. 0.5 ; 0.3 0.1], ρ = 0.95, σξ = 0.5, σϵ = 0.2)
Random.seed!(1234)
data = simulate_panel(1000,10,pars);

A natural first idea is to cluster units based on their outcomes \(y_{nt}\), using k-means. Since types are only identified up to a relabeling, we measure the classification accuracy as the share of units correctly classified under the best matching of estimated groups to true types.

# fraction correctly classified, allowing for relabeling of the two groups
accuracy(g,type) = max(mean(g .== type), mean(g .== 3 .- type))

res = kmeans(data.Y,2)
accuracy(res.assignments,data.type)

0.797

We misclassify many units. The problem is that the outcomes \(y_{nt}\) depend on the observed covariate \(x_{nt}\), which is very persistent. A unit with high values of \(x\) throughout the sample will look like a “high” type, regardless of its true type. In the language of Bonhomme, Lamadon, and Manresa (2022), the moments we’re using to classify are not informative enough about the unobserved type. The histograms of average outcomes by type make this clear:

using Plots
ȳ = mean(data.Y,dims=1)[:]
histogram(ȳ[data.type.==1],alpha=0.5,label="Type 1",normalize=:pdf)
histogram!(ȳ[data.type.==2],alpha=0.5,label="Type 2",normalize=:pdf,xlabel="Mean of y")

A better approach is to use the structure of the model and classify based on the least squares criterion itself:

\[ \hat{\mathcal{K}},\hat{\beta} = \arg\min\sum_{n}\sum_{t}(y_{nt} - \beta_{0,k(n)} - \beta_{1,k(n)}x_{nt})^2 \]

which we can solve by iterating in the same way as k-means. Here’s code to do that, using the k-means classification as a starting point.

# assign each unit to the type that minimizes its sum of squared residuals
function assign!(type,data,β)
    K = size(β,2)
    ssr = zeros(K)
    for n in axes(data.X,2)
        x = view(data.X,:,n)
        y = view(data.Y,:,n)
        for k in 1:K
            ssr[k] = sum((y .- β[1,k] .- β[2,k] .* x).^2)
        end
        type[n] = argmin(ssr)
    end
end

# estimate the coefficients for each type by OLS given the assignments
function estimate!(type,data,β)
    for k in axes(β,2)
        Ik = type.==k
        y = data.Y[:,Ik][:]
        X = [ones(length(y)) data.X[:,Ik][:]]
        β[:,k] = X \ y
    end
end

function classify_least_squares!(type,β,data; max_iter = 100)
    for iter in 1:max_iter
        type_old = copy(type)
        estimate!(type,data,β)
        assign!(type,data,β)
        if type==type_old #<- stop when the assignments no longer change
            return iter
        end
    end
    println("Warning: maximum iterations reached")
    return max_iter
end

β̂ = zeros(2,2)
type = copy(res.assignments)
iterations = classify_least_squares!(type,β̂,data)
println("Converged in $iterations iterations")
println("Classification accuracy: $(accuracy(type,data.type))")
β̂

Converged in 3 iterations
Classification accuracy: 0.99

2×2 Matrix{Float64}:
 0.499351   -0.00332652
 0.0998685   0.297928

We do much better now, and the estimated coefficients are close to the truth (up to the labeling of types). The lesson is that what we use to classify matters: we want statistics that vary with the unobserved type, and only with the unobserved type.

22.2 A Two-Step Approach

The GFE estimator requires us to solve a joint classification and estimation problem, which can be costly in structural models. Bonhomme, Lamadon, and Manresa (2022) propose a simpler alternative:

Classification: for each unit, compute a vector of moments \(\mathbf{h}_{n} = \frac{1}{T}\sum_{t}h(\mathbf{w}_{nt})\) and classify units into \(K\) groups using k-means: \[ (\hat{k},\hat{\mathbf{h}}) = \arg\min_{k(1),\ldots,k(N),\tilde{\mathbf{h}}_{1},\ldots,\tilde{\mathbf{h}}_{K}}\sum_{n=1}^{N}\|\mathbf{h}_{n} - \tilde{\mathbf{h}}_{k(n)}\|^{2}. \]
Estimation: estimate the structural model treating \(\hat{k}(n)\) as an observed type, allowing parameters to vary by group.

The approach rests on an injectivity condition: the population moments \(\text{plim}_{T\rightarrow\infty}\mathbf{h}_{n} = \varphi(\alpha_{n})\) should be a one-to-one function of the unit’s unobserved heterogeneity \(\alpha_{n}\). This rules out moments, like the mean of \(y\) in Example 22.1, that also vary systematically with observables. A useful strategy is to classify on residuals from a first-stage model that conditions on observables.

Bonhomme, Lamadon, and Manresa (2022) show that this approach has an attractive interpretation even when the true heterogeneity is continuous. In that case, the groups provide a discrete approximation to the heterogeneity, and the approximation error shrinks as \(K\) grows. The resulting estimators have a bias of order \(1/T\) (plus an approximation error that is negligible if \(K\) grows fast enough), which is similar to fixed effects estimators in nonlinear panel models. They also propose a bias-correction method.

This last point deserves emphasis in dynamic models. Many of the observed states in a dynamic model, like accumulated work experience, are themselves the outcome of past choices, and so they depend on the unobserved type. Residuals from a first stage that conditions on these endogenous states will partly purge the very variation in behavior that reveals the type. We will see this concretely in the example below.

The two-step approach has been particularly influential in labor economics. Bonhomme, Lamadon, and Manresa (2022) originally developed it to estimate models of the labor market in which both workers and firms are heterogeneous. In that setting, firms are classified based on the distribution of wages they pay, which can be estimated precisely because each firm employs many workers, and the firm classes are then treated as known when estimating worker heterogeneity and the returns to firm type.

Example: Two-Step Estimation of the Labor Force Participation Model

Example 22.2 Let’s return to the female labor force participation model and suppose that women belong to one of two types that differ in their taste for work, \(\alpha_{0}\), and in their earnings potential, through the intercept of the wage equation, \(\gamma_{0}\):

\[ \log(w) = \gamma_{0\ell} + \gamma_{1}h - \gamma_{2}h^2,\qquad u_{0}(h,y,k;\ell) = \log(y) + \alpha_{0\ell} + \alpha_{1}k, \]

where \(\ell\in\{1,2\}\) indexes the type and, as before, \(k\) indicates the presence of a young child.

Type 1 women have lower wages and a higher value of time at home, while type 2 women have higher wages and a lower value of time at home. Both differences make type 2 women more likely to work. We’ll solve the model for each type and simulate a panel of 4,000 women over their full life-cycle.

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

β = 0.95
γ_types = [(;γ0 = log(1.5) - 0.15, γ1 = 0.08, γ2 = 0.0015),
           (;γ0 = log(1.5) + 0.15, γ1 = 0.08, γ2 = 0.0015)]
θ_types = [(;α0 = 0.45, α1 = 0.3, σ = 0.25),
           (;α0 = 0.15, α1 = 0.3, σ = 0.25)]
env_types = [flfp_environment(γ = γ_types[k]) for k in 1:2]
sol_types = [solve_flfp(θ_types[k],env_types[k],β) for k in 1:2]

# simulate N women, each drawn to be of type k with probability πk[k]
function simulate_flfp_types(sol_types,env_types,πk,N; seed = nothing)
    if !isnothing(seed)
        Random.seed!(seed)
    end
    type = rand(Categorical(πk),N)
    T = env_types[1].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)
    for k in eachindex(πk)
        Ik = findall(type.==k)
        dk = simulate_flfp(sol_types[k].logp1,env_types[k],length(Ik))
        H[:,Ik] .= dk.H; Y[:,Ik] .= dk.Y; K[:,Ik] .= dk.K; D[:,Ik] .= dk.D
        logW[:,Ik] .= dk.logW
    end
    return (;H,Y,K,D,logW,type)
end

# keep the first T periods of the panel, or a subset I of women
window(data,T) = (;H = data.H[1:T,:], Y = data.Y[1:T,:], K = data.K[1:T,:], D = data.D[1:T,:], logW = data.logW[1:T,:], type = data.type)
subset(data,I) = (;H = data.H[:,I], Y = data.Y[:,I], K = data.K[:,I], D = data.D[:,I], logW = data.logW[:,I])

data_flfp = simulate_flfp_types(sol_types,env_types,[0.5,0.5],4000; seed = 1)
println("Participation rate by type: $(round.([mean(data_flfp.D[:,data_flfp.type.==k]) for k in 1:2],digits=3))")

Participation rate by type: [0.363, 0.801]

Before we classify, notice that heterogeneity contaminates even the simplest part of the model. If we estimate the wage equation by OLS on workers, ignoring types, we get:

estimate_wages(data_flfp)

(γ0 = 0.38768502030470703, γ1 = 0.08680077853043323, γ2 = 0.0014476043237247429)

The return to experience is overstated (the truth is \(\gamma_{1}=0.08\)). This is a selection problem: high-wage women work more and so accumulate more experience, so experience is positively correlated with the unobserved wage intercept. The measurement error assumption that justified OLS in the original model is no longer enough once there is persistent heterogeneity in wages.

For classification, we’ll use two moments for each woman. The first is her average participation residual, \(\frac{1}{T}\sum_{t}(D_{nt}-\hat{p}_{nt})\), where \(\hat{p}_{nt}\) comes from a flexible logit first stage that ignores types. The second is her average wage residual from the pooled wage regression, across the periods in which she works (women who never work are assigned the average residual of zero). Following the discussion above, we will compare a first stage that conditions only on the exogenous states (husband’s income, children, and age) to one that also conditions on experience.

# flexible logit for the probability of working, ignoring types
function first_stage(data; endogenous = false)
    (;H,Y,K,D) = data
    T,N = size(D)
    features(h,iy,k,t) = endogenous ? [1., h, h^2/10, iy, k, t, t^2/10, h*t/10, k*t] : [1., iy, k, t, t^2/10, k*t]
    Z = reduce(hcat,[features(H[t,n],Y[t,n],K[t,n],t) for t in 1:T, n in 1:N][:])'
    d = D[:]
    ll(b) = (v = Z*b; sum(d .* v .- log1p.(exp.(v))))
    b = optimize(b->-ll(b),zeros(size(Z,2)),Newton(),autodiff=:forward).minimizer
    return reshape(1 ./ (1 .+ exp.(-Z*b)),T,N)
end

# participation and wage residual moments for each woman
function moments(data; endogenous = false)
    (;H,D,logW) = data
    N = size(D,2)
    p̂ = first_stage(data; endogenous)
    hD = vec(mean(D .- p̂,dims=1))
    γ̂ = estimate_wages(data)
    r = logW .- (γ̂.γ0 .+ γ̂.γ1 .* H .- γ̂.γ2 .* H.^2)
    hW = [any(D[:,n].==1) ? mean(r[D[:,n].==1,n]) : 0. for n in 1:N]
    return hD,hW
end

# k-means on standardized moments (best of 10 starts), groups ordered by participation residual
function classify(h; K = 2)
    hs = (h .- mean(h,dims=2)) ./ std(h,dims=2)
    best = kmeans(hs,K)
    for r in 2:10
        res = kmeans(hs,K)
        if res.totalcost < best.totalcost
            best = res
        end
    end
    g = best.assignments
    order = sortperm([mean(h[1,g.==k]) for k in 1:K])
    return invperm(order)[g]
end

Random.seed!(2)
println("Classification accuracy by moments used:")
for T in [5,10,20,30]
    d = window(data_flfp,T)
    hD,hW = moments(d)
    hD_endog,_ = moments(d; endogenous = true)
    acc = [accuracy(classify(reshape(hD,1,:)),d.type),
           accuracy(classify(reshape(hW,1,:)),d.type),
           accuracy(classify([hD hW]'),d.type),
           accuracy(classify([hD_endog hW]'),d.type)]
    println("T = $T: participation = $(round(acc[1],digits=3)), wages = $(round(acc[2],digits=3)), both = $(round(acc[3],digits=3)), both (conditioning on experience) = $(round(acc[4],digits=3))")
end

Classification accuracy by moments used:
T = 5: participation = 0.847, wages = 0.682, both = 0.88, both (conditioning on experience) = 0.87
T = 10: participation = 0.918, wages = 0.784, both = 0.956, both (conditioning on experience) = 0.943
T = 20: participation = 0.97, wages = 0.885, both = 0.991, both (conditioning on experience) = 0.956
T = 30: participation = 0.984, wages = 0.93, both = 0.997, both (conditioning on experience) = 0.972

A few lessons emerge. First, both moments are informative, and combining them improves classification at every panel length. Wages are less informative than participation in short panels, since we only observe wages in periods of work and they are measured with considerable noise, but they become very informative as \(T\) grows. Second, classification improves quickly with \(T\) when we condition only on exogenous states. Third, conditioning on experience in the first stage hurts, and increasingly so as \(T\) grows. Experience is a summary of a woman’s past participation, which is driven by her type. A woman with a strong taste for work accumulates experience and so her first-stage predicted probability of working is high: the first stage “explains” her persistent participation with her experience, leaving a residual that says little about her type. This is a violation of the injectivity condition.

Now let’s implement the second step. Given the classification, we estimate the wage equation with group-specific intercepts, and then estimate the preference parameters by NFXP with a group-specific \(\alpha_{0}\), solving the model separately for each group.

# wage equation with group-specific intercepts and common returns to experience
function estimate_wages_groups(data,g,G)
    (;H,D,logW) = data
    idx = findall(D.==1)
    n = [i[2] for i in idx] #<- the woman in each observation
    X = [reduce(hcat,[Float64.(g[n] .== k) for k in 1:G]) H[idx] -H[idx].^2]
    b = X \ logW[idx]
    return [(;γ0 = b[k], γ1 = b[G+1], γ2 = b[G+2]) for k in 1:G]
end

# log-likelihood with x = (α0 for each group..., α1, log σ)
function log_likelihood_groups(x,data,g,γ_groups,β)
    G = length(γ_groups)
    ll = zero(eltype(x))
    for k in 1:G
        θ = (;α0 = x[k], α1 = x[G+1], σ = exp(x[G+2]))
        env = flfp_environment(γ = γ_groups[k])
        (;logp1,logp0) = solve_flfp(θ,env,β)
        ll += choice_loglikelihood(logp1,logp0,subset(data,g.==k))
    end
    return ll
end

function estimate_given_groups(data,g,β; G = 2)
    γ_groups = estimate_wages_groups(data,g,G)
    x0 = [zeros(G); 0.; log(0.5)]
    x̂ = optimize(x->-log_likelihood_groups(x,data,g,γ_groups,β),x0,Newton(),autodiff=:forward).minimizer
    return (;γ0 = [γ.γ0 for γ in γ_groups], γ1 = γ_groups[1].γ1, α0 = x̂[1:G], α1 = x̂[G+1], σ = exp(x̂[G+2]))
end

function two_step(data,β; G = 2)
    hD,hW = moments(data)
    g = classify([hD hW]'; K = G)
    return estimate_given_groups(data,g,β; G), g
end

two_step (generic function with 1 method)

We’ll compare the two-step estimates to an infeasible estimator that knows each woman’s true type. The difference between the two isolates the effect of classification error, since both estimators share the same sampling variation.

show_est(e) = "γ0 = $(round.(e.γ0 .- log(1.5),digits=3)), γ1 = $(round(e.γ1,digits=4)), α0 = $(round.(e.α0,digits=3)), α1 = $(round(e.α1,digits=3)), σ = $(round(e.σ,digits=3))"
Random.seed!(3)
for T in [5,10,20,30]
    d = window(data_flfp,T)
    est,g = two_step(d,β)
    est_inf = estimate_given_groups(d,d.type,β)
    println("T = $T (accuracy = $(round(accuracy(g,d.type),digits=3)))")
    println("   two-step:   ",show_est(est))
    println("   infeasible: ",show_est(est_inf))
end

T = 5 (accuracy = 0.88)
   two-step:   γ0 = [-0.182, 0.17], γ1 = 0.0665, α0 = [0.441, 0.09], α1 = 0.33, σ = 0.269
   infeasible: γ0 = [-0.141, 0.151], γ1 = 0.0803, α0 = [0.422, 0.07], α1 = 0.304, σ = 0.244
T = 10 (accuracy = 0.956)
   two-step:   γ0 = [-0.158, 0.159], γ1 = 0.0759, α0 = [0.453, 0.166], α1 = 0.328, σ = 0.267
   infeasible: γ0 = [-0.145, 0.155], γ1 = 0.0771, α0 = [0.455, 0.159], α1 = 0.324, σ = 0.266
T = 20 (accuracy = 0.991)
   two-step:   γ0 = [-0.149, 0.154], γ1 = 0.0796, α0 = [0.452, 0.156], α1 = 0.304, σ = 0.251
   infeasible: γ0 = [-0.147, 0.154], γ1 = 0.0794, α0 = [0.454, 0.143], α1 = 0.309, σ = 0.255
T = 30 (accuracy = 0.997)
   two-step:   γ0 = [-0.146, 0.152], γ1 = 0.0795, α0 = [0.456, 0.149], α1 = 0.304, σ = 0.253
   infeasible: γ0 = [-0.145, 0.153], γ1 = 0.0795, α0 = [0.456, 0.143], α1 = 0.307, σ = 0.255

The wage intercepts are reported relative to \(\log(1.5)\), so the true values are \(\gamma_{0}=(-0.15,0.15)\), \(\gamma_{1}=0.08\), \(\alpha_{0}=(0.45,0.15)\), \(\alpha_{1}=0.3\), and \(\sigma=0.25\). With \(T=20\) or more, the two-step estimator is essentially indistinguishable from the infeasible one. With \(T=5\), about one in eight women is misclassified, and the effect on the wage equation is instructive. The estimated wage intercepts are further apart than the true ones, and the return to experience is understated. The reason is that classification is based on realized outcomes. A woman who happened to draw high wage residuals (or who happened to work a lot) in a short panel is likely to be assigned to the high type. Classification therefore attributes too much of the variation in wages to persistent heterogeneity and too little to experience. This is the discrete analogue of the well-known bias of fixed effects estimators in nonlinear and dynamic panel models with short \(T\). Notice that it pushes the return to experience in the opposite direction to the selection bias from ignoring heterogeneity altogether.

22.3 EM or Clustering?

How should we choose between these approaches? Here are some considerations.

Panel length. The EM algorithm delivers consistent estimates for fixed \(T\), provided the model is correctly specified and the mixture is identified (see the next chapter). Clustering methods rely on \(T\) growing to control misclassification. With short panels, EM is generally preferable.
Computation. Clustering is cheap and only needs to be done once. EM requires solving (or approximating) the model for each type in every iteration, and many iterations may be needed. Even when EM is the final estimator, clustering often provides good starting values.
Robustness. Classification in the two-step approach does not depend on the structural model, which makes the classification step robust to misspecification of the model. On the other hand, the choice of moments matters, and poorly chosen moments lead to poor classification.
Nature of heterogeneity. If heterogeneity is continuous, neither approach is exactly correct. Bonhomme, Lamadon, and Manresa (2022) provide a framework for interpreting the two-step estimator as an approximation in this case.

Discussion

Suppose you are working with a panel of welfare recipients observed monthly for five years, and you want to account for unobserved differences in the cost of working.

What moments would you use to classify individuals? How would you check that they satisfy the injectivity condition?
Which of the two approaches would you choose? Would it matter if you had annual rather than monthly data?
In Example 22.2, wages were informative about type because types differed in their wage intercept. In the welfare setting, many individuals rarely work. What does this imply for the quality of classification, and what other outcomes might you use?
In the matched employer-employee setting of Bonhomme, Lamadon, and Manresa (2022), firms are classified using the wages of their workers. Why is classification likely to be much more accurate for firms than for workers? What plays the role of \(T\)?

22.4 Exercises

Exercise

Exercise 22.1 Return to the labor supply model with two types from Exercise 21.1.

Propose a set of moments to classify individuals into types. Explain why your moments should be informative about the type.
Implement the two-step estimator of Bonhomme, Lamadon, and Manresa (2022) using these moments and the CCP estimator from Exercise 20.1.
Run a Monte Carlo simulation comparing the two-step estimator to the EM estimator from Exercise 21.1, for \(T=5\), \(T=10\), and \(T=30\). Comment on the bias and variance of each estimator.
Use the clustering results to provide starting values for the EM algorithm. How much does this reduce the number of iterations needed for convergence?

Exercise

Exercise 22.2 Using the PSID data from Example 10.1, consider the income process from Example 12.2, but now suppose that individuals belong to one of \(K\) groups that differ in the deterministic component of their life-cycle income profile.

Construct residuals from a regression of log income on a quadratic in age. For each individual, compute the mean residual and the slope of residuals with respect to age.
Classify individuals into \(K=3\) groups using k-means on these moments. Plot the average life-cycle income profile for each group.
Re-estimate the parameters of the income process by minimum distance, using residuals from group-specific income profiles. How do your estimates of the variance of persistent shocks change? Why?
Explain how you would conduct inference that accounts for the classification step.