A Simple Labor Supply Model

Model Setup and Solution

Consider a dynamic labor supply model (with no uncertainty) where each agent \(n\) chooses a sequence of consumption and hours, \(\{c_{t},h_{t}\}_{t=1}^{\infty}\), to solve: \[ \max \sum_{t=0}^\infty \beta^{t} \left(\frac{c_{t}^{1-\sigma}}{1-\sigma} - \frac{\alpha_{n}^{-1}}{1 + 1/\psi}h_{t}^{1+1/\psi}\right)\] subject to the intertemporal budget constraint: \[ \sum_{t}q_{t}c_{t} \leq A_{n,0} + \sum_{t}q_{t}W_{n,t}h_{t},\qquad q_{t} = (1+r)^{-t}.\] Let \(H_{n,t}\) and \(C_{n,t}\) be the realizations of labor supply for agent \(n\) at time \(t\). Labor supply in this model obeys: \[H_{n,t}^{1/\psi} = (\alpha_{n}W_{n,t})C^{-\sigma}_{n,t}.\] To simplify below, assume that \(\beta=(1+r)^{-1}\), so that the optimal solution features perfectly smoothed consumption, \(C^*_{n}\). Making appropriate substitutions gives \(C^*_{n}\) as the solution to: \[ \left(\sum_{t}q_{t}\right)C^*_{n} = \sum_{t}\left(q_{t}W_{n,t}^{1+\psi}\right)\alpha_{n}^{\psi}(C_{n}^*)^{-\psi\sigma} + A_{n,0}.\]

Code to solve the model

There is only one object to solve here which is consumption given a sequence of net wages. If one were to assume also constant wages (or net ages in the case of Assignment 2, the function below solves optimal consumption.

using Optim
function solve_consumption(r,α,W,A,σ,ψ)
    Q = 1/ (1 - 1/(1+r))
    f(c) = (Q * c - Q * W^(1 + ψ) * α^ψ * c^(-σ*ψ) - A)^2
    r = Optim.optimize(f,0.,A+W)
    return r.minimizer
end
solve_consumption (generic function with 1 method)

Code to simulate a cross-section

Here we’ll assume that wages, tastes for work, and assets co-vary systematically. For simplicity we’ll use a multivariate log-normal distribution.

Below is code to simulate a cross-section of 1,000 observations.

using Distributions
function simulate_data(σ,ψ,r,N)
    ch = [0.3 0. 0.; 0.5 0.5 0.; 0.4 0.8 1.8]
    Σ = ch * ch'
    X = rand(MvNormal(Σ),N)
    α = exp.(X[1,:])
    W = exp.(X[2,:])
    A = exp.(X[3,:])
    C = [solve_consumption(r,α[i],W[i],A[i],σ,ψ) for i in eachindex(A)]
    @views H = exp.( X[1,:] .+ ψ .* X[2,:] .- ψ * σ .* log.(C) )
    return (;α,W,A,C,H)
end

# assume risk-aversion of 2 and frisch of 0.5
σ = 2.
ψ = 0.5
r = 0.05

dat = simulate_data(σ,ψ,r,1_000)
(α = [1.2009314280657941, 1.3840635001194714, 0.7932036143727019, 0.6496050659926231, 1.1763954044345672, 0.757638882397994, 0.9044004976586716, 0.6364941508058962, 0.9848492547280477, 0.7665715755742666  …  0.7913335134653614, 0.9185266909357928, 1.1221059739006543, 0.784576507966877, 0.9085786433280679, 1.481766688920262, 1.270695828772059, 1.186733064198457, 1.3704311810289025, 1.2751060776129926], W = [2.0926843882464046, 1.2294345771934665, 0.7839372444902444, 0.4932246921798898, 1.4476323517690834, 0.5990320722671819, 1.3574418855472676, 0.4501202284057226, 0.648350821898256, 0.4438307015980542  …  0.9887124774895775, 0.9792717737257789, 1.374143833129326, 0.5262309716779926, 0.7540809731942225, 1.9551883554551543, 1.5747671789102187, 0.5076303709683835, 0.6007997162881663, 0.8540316558671761], A = [0.42717686715183817, 0.9004583983324581, 2.0549539924104545, 0.18263287570010597, 0.1811592221765357, 2.8252890589639694, 3.002328213654121, 0.23715153542691714, 0.01948962388559997, 0.5256556088827199  …  0.22146690771594454, 0.8597978798731503, 1.9835320101878544, 2.6752159614019377, 0.7087343049240805, 3.4831868622571656, 8.648948202054312, 0.04856734023740389, 0.02574314603678183, 0.008057502896994183], C = [1.8316077596136042, 1.2880126944326615, 0.8366920060842069, 0.5327457924842408, 1.3787798723691504, 0.7060828047981375, 1.2999627020592937, 0.49652465678494007, 0.6678404265684732, 0.5214742973789857  …  0.9404618377737518, 0.9844076729348431, 1.3543493773467359, 0.6486622886044833, 0.8071039843110066, 1.9090757912913026, 1.7125950958734384, 0.5561976954282956, 0.6265428447358602, 0.8620891340767702], H = [0.9485005843996103, 1.1914850131607146, 0.8393821198757727, 0.8563509422325061, 1.026567532201603, 0.8304847406162642, 0.8105699954999359, 0.8600384392333683, 1.1874137910752485, 0.9793286828309261  …  0.8366684431108959, 0.9233543700260457, 0.971223671050765, 0.8774145921152507, 0.9775568526887064, 1.0853028860714815, 0.9310973392526273, 1.5201892474294891, 1.695396632883908, 1.3668823286993619])

Simulating bias in OLS

We know that OLS will be biased here, but let’s run a simulation anyway to hammer the point home.

function simulate_bias(B,N,σ,ψ,r)
    ψ_est = zeros(B)
    σ_est = zeros(B)
    for b in 1:B
        dat = simulate_data(σ,ψ,r,N)
        X = [ones(N) log.(dat.W) log.(dat.C)]
        Y = log.(dat.H)
        b_est = inv(X' * X) * X' * Y
        ψ_est[b] = b_est[2]
        σ_est[b] = -b_est[3] / b_est[2]
    end
    return ψ_est,σ_est
end

ψ_est,σ_est = simulate_bias(500,1_000,σ,ψ,r)
([0.7223802219229686, 0.7139442086800205, 0.7396470356525909, 0.7246208926992348, 0.7726181461197881, 0.7276328710547488, 0.7317515458346843, 0.7404793482805629, 0.7798694647971162, 0.6973219176194752  …  0.6665532380154393, 0.7734667272980555, 0.7710703160005671, 0.739778694895188, 0.7414361916136651, 0.7390937737267499, 0.7640303196929837, 0.7342307577101521, 0.7427383479952441, 0.776025027636072], [1.2757790264104605, 1.2804478804796597, 1.259854310253586, 1.2461437689542316, 1.2548669801835217, 1.2453144000610175, 1.2734763370834474, 1.252898187072911, 1.2404787546913614, 1.2801763411638727  …  1.2906778550433584, 1.2431398509674618, 1.2427656037727224, 1.2365281668327026, 1.244401467239429, 1.262850907861914, 1.251095593857676, 1.2717938681082661, 1.2662647469210373, 1.2497172792363864])

Given the covariance structure, we can see a clear bias in the estimates.

using Plots
histogram(ψ_est)