flowchart TD
A(("Period t<br/>Offer w̄ⱼ"))
A -- "Accept → earn w̄ⱼ" --> B(["Period t+1<br/>Keep w̄ⱼ or quit?"])
A -- "Reject → earn c" --> C(["Period t+1<br/>New offer w̄ₖ"])
B -- "Keep → w̄ⱼ" --> D([t+2])
B -- "Quit → c" --> E([t+2])
C -- "Accept → w̄ₖ" --> F([t+2])
C -- "Reject → c" --> G([t+2])
D --- D1((·)) & D2((·))
E --- E1((·)) & E2((·))
F --- F1((·)) & F2((·))
G --- G1((·)) & G2((·))
style A fill:#154733,color:#fff,stroke:#154733
style B fill:#4a8c6f,color:#fff,stroke:#4a8c6f
style C fill:#4a8c6f,color:#fff,stroke:#4a8c6f
style D fill:#8bb8a0,color:#000,stroke:#8bb8a0
style E fill:#8bb8a0,color:#000,stroke:#8bb8a0
style F fill:#8bb8a0,color:#000,stroke:#8bb8a0
style G fill:#8bb8a0,color:#000,stroke:#8bb8a0
style D1 fill:#c8ddd0,stroke:#8bb8a0
style D2 fill:#c8ddd0,stroke:#8bb8a0
style E1 fill:#c8ddd0,stroke:#8bb8a0
style E2 fill:#c8ddd0,stroke:#8bb8a0
style F1 fill:#c8ddd0,stroke:#8bb8a0
style F2 fill:#c8ddd0,stroke:#8bb8a0
style G1 fill:#c8ddd0,stroke:#8bb8a0
style G2 fill:#c8ddd0,stroke:#8bb8a0
The McCall Search Model
1 Introduction
These notes develop the McCall search model, one of the foundational models in labor economics. The model formalizes a simple but powerful idea: an unemployed worker who receives random wage offers must decide each period whether to accept the current offer or reject it and wait for something better. This trade-off between taking a sure thing now and holding out for a potentially better future outcome is at the heart of search theory.
We will build the model in stages:
- Finite horizon — solve the worker’s problem when there are a fixed number of periods remaining, using backward induction.
- Solving with functions — package our code into reusable Julia functions with named tuple returns.
- Infinite horizon — take the number of periods to infinity and solve the resulting fixed point problem with value function iteration.
- Organizing with structs — group model parameters into a Julia
structfor cleaner code. - Firing risk — extend the model so that employed workers can lose their jobs.
- Comparative statics — examine how changes in unemployment benefits, firing probability, and wage distributions affect outcomes.
By the end of these notes you should be comfortable setting up and solving dynamic programming problems in Julia, and you should have a solid understanding of the economics of job search.
2 Model Setup
2.1 The Environment
A worker lives for \(T\) periods (possibly infinite) and can be either unemployed or employed. The worker discounts the future at rate \(\beta \in (0, 1)\). Each period, an unemployed worker receives a random wage offer \(\bar{w}_j\) with probability \(p_j\), where \(j \in \{1, \ldots, S\}\) indexes the possible wage levels. The worker then faces a binary choice:
- Accept the offer: receive wage \(\bar{w}_j\) and become employed.
- Reject the offer: receive unemployment benefits \(c\) and remain unemployed, drawing a new offer next period.
This binary choice creates a branching decision tree. Each period the worker faces the same accept-or-reject decision, so the number of possible life paths grows exponentially with the horizon:
At the root node (period \(t\)) the worker receives offer \(\bar{w}_j\) and makes one choice, producing two branches. Each of those branches leads to a new decision in period \(t+1\)—four possible states. By period \(t+2\) there are eight, and by period \(t+3\) there are sixteen. The tree grows exponentially: after \(n\) periods there are \(2^n\) possible paths. The Bellman equation will allow us to solve this tree efficiently by working backwards from the terminal period, without having to enumerate every possible path.
The worker’s objective is to maximize expected discounted lifetime income:
\[ \mathbb{E} \sum_{t=0}^{T} \beta^t y_t \]
where \(y_t = c\) if unemployed and \(y_t = w_t\) if employed. We initially assume there is no firing (\(\alpha = 0\)), so once a worker accepts a job, they keep it forever (or until they choose to quit).
An important observation simplifies the analysis considerably: an employed worker with wage \(\bar{w}_j\) can quit, collect \(c\) for the current period, and draw a new offer next period. This means that an employed worker deciding whether to keep wage \(\bar{w}_j\) faces exactly the same problem as an unemployed worker who has just received offer \(\bar{w}_j\). We can therefore use a single value function \(V_{t,s}\) to represent the value of having the option of wage \(\bar{w}_s\) with \(t\) periods remaining, regardless of whether the worker is currently employed or unemployed.
2.2 Setup Parameters
We begin by loading the necessary Julia packages and defining the model parameters:
using Plots
using LinearAlgebraβ = 0.95 # discount factor
S = 50 # number of possible wage offers
w̄ = LinRange(1., 10., S) # wage offers (evenly spaced from 1 to 10)
p = ones(S)/S # equal probability of each wage
c = 3. # unemployment benefit3.0
The wage grid w̄ contains 50 evenly spaced values from 1 to 10, and the probability vector p assigns equal probability \(1/S\) to each wage. The discount factor \(\beta = 0.95\) means the worker values a dollar next period at 95 cents today, and the unemployment benefit is \(c = 3\).
3 Finite Horizon
3.1 The One-Period Problem
We begin with the simplest case: a worker with only one period remaining. If the worker receives wage offer \(\bar{w}_j\), they simply choose whichever is larger—the wage or the unemployment benefit:
\[ V_{1,j} = \max\left\{ \bar{w}_j,\; c \right\} \]
Here \(V_{1,j}\) denotes the value with one period left and wage offer \(\bar{w}_j\). In Julia, we can compute this for all wage offers at once using the dot syntax for element-wise operations:
# V1[j] = max(w̄[j], c) for each j
V1 = max.(w̄, c)50-element Vector{Float64}:
3.0
3.0
3.0
3.0
3.0
3.0
3.0
3.0
3.0
3.0
⋮
8.53061224489796
8.714285714285714
8.89795918367347
9.081632653061224
9.26530612244898
9.448979591836736
9.63265306122449
9.816326530612246
10.0
scatter(w̄, V1, label="V₁", legend=:topleft)
xlabel!("Wage")
ylabel!("Value (1 period left)")The plot shows the familiar kinked shape: for wages below \(c = 3\), the worker rejects and receives \(c\); for wages above \(c\), the worker accepts and receives the wage.
We also define \(Q_1\), the expected value of receiving a random wage offer with one period left:
\[ Q_1 = \sum_{s=1}^{S} p_s V_{1,s} = p \cdot V_1 \]
Q1 = dot(p, V1)5.737959183673469
3.2 The Two-Period Problem
Now consider a worker with two periods remaining who receives wage offer \(\bar{w}_s\). The value of accepting is the wage this period plus the discounted continuation value of being employed at that wage next period:
\[ \bar{w}_s + \beta V_{1,s} \]
The value of rejecting is the unemployment benefit this period plus the discounted expected value of a random offer next period:
\[ c + \beta Q_1 \]
The worker chooses whichever is larger:
\[ V_{2,s} = \max\left\{ \bar{w}_s + \beta V_{1,s},\; c + \beta Q_1 \right\} \]
# value of accepting for each wage
V2accept = w̄ + β*V1
# value of rejecting (a scalar, same for all wages)
V2reject = c + β*Q1
V2 = max.(V2accept, V2reject)50-element Vector{Float64}:
8.451061224489795
8.451061224489795
8.451061224489795
8.451061224489795
8.451061224489795
8.451061224489795
8.451061224489795
8.451061224489795
8.451061224489795
8.451061224489795
⋮
16.634693877551022
16.99285714285714
17.351020408163265
17.709183673469386
18.067346938775508
18.425510204081633
18.783673469387757
19.14183673469388
19.5
scatter(w̄, V2, label="V₂", legend=:topleft)
scatter!(w̄, V1, label="V₁")
xlabel!("Wage")
ylabel!("Value")Notice that \(V_2 \geq V_1\) for all wages: having more time remaining is always weakly better, because the worker has more opportunities to find a good wage. We can also compute the expected value of a random offer with two periods left:
Q2 = dot(p, V2)11.970484897959183
The pattern is now clear. We can use \(Q_2\) to solve the three-period problem, \(Q_3\) for the four-period problem, and so on.
3.3 The General \(t\)-Period Problem: The Bellman Equation
For a worker with \(t\) periods remaining and wage offer \(\bar{w}_s\), the same logic applies:
- Value of accepting: \(\bar{w}_s + \beta V_{t-1,s}\)
- Value of rejecting: \(c + \beta Q_{t-1}\)
The worker’s optimal value is:
\[ \boxed{V_{t,s} = \max\left\{ \bar{w}_s + \beta V_{t-1,s},\; c + \beta Q_{t-1} \right\}} \]
This is the Bellman equation for the McCall search model. It expresses the value today as a function of the value tomorrow, capturing the recursive nature of the decision problem.
We solve backwards from the last period: compute \(V_1\), then use it to compute \(V_2\), then \(V_3\), and so on up to \(V_T\). This approach is called backward induction.
We can implement this with a for loop:
T = 50 # number of periods the worker lives
S = length(w̄) # number of possible wages
# V[t,s] = value with t periods left, offer s
V = zeros(T, S)
# Q[t] = expected value of random offer
Q = zeros(T)
# base case: last period, just pick max
V[1, :] = max.(c, w̄)
Q[1] = dot(p, V[1, :])
# build up from t=2 to t=T
for t in 2:T
# wage today + discounted continuation
Vaccept = w̄ + β*V[t-1, :]
# benefits today + discounted random offer
Vreject = c + β*Q[t-1]
# optimal choice for each wage
V[t, :] = max.(Vaccept, Vreject)
# update expected value
Q[t] = dot(p, V[t, :])
endLet us plot the value functions for several time horizons:
plot()
for t in [1, 10, 20, 30, 50]
scatter!(w̄, V[t, :], label="t = $t")
end
xlabel!("Wage")
ylabel!("Value")As the number of remaining periods increases, the value function shifts upward. This reflects the fact that a worker with more time ahead has more opportunities to search for a high wage. Notice, however, that the difference between \(V_{30}\) and \(V_{50}\) is quite small—the value function appears to be converging as \(t\) grows.
4 Solving with Functions
4.1 The Bellman Iteration Function
Rather than repeating the Bellman update code every time we need it, we wrap it in a reusable function. The function iterateBellman takes the previous period’s value function \(V\), expected value \(Q\), and the model parameters, and returns the updated values along with the optimal acceptance policy:
function iterateBellman(V, Q, β, p, w̄, c)
# value of accepting each wage
V_accept = w̄ + β * V
# value of rejecting (same for all wages)
V_reject = c + β * Q
# pick the better option
V_new = max.(V_accept, V_reject)
# expected value of random offer
Q_new = dot(p, V_new)
# 1 = accept, 0 = reject
C = V_accept .>= V_reject
# return as named tuple
return (V=V_new, Q=Q_new, C=C)
enditerateBellman (generic function with 1 method)
Notice that the function returns a named tuple (V=V_new, Q=Q_new, C=C) rather than a plain tuple. Named tuples are a lightweight way to return multiple values with descriptive labels. The caller can access results either by name or by destructuring:
ret = iterateBellman(V, Q, β, p, w̄, c)
ret.V # access by name — self-documenting
ret.C # no need to remember the order# destructuring still works
V, Q, C = iterateBellman(V, Q, β, p, w̄, c)With three or four return values, it is easy to forget the order. Named tuples let you access results by name, making the code self-documenting and less error-prone. We will use this pattern throughout these notes.
4.2 The Finite Horizon Solver
We can now write a complete solver for the finite horizon problem. The function solveMcCall takes the model parameters and the number of periods \(T\), and returns matrices of values and policies:
function solveMcCall(β, p, w̄, c, T)
S = length(w̄)
V = zeros(T, S); Q = zeros(T); C = zeros(Int, T, S)
# base case
V[1, :] = max.(w̄, c)
Q[1] = dot(p, V[1, :])
C[1, :] = w̄ .>= c
for t in 2:T
V[t, :], Q[t], C[t, :] =
iterateBellman(V[t-1, :], Q[t-1], β, p, w̄, c)
end
return (V=V, Q=Q, C=C)
endsolveMcCall (generic function with 1 method)
Let us use the function to solve for 100 periods and examine how the value function evolves:
V, Q, C = solveMcCall(β, p, w̄, c, 100)
plot()
for t in [1, 2, 3, 10, 20, 50, 100]
scatter!(w̄, V[t, :], label="t = $t")
end
xlabel!("Wage")
ylabel!("Value of Offer")An important observation: there is very little difference between \(V_{50}\) and \(V_{100}\). The value function is converging as the horizon grows, which motivates us to consider the infinite horizon case directly.
5 Infinite Horizon
5.1 From Finite to Infinite
As the number of remaining periods \(t\) grows large, the value function converges: \(V_{t,s} \approx V_{t+1,s}\). Intuitively, the problem faced today becomes indistinguishable from the problem faced tomorrow when the horizon is sufficiently long. Taking the limit of the Bellman equation as \(t \to \infty\), the time subscripts drop out:
\[ V_s = \max\left\{ \bar{w}_s + \beta V_s,\; c + \beta Q \right\} \]
where \(Q = p \cdot V\). This is a fixed point equation: the value function \(V\) appears on both sides.
The infinite horizon value function satisfies a fixed point equation. We cannot solve for \(V\) directly by plugging in known quantities from the previous period, because there is no “previous period”—the problem is stationary.
Instead, we solve by value function iteration (VFI): start with an initial guess for \(V\), apply the Bellman equation to get an updated \(V\), and repeat until the function stops changing. This is the computational workhorse for solving dynamic programming problems.
5.2 Value Function Iteration
The algorithm is straightforward:
- Initialize \(V\) with a reasonable guess (e.g., \(V_s = \max(\bar{w}_s, c)\)).
- Apply the Bellman operator to get a new \(V\).
- Check convergence: if \(\|V_{\text{new}} - V\|_\infty < \varepsilon\), stop.
- Otherwise, replace \(V\) with \(V_{\text{new}}\) and go to step 2.
# initial guesses
V = max.(w̄, c)
C = w̄ .>= c
Q = dot(p, V)
dist = 1.0
# loop until convergence
while dist > 1e-10
# one Bellman step
V_new, Q_new, C = iterateBellman(V, Q, β, p, w̄, c)
# max |V_new - V| across wages
dist = norm(V - V_new, Inf)
V = V_new
Q = Q_new
endThe Bellman operator is a contraction mapping with modulus \(\beta\). This means that after \(n\) iterations, the error is at most \(\beta^n \|V_0 - V^*\|\), where \(V^*\) is the true solution. With \(\beta = 0.95\), each iteration shrinks the error by 5%. Convergence is geometric—fast enough for most applications, though models with \(\beta\) very close to 1 may require many iterations.
5.3 The Infinite Horizon Solution
Let us examine the converged value function and the associated optimal policy:
scatter(w̄, V, label="Value", legend=:topleft)
xlabel!("Wage")
ylabel!("Value")scatter(w̄, C, label="Accept (1) / Reject (0)", legend=:topleft)
xlabel!("Wage")
ylabel!("Choice")5.4 The Reservation Wage
The optimal policy has a clean structure: it is a threshold rule. The worker accepts any wage above a cutoff value and rejects wages below it. This cutoff is called the reservation wage.
The optimal policy is a threshold rule: accept any wage above the reservation wage and reject wages below it. The reservation wage is the wage at which the worker is exactly indifferent between accepting and rejecting.
We can compute the reservation wage from our solution:
# first wage where worker accepts
w_res = w̄[findfirst(C .== 1)]
println("Reservation wage: ", round(w_res, digits=2))Reservation wage: 7.98
The worker rejects any offer below this threshold and waits for a better one. The reservation wage reflects the balance between the cost of waiting (forgoing current income) and the benefit of waiting (the chance of receiving a higher offer in the future).
5.5 Wrapping in a Function
We package the infinite horizon solver into a function. Notice that this version of solveMcCall takes only four arguments (no \(T\)), which distinguishes it from the finite horizon version through Julia’s multiple dispatch:
function solveMcCall(β, p, w̄, c)
V = max.(w̄, c); C = w̄ .>= c; Q = dot(p, V)
dist = 1.0
while dist > 1e-10
V_new, Q_new, C = iterateBellman(V, Q, β, p, w̄, c)
dist = norm(V - V_new, Inf)
V = V_new; Q = Q_new
end
return (V=V, Q=Q, C=C)
endsolveMcCall (generic function with 2 methods)
Julia’s multiple dispatch system automatically selects the correct method based on the number and types of arguments: calling solveMcCall(β, p, w̄, c, T) invokes the finite horizon solver, while solveMcCall(β, p, w̄, c) invokes the infinite horizon version.
6 Organizing with Structs
6.1 The Problem with Many Arguments
So far, our solver functions take four or five separate arguments: β, p, w̄, c, and possibly T. This works, but as models grow in complexity—adding firing probability, different wage distributions, or additional features—the number of parameters proliferates. Passing many individual arguments to every function is error-prone (it is easy to mix up the order) and tedious to maintain (changing a parameter requires updating every function call).
The solution is to bundle all model parameters into a single object using a Julia struct.
6.2 Structs as Blueprints
A struct defines a new data type that groups related data together. Think of it as a blueprint or template: the struct declaration specifies what fields the object needs, and creating an instance fills in specific values for those fields.
The analogy is useful: McCallModel is like a blank form that says “every McCall model needs a \(\beta\), \(S\), \(\bar{w}\), \(p\), \(c\), and \(\alpha\).” Creating an instance fills in the blanks with specific numbers. You can create as many different instances as you like, each with different parameter values—this is especially handy for comparative statics.
6.3 The @kwdef Macro
Julia’s built-in @kwdef macro lets us define structs with default values for each field. Combined with Julia’s destructuring syntax (; field1, field2) = struct, we can extract fields into local variables so our code reads cleanly.
6.4 Defining the McCallModel Struct
@kwdef struct McCallModel
β::Float64 = 0.95 # discount factor
S::Int = 50 # number of wage offers
w̄::Vector{Float64} =
collect(LinRange(1., 10., S)) # wage grid
p::Vector{Float64} = ones(S)/S # probabilities (uniform)
c::Float64 = 3.0 # unemployment benefit
α::Float64 = 0.0 # firing probability
endMcCallModel
The @kwdef macro lets every field have a default value. This means you can create a “baseline” model with no arguments at all, or customize only the parameters you want to change:
# all defaults
m = McCallModel()McCallModel(0.95, 50, [1.0, 1.183673469387755, 1.3673469387755102, 1.5510204081632653, 1.7346938775510203, 1.9183673469387754, 2.1020408163265305, 2.2857142857142856, 2.4693877551020407, 2.6530612244897958 … 8.346938775510203, 8.53061224489796, 8.714285714285714, 8.89795918367347, 9.081632653061224, 9.26530612244898, 9.448979591836736, 9.63265306122449, 9.816326530612246, 10.0], [0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02 … 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02], 3.0, 0.0)
# change only c
m_generous = McCallModel(c = 5.0)McCallModel(0.95, 50, [1.0, 1.183673469387755, 1.3673469387755102, 1.5510204081632653, 1.7346938775510203, 1.9183673469387754, 2.1020408163265305, 2.2857142857142856, 2.4693877551020407, 2.6530612244897958 … 8.346938775510203, 8.53061224489796, 8.714285714285714, 8.89795918367347, 9.081632653061224, 9.26530612244898, 9.448979591836736, 9.63265306122449, 9.816326530612246, 10.0], [0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02 … 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02], 5.0, 0.0)
# change β and c
m_patient = McCallModel(β = 0.99, c = 4.0)McCallModel(0.99, 50, [1.0, 1.183673469387755, 1.3673469387755102, 1.5510204081632653, 1.7346938775510203, 1.9183673469387754, 2.1020408163265305, 2.2857142857142856, 2.4693877551020407, 2.6530612244897958 … 8.346938775510203, 8.53061224489796, 8.714285714285714, 8.89795918367347, 9.081632653061224, 9.26530612244898, 9.448979591836736, 9.63265306122449, 9.816326530612246, 10.0], [0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02 … 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02], 4.0, 0.0)
You access fields using dot notation:
m.β0.95
m.c3.0
m.w̄[1:5]5-element Vector{Float64}:
1.0
1.183673469387755
1.3673469387755102
1.5510204081632653
1.7346938775510203
Inside functions, Julia’s destructuring syntax extracts fields into local variables so that mathematical expressions remain clean:
(; β, p, w̄, c) = m
β0.95
Instead of writing m.β, m.p, m.w̄, m.c everywhere in your equations, destructuring gives you clean local variable names. The code looks the same as the mathematical notation.
6.5 Rewriting Functions with Structs
We now rewrite iterateBellman and solveMcCall to accept a McCallModel as their first argument. The function signatures become simpler, and the body uses destructuring to extract the parameters:
function iterateBellman(m::McCallModel, V, Q)
(; β, p, w̄, c) = m
V_accept = w̄ + β * V
V_reject = c + β * Q
V_new = max.(V_accept, V_reject)
Q_new = dot(p, V_new)
C = V_accept .>= V_reject
return (V=V_new, Q=Q_new, C=C)
enditerateBellman (generic function with 2 methods)
The function now takes three arguments instead of six—the model struct carries the parameters.
function solveMcCall(m::McCallModel)
(; β, p, w̄, c) = m
V = max.(w̄, c); C = w̄ .>= c; Q = dot(p, V)
dist = 1.0
while dist > 1e-10
V_new, Q_new, C = iterateBellman(m, V, Q)
dist = norm(V - V_new, Inf)
V = V_new; Q = Q_new
end
return (V=V, Q=Q, C=C)
endsolveMcCall (generic function with 3 methods)
Using the struct version is straightforward:
m = McCallModel()
sol = solveMcCall(m)(V = [158.24988988232073, 158.24988988232073, 158.24988988232073, 158.24988988232073, 158.24988988232073, 158.24988988232073, 158.24988988232073, 158.24988988232073, 158.24988988232073, 158.24988988232073 … 166.93877550863462, 170.61224489635495, 174.28571428407582, 177.9591836717961, 181.63265305951697, 185.3061224472373, 188.97959183495811, 192.65306122267842, 196.32653061039926, 199.9999999981196], Q = 163.42093671832137, C = Bool[0, 0, 0, 0, 0, 0, 0, 0, 0, 0 … 1, 1, 1, 1, 1, 1, 1, 1, 1, 1])
# access value function by name
sol.V[1:5]5-element Vector{Float64}:
158.24988988232073
158.24988988232073
158.24988988232073
158.24988988232073
158.24988988232073
scatter(m.w̄, sol.V, label="Value", legend=:topleft)
xlabel!("Wage")
ylabel!("Value")Julia picks which version of solveMcCall to call based on the argument types:
solveMcCall(m)wherem::McCallModelcalls the struct version.solveMcCall(β, p, w̄, c)calls the original version with separate arguments.
Both coexist and Julia dispatches to the correct method automatically.
7 Hazard Rates and Unemployment Duration
7.1 The Hazard Rate
The hazard rate is the probability of leaving unemployment in a given period. Since the worker accepts wage \(\bar{w}_s\) with probability \(p_s\) whenever \(C_s = 1\), the hazard rate is:
\[ h = \sum_{s=1}^{S} p_s C_s = p \cdot C \]
m = McCallModel()
sol = solveMcCall(m)
# hazard rate
dot(m.p, sol.C)0.23999999999999996
7.2 Expected Unemployment Duration
If the hazard rate is \(h\), the number of periods until the worker finds a job follows a geometric distribution. The expected duration of unemployment is therefore:
\[ \mathbb{E}[\text{duration}] = \frac{1}{h} \]
h = dot(m.p, sol.C)
println("Hazard rate: ", round(h, digits=3))
println("Expected duration: ", round(1/h, digits=1), " periods")Hazard rate: 0.24
Expected duration: 4.2 periods
A higher hazard rate means the worker is less selective and finds a job more quickly. A lower hazard rate means the worker is more patient or more selective, leading to longer unemployment spells but (on average) higher accepted wages.
7.3 Hazard Rate and Unemployment Benefits
How does the hazard rate change as unemployment benefits \(c\) increase? We can compute the hazard rate for a grid of \(c\) values:
cvalues = LinRange(0, 5, 100)
hvalues = zeros(length(cvalues))
for i in 1:length(cvalues)
sol = solveMcCall(McCallModel(c=cvalues[i]))
hvalues[i] = dot(p, sol.C)
end
plot(cvalues, hvalues, linewidth=2, label="Hazard Rate", legend=:topright)
xlabel!("Unemployment Benefits (c)")
ylabel!("Hazard Rate")The relationship is clear: higher unemployment benefits make the worker more selective (raising the reservation wage), which lowers the hazard rate and increases expected unemployment duration. This is a central prediction of search theory and has important implications for the design of unemployment insurance programs.
8 Firing
8.1 Adding Firing Risk
We now extend the model to allow for firing. An employed worker gets fired with probability \(\alpha\) at the end of each period. When fired, the worker immediately draws a new wage offer. This changes the value of accepting a job, because employment is no longer permanent.
Let \(U_t\) denote the value of being unemployed with \(t\) periods left:
\[ U_t = c + \beta Q_{t-1} \]
The value of accepting wage offer \(\bar{w}_s\) now accounts for the risk of being fired:
\[ \bar{w}_s + \beta\left[(1-\alpha) V_{t-1,s} + \alpha U_{t-1}\right] \]
With probability \(1 - \alpha\) the worker keeps the job and has continuation value \(V_{t-1,s}\); with probability \(\alpha\) they are fired and return to unemployment with value \(U_{t-1}\).
The Bellman equation with firing is:
\[ \boxed{V_{t,s} = \max\left\{ \bar{w}_s + \beta\left[(1-\alpha) V_{t-1,s} + \alpha U_{t-1}\right],\; U_t \right\}} \]
Firing risk reduces the value of accepting a job, since employment is no longer permanent. This makes the worker less selective—they are more willing to accept lower wages because even a good job might not last.
8.2 Bellman Iteration with Firing
function iterateBellmanFiring(m::McCallModel, V, Q, U)
(; β, p, w̄, c, α) = m
# keep job w.p. 1-α, fired w.p. α
V_accept = w̄ + β * ((1-α)*V .+ α*U)
# unemployment value
U_new = c + β * Q
# optimal choice
V_new = max.(V_accept, U_new)
# expected value of random offer
Q_new = dot(p, V_new)
# acceptance policy
C = V_accept .>= U_new
return (V=V_new, Q=Q_new, U=U_new, C=C)
enditerateBellmanFiring (generic function with 1 method)
The function signature is clean: the McCallModel carries all parameters including \(\alpha\), so we only need to pass the state variables \(V\), \(Q\), and \(U\).
8.3 Finite Horizon Solver with Firing
function solveMcCallFiring(m::McCallModel, T::Int)
(; β, p, w̄, c, α) = m
S = length(w̄)
V = zeros(T, S); Q = zeros(T); U = zeros(T)
C = zeros(Int, T, S)
V[1, :] = max.(w̄, c); Q[1] = dot(p, V[1, :])
U[1] = c; C[1, :] = w̄ .>= c
for t in 2:T
V[t,:], Q[t], U[t], C[t,:] =
iterateBellmanFiring(m, V[t-1,:], Q[t-1], U[t-1])
end
return (V=V, Q=Q, U=U, C=C)
endsolveMcCallFiring (generic function with 1 method)
8.4 Infinite Horizon Solver with Firing
function solveMcCallFiring(m::McCallModel)
(; β, p, w̄, c, α) = m
V = max.(w̄, c); C = w̄ .>= c
Q = dot(p, V); U = c
dist = 1.0
while dist > 1e-10
V_new, Q_new, U_new, C = iterateBellmanFiring(m, V, Q, U)
dist = norm(V - V_new, Inf)
V = V_new; Q = Q_new; U = U_new
end
return (V=V, Q=Q, U=U, C=C)
endsolveMcCallFiring (generic function with 2 methods)
Again, Julia’s multiple dispatch distinguishes between the finite horizon version solveMcCallFiring(m, T) and the infinite horizon version solveMcCallFiring(m) based on the number of arguments.
9 Comparative Statics
With our solvers in hand, we can now examine how changes in model parameters affect the worker’s value function and acceptance policy.
9.1 Effect of Firing Probability (\(\alpha\))
How does increasing the firing probability affect the value of wage offers?
plot()
for α in LinRange(0., 0.15, 6)
sol = solveMcCallFiring(McCallModel(α=α))
scatter!(w̄, sol.V, label="α = $(round(α, digits=2))")
end
xlabel!("Wage")
ylabel!("Value")Higher firing risk reduces the value of employment (since jobs are less stable), which in turn lowers the overall value function. The worker responds by becoming less selective—accepting lower wages because even an imperfect job provides income while it lasts.
9.2 Effect of Unemployment Benefits (\(c\))
Fixing the firing probability at \(\alpha = 0.02\), we can examine how changes in unemployment benefits affect outcomes:
α = 0.020.02
plot()
for c in LinRange(2., 4., 5)
sol = solveMcCallFiring(McCallModel(α=α, c=c))
scatter!(w̄, sol.V, label="c = $(round(c, digits=1))")
end
xlabel!("Wage")
ylabel!("Value")Higher unemployment benefits increase the value function (since the outside option is more attractive) and raise the reservation wage, making the worker more selective.
9.3 Mean Preserving Spread
A particularly interesting comparative static involves changing the variance of the wage distribution while keeping the mean constant. This is called a mean preserving spread.
phat = 0.04 * ((w̄ .- 5.5) / 4.5).^2
phat .-= dot(phat, w̄) / sum(w̄)
p2 = p + phatplot(w̄, p, label="p (original)", linewidth=2)
plot!(w̄, p2, label="p₂ (more spread)", linewidth=2)
xlabel!("Wage")
ylabel!("Probability")We can verify that both distributions have the same mean:
println("Mean under p: ", dot(p, w̄))
println("Mean under p₂: ", dot(p2, w̄))Mean under p: 5.500000000000001
Mean under p₂: 5.5
The distribution \(p_2\) places more weight on both high and low wages relative to \(p\). How does this affect the worker?
sol = solveMcCallFiring(McCallModel(α=α))
sol2 = solveMcCallFiring(McCallModel(α=α, p=p2))
scatter(w̄, sol.V, label="p = original")
scatter!(w̄, sol2.V, label="p = more spread")
xlabel!("Wage")
ylabel!("Value")scatter(w̄, sol.C, label="p = original")
scatter!(w̄, sol2.C, label="p = more spread")
xlabel!("Wage")
ylabel!("Choice (1 = accept)")The worker is strictly better off with the more dispersed wage distribution, even though the average wage is unchanged. The worker also becomes more selective, raising the reservation wage.
Why does more variance help the worker? Because the worker has the option to reject bad offers. More variance means more very high offers (which the worker accepts) and more very low offers (which the worker rejects at no cost). The ability to reject insulates the worker from downside risk while letting them capture upside risk. This is a manifestation of option value—the value of having the right, but not the obligation, to take an action.
10 Contraction Mapping Theorem
10.1 Why Does Value Function Iteration Converge?
In the infinite horizon section, we claimed that value function iteration converges to the true solution. We also noted in a callout that the Bellman operator is a contraction mapping with modulus \(\beta\). This section makes that claim precise and explores its computational implications.
Define the Bellman operator \(T\) as the mapping that takes a value function \(V\) and produces an updated value function:
\[ (TV)_s = \max\left\{ \frac{\bar{w}_s}{1 - \beta},\; c + \beta \sum_j p_j V_j \right\} \]
Value function iteration is the process of starting with some initial guess \(V_0\) and computing the sequence \(V_1 = TV_0\), \(V_2 = TV_1\), and so on. The question is: does this sequence converge, and if so, does it converge to the right answer?
10.2 Contraction Mappings
A mapping \(T\) on a space of functions is a contraction if it brings any two functions closer together. Formally, \(T\) is a contraction with modulus \(\beta\) if for any two functions \(V\) and \(V'\):
\[ \|TV - TV'\|_\infty \leq \beta \|V - V'\|_\infty \]
where \(\|f\|_\infty = \max_s |f_s|\) is the sup norm (the largest absolute value across all states). The key requirement is that \(\beta < 1\): each application of \(T\) shrinks the distance between any two functions by at least a factor of \(\beta\).
The Bellman operator satisfies this condition because it is a weighted average of future values discounted by \(\beta\). Intuitively, the max operator preserves contractions (you can verify that \(|\max(a, b) - \max(a', b')| \leq \max(|a - a'|, |b - b'|)\)), and the discounting by \(\beta\) provides the contraction.
10.3 The Theorem and Its Implications
The Contraction Mapping Theorem guarantees three things when \(T\) is a contraction:
- Existence and uniqueness: there is exactly one fixed point \(V^*\) satisfying \(TV^* = V^*\).
- Global convergence: starting from any initial guess \(V_0\), the sequence \(T^n V_0\) converges to \(V^*\).
- Geometric convergence rate: after \(n\) iterations, \(\|T^n V_0 - V^*\|_\infty \leq \beta^n \|V_0 - V^*\|_\infty\).
Point 3 is particularly useful in practice. It tells us that the error decreases geometrically at rate \(\beta\). With \(\beta = 0.95\), each iteration shrinks the error by 5%. After 100 iterations, the error is at most \(0.95^{100} \approx 0.006\) times the initial error. After 200 iterations, it is at most \(0.95^{200} \approx 0.00004\) times the initial error.
10.4 Tracking Convergence in Julia
We can verify the geometric convergence rate empirically by tracking \(\|V_{n+1} - V_n\|_\infty\) across iterations. If the theory is correct, this distance should decrease by a factor of approximately \(\beta\) each iteration:
function solveMcCallTracked(m::McCallModel)
(; β, p, w̄, c) = m
V = max.(w̄, c); C = w̄ .>= c; Q = dot(p, V)
distances = Float64[]
dist = 1.0
while dist > 1e-10
V_new, Q_new, C = iterateBellman(m, V, Q)
dist = norm(V - V_new, Inf)
push!(distances, dist)
V = V_new; Q = Q_new
end
return (V=V, Q=Q, C=C, distances=distances)
endsolveMcCallTracked (generic function with 1 method)
m = McCallModel()
sol_tracked = solveMcCallTracked(m)
plot(sol_tracked.distances, yscale=:log10, linewidth=2,
label="‖Vₙ₊₁ - Vₙ‖∞", legend=:topright)
xlabel!("Iteration")
ylabel!("Distance (log scale)")On a log scale, the convergence is nearly a straight line, confirming geometric convergence. The slope of this line corresponds to \(\log_{10}(\beta)\)—steeper for lower \(\beta\) (faster convergence) and flatter for higher \(\beta\) (slower convergence).
We can also compare the convergence rates for different values of \(\beta\):
plot()
for β_val in [0.85, 0.90, 0.95, 0.99]
sol_t = solveMcCallTracked(McCallModel(β = β_val))
plot!(sol_t.distances, yscale=:log10, linewidth=2,
label="β = $β_val")
end
xlabel!("Iteration")
ylabel!("Distance (log scale)")The pattern is clear: higher \(\beta\) means slower convergence. When \(\beta = 0.99\), the Bellman operator barely contracts—each iteration only reduces the error by 1%. This can require hundreds of iterations to reach a tight tolerance. When \(\beta = 0.85\), convergence is rapid, reaching machine precision in a few dozen iterations.
A tolerance of \(10^{-10}\) is standard in most applications. This is far below any economically meaningful threshold, so the numerical solution is effectively exact. For models with \(\beta\) very close to 1 (e.g., \(\beta = 0.999\) at a monthly frequency), you may need to increase the maximum number of iterations or use acceleration techniques like Howard policy iteration, which can converge much faster than plain value function iteration.
11 Model Extensions
The baseline McCall model captures the essential trade-off of job search, but it makes two simplifying assumptions that limit its applicability. First, it assumes that employed workers cannot search for better jobs—once a worker accepts an offer, search stops. Second, it assumes that the worker always receives exactly one offer per period. Relaxing these assumptions leads to richer and more realistic models of labor markets.
11.1 On-the-Job Search
In the baseline model, accepting a job ends the search process. But in reality, employed workers frequently receive outside offers and switch jobs to climb the wage ladder. This section extends the model to allow on-the-job search: employed workers continue to draw offers each period and can switch to any job that pays more than their current wage.
The key change is to the value of employment. An employed worker with wage \(\bar{w}_s\) now draws a new offer \(\bar{w}_j\) each period (with probability \(p_j\)) and takes whichever job pays more:
\[ V_s^{\text{OJS}} = \bar{w}_s + \beta \sum_{j=1}^{S} p_j \max\left( V_j^{\text{OJS}},\; V_s^{\text{OJS}} \right) \]
The worker keeps the current job (\(V_s^{\text{OJS}}\)) unless the new offer is strictly better (\(V_j^{\text{OJS}} > V_s^{\text{OJS}}\)). Since \(V\) is increasing in the wage, this simplifies: the worker accepts any new offer above \(\bar{w}_s\) and rejects any offer below it. Workers ratchet upward through the wage distribution over time.
The value of unemployment is unchanged in structure—the worker still collects \(c\) and draws an offer:
\[ U^{\text{OJS}} = c + \beta \sum_{j=1}^{S} p_j \max\left( V_j^{\text{OJS}},\; U^{\text{OJS}} \right) \]
We can solve this system with a modified version of value function iteration. The Bellman operator for on-the-job search iterates on the value function \(V^{\text{OJS}}\) and the unemployment value \(U^{\text{OJS}}\) simultaneously:
function solveMcCallOJS(m::McCallModel)
(; β, p, w̄, c) = m
S = length(w̄)
# initial guess: value of permanent job at each wage
V = w̄ ./ (1 - β)
U = c / (1 - β)
dist = 1.0
while dist > 1e-10
# employed workers compare current V_s to all V_j
V_new = similar(V)
for s in 1:S
V_new[s] = w̄[s] + β * sum(p[j] * max(V[j], V[s]) for j in 1:S)
end
# unemployed workers compare all V_j to U
U_new = c + β * sum(p[j] * max(V[j], U) for j in 1:S)
dist = max(norm(V - V_new, Inf), abs(U - U_new))
V = V_new
U = U_new
end
# acceptance policy for unemployed: accept if V_s >= U
C = V .>= U
return (V=V, U=U, C=C)
endsolveMcCallOJS (generic function with 1 method)
Let us compare the value functions with and without on-the-job search:
m = McCallModel()
sol_baseline = solveMcCall(m)
sol_ojs = solveMcCallOJS(m)
scatter(m.w̄, sol_baseline.V, label="Baseline", legend=:topleft, alpha=0.7)
scatter!(m.w̄, sol_ojs.V, label="On-the-Job Search", alpha=0.7)
xlabel!("Wage")
ylabel!("Value")On-the-job search raises the value function for all wages because even a mediocre job is now a stepping stone to a better one. The worker can accept a low wage today while continuing to search, rather than being locked in forever. This has two important consequences for the reservation wage:
wstar_baseline = m.w̄[findfirst(sol_baseline.C .== 1)]
wstar_ojs = m.w̄[findfirst(sol_ojs.C .== 1)]
println("Reservation wage (baseline): ", round(wstar_baseline, digits=2))
println("Reservation wage (on-the-job): ", round(wstar_ojs, digits=2))Reservation wage (baseline): 7.98
Reservation wage (on-the-job): 3.02
On-the-job search lowers the reservation wage. This is perhaps surprising at first: if the worker can keep searching after accepting, why not be more selective? The answer is that accepting a low-wage job is no longer costly—the worker can always upgrade later. The option to search on the job reduces the cost of accepting a “wrong” offer, which makes the worker more willing to take the first reasonable job and improve from there.
11.2 Offer Arrival Rate
The baseline model assumes that the worker receives exactly one offer per period. In reality, job offers do not arrive like clockwork—some weeks a worker may receive multiple leads, while other weeks nothing comes through. We can capture this by introducing an offer arrival probability \(\lambda \in (0, 1]\): each period, the worker receives an offer with probability \(\lambda\) and receives no offer (staying unemployed for another period) with probability \(1 - \lambda\).
The Bellman equation for unemployment becomes:
\[ U = c + \beta \left[ \lambda Q + (1 - \lambda) U \right] \]
With probability \(\lambda\) the worker draws an offer and gets expected value \(Q = \sum_s p_s V_s\); with probability \(1 - \lambda\) no offer arrives and the worker remains unemployed with continuation value \(U\). The value of employment is unchanged: an employed worker still earns \(\bar{w}_s\) forever (in the no-firing version).
The modified solver adjusts the unemployed worker’s continuation value to account for the probability of not receiving an offer:
function solveMcCallArrival(m::McCallModel; λ = 1.0)
(; β, p, w̄, c) = m
V = w̄ ./ (1 - β)
U = c / (1 - β)
dist = 1.0
while dist > 1e-10
Q = sum(p[s] * max(V[s], U) for s in 1:length(w̄))
U_new = c + β * (λ * Q + (1 - λ) * U)
V_new = max.(w̄ ./ (1 - β), U_new)
dist = max(norm(V - V_new, Inf), abs(U - U_new))
V = V_new
U = U_new
end
C = (w̄ ./ (1 - β)) .>= U
return (V=V, U=U, C=C)
endsolveMcCallArrival (generic function with 1 method)
How does the arrival rate affect the worker’s behavior? Let us solve the model for several values of \(\lambda\):
m = McCallModel()
plot()
for λ in [0.25, 0.50, 0.75, 1.0]
sol_arr = solveMcCallArrival(m, λ = λ)
scatter!(m.w̄, sol_arr.V, label="λ = $λ", alpha=0.7)
end
xlabel!("Wage")
ylabel!("Value")Lower \(\lambda\) reduces the value function because the worker spends more time waiting without any offer to consider. However, the reservation wage is largely unaffected by \(\lambda\)—conditional on receiving an offer, the accept/reject decision depends on the same trade-off between the flow value of the job and the option value of search. What changes is the expected duration of unemployment. If the hazard rate conditional on receiving an offer is \(h\), then the overall hazard rate including offer arrival is \(\lambda h\), and expected unemployment duration is \(1 / (\lambda h)\).
println("Expected duration by arrival rate:")
for λ in [0.25, 0.50, 0.75, 1.0]
sol_arr = solveMcCallArrival(m, λ = λ)
h = sum(m.p[s] * (sol_arr.C[s] == 1) for s in 1:length(m.w̄))
duration = 1.0 / (λ * h)
println(" λ = $λ: duration = ", round(duration, digits=1), " periods")
endExpected duration by arrival rate:
λ = 0.25: duration = 10.0 periods
λ = 0.5: duration = 6.2 periods
λ = 0.75: duration = 5.1 periods
λ = 1.0: duration = 4.2 periods
Halving the arrival rate approximately doubles the expected unemployment duration, as the worker waits twice as long between offers on average.
The on-the-job search and offer arrival extensions introduced here are building blocks for more advanced models. The Burdett-Mortensen (1998) model combines on-the-job search with wage-posting firms to generate equilibrium wage dispersion. The Diamond-Mortensen-Pissarides (DMP) framework adds a matching function and vacancy creation by firms, producing a full general equilibrium model of unemployment. Both frameworks use the same option-value logic that drives the McCall reservation wage.
12 Summary
These notes developed the McCall search model from the ground up, covering both finite and infinite horizon formulations, the extension to firing risk, and several comparative statics exercises.
Economic Concepts:
- The Bellman equation provides a recursive characterization of the worker’s optimal value function.
- Backward induction solves the finite horizon problem by starting from the last period and working backwards.
- Value function iteration solves the infinite horizon problem by repeatedly applying the Bellman operator until convergence, exploiting the contraction mapping property.
- The optimal policy is a reservation wage rule: accept any offer above the threshold, reject below.
- The hazard rate \(h = p \cdot C\) gives the probability of leaving unemployment each period, and \(1/h\) gives the expected unemployment duration.
- Higher unemployment benefits raise the reservation wage and lower the hazard rate (longer unemployment spells).
- Higher firing risk lowers the reservation wage (workers are less selective when jobs are fragile).
- A mean preserving spread in the wage distribution makes workers strictly better off due to option value.
Julia Tools:
| Tool | Purpose |
|---|---|
max.(x, y) |
Element-wise maximum |
dot(x, y) |
Inner product (expected values) |
norm(x, Inf) |
Maximum absolute value (convergence check) |
@kwdef struct |
Define structs with default values |
(; β, p, w̄, c) = m |
Destructure struct fields into local variables |
Named tuples (V=V, Q=Q, C=C) |
Self-documenting function returns |
| Multiple dispatch | Same function name, different methods based on argument types |
The McCall model shows that search is valuable: a rational worker will reject low offers and wait for better ones, especially when the future is long (\(T\) large) and they are patient (\(\beta\) high). The model provides a rigorous framework for thinking about unemployment duration, the effects of unemployment insurance, and the role of wage dispersion in labor markets.