An unemployed worker receives in each period a job offer at wage $ W_t $.
At time $ t $, our worker has two choices:
The wage sequence $ \{W_t\} $ is assumed to be iid with probability mass function $ p_1, \ldots, p_n $.
Here $ p_i $ is the probability of observing wage offer $ W_t = w_i $ in the set $ w_1, \ldots, w_n $.
The worker is infinitely lived and aims to maximize the expected discounted sum of earnings.
$$ \mathbb{E} \sum_{t=0}^{\infty} \beta^t Y_t $$The discount factor is $ \beta $ lies in $ (0, 1) $. The smaller is $ \beta $, the more the worker discounts future utility relative to current utility.
The variable $ Y_t $ is income, equal to
he worker faces a trade-off:
Let $ V $ be a function that assigns to each possible wage $ w $ the maximal lifetime value that can be obtained with that offer in hand.
$ V $ must satisfy the Bellman Equation:
$$ V(w) = \max \left\{ \frac{w}{1 - \beta}, \, c + \beta \sum_{i=1}^n V(w_i) p_i \right\} \tag{1} $$
for every possible $ w_i $ in $ w_1, \ldots, w_n $, where:
Given any $ w $, we can read off the corresponding best choice (accept or reject) by picking the max on the r.h.s. of (1).
Thus, we have a map from $ \mathbb{R} $ to $ \{0, 1\} $, with 1 meaning accept and zero meaning reject.
We can write the policy as follows
$$ \sigma(w) := \mathbf{1} \left\{ \frac{w}{1 - \beta} \geq c + \beta \sum_{i=1}^n V(w_i) p_i \right\} $$Here $ \mathbf{1}\{ P \} = 1 $ if statement $ P $ is true and equals zero otherwise.
We can also write this as
$$ \sigma(w) := \mathbf{1} \{ w \geq \bar w \} $$where
$$ \bar w := (1 - \beta) \left\{ c + \beta \sum_{i=1}^n V(w_i) p_i \right\} \tag{2} $$
Here $ \bar w $ is a constant depending on $ \beta, c $ and the wage distribution, called the reservation wage.
The agent should accept if and only if the current wage offer exceeds the reservation wage.
Clearly, we can compute this reservation wage if we can compute the value function.
Define that $v_i:=V(w_i)$, we proceed as follows:
Step 1: pick an arbitrary initial guess $ v \in \mathbb R^n $.
Step 2: compute a new vector $ v' \in \mathbb R^n $ via
$$ v'_i = \max \left\{ \frac{w_i}{1 - \beta}, \, c + \beta \sum_{i=1}^n v_i p_i \right\} \quad \text{for } i = 1, \ldots, n \tag{4} $$
Step 3: calculate a measure of the deviation between $ v $ and $ v' $, such as $ \max_i |v_i - v_i'| $.
Step 4: if the deviation is larger than some fixed tolerance, set $ v = v' $ and go to step 2, else continue.
Step 5: return $ v $.
This algorithm returns an arbitrarily good approximation to the true solution to (3), which represents the value function.
(Arbitrarily good means here that the approximation converges to the true solution as the tolerance goes to zero)
We can show that the mapping $T$ from $\mathbb{R}^n$ to itself defined as: $$ Tv_i = \max \left\{ \frac{w_i}{1 - \beta}, \, c + \beta \sum_{i=1}^n v_i p_i \right\} \quad \text{for } i = 1, \ldots, n \tag{5} $$ is a contraction mapping.
#set up (This takes relatively long time for first compiling)
using LinearAlgebra, Statistics
using Distributions, Expectations, NLsolve, Roots, Random, Plots, Parameters, StatsPlots
┌ Info: Precompiling Distributions [31c24e10-a181-5473-b8eb-7969acd0382f] └ @ Base loading.jl:1278 ┌ Info: Precompiling Expectations [2fe49d83-0758-5602-8f54-1f90ad0d522b] └ @ Base loading.jl:1278 ┌ Info: Precompiling NLsolve [2774e3e8-f4cf-5e23-947b-6d7e65073b56] └ @ Base loading.jl:1278 ┌ Info: Precompiling Plots [91a5bcdd-55d7-5caf-9e0b-520d859cae80] └ @ Base loading.jl:1278 ┌ Info: Precompiling StatsPlots [f3b207a7-027a-5e70-b257-86293d7955fd] └ @ Base loading.jl:1278
#set the graph format
gr(fmt = : svg) #png, jpg, etc.
syntax: whitespace not allowed after ":" used for quoting
Stacktrace:
[1] top-level scope at In[2]:2
[2] include_string(::Function, ::Module, ::String, ::String) at .\loading.jl:1091
We assume that the wage offer is drawn from a Beta-Binomial Distribution.
The beta-binomial distribution is the binomial distribution in which the probability of success at each of n trials is not fixed but randomly drawn from a beta distribution.
# The distribution of wage offers
n = 50
dist = BetaBinomial(n, 200, 100) # probability distribution Beta-Binomial(n, α,β)
@show support(dist)
w = range(10.0, 60.0, length = n+1) # linearly space wages
support(dist) = 0:50
10.0:1.0:60.0
#plots of density function
plt = plot(w, pdf.(dist, support(dist)), xlabel = "wages", lw = 2, ylabel = "probabilities",label= "α=200")
# changing parameters
plot!(w,pdf.(BetaBinomial(n,100,100), support(BetaBinomial(n,100,100))), lw = 2, label = "α=100")
plot!(w,pdf.(BetaBinomial(n,10,100), support(BetaBinomial(n,100,100))), lw = 2, label = "α=10")
plot!(w,pdf.(BetaBinomial(n,1,100), support(BetaBinomial(n,100,100))), lw = 2, label = "α=1")
┌ Info: Precompiling GR_jll [d2c73de3-f751-5644-a686-071e5b155ba9] └ @ Base loading.jl:1278
We can use the function expectation
to calculate moments.
#expectations and variances
E = expectation(dist) # operator of expectation
# exploring the properties of the operator
wage(i) = w[i+1] # +1 to map from support of 0
E_w = E(wage)
E_w_2 = E(i -> wage(i)^2) - E_w^2 # variance
@show E_w, E_w_2
# use operator with left-multiply
@show E * w # the `w` are values assigned for the discrete states
@show dot(pdf.(dist, support(dist)), w); # identical calculation
(E_w, E_w_2) = (43.33333333333569, 12.919896640724573) E * w = 43.33333333333569 dot(pdf.(dist, support(dist)), w) = 43.33333333333569
Our intial guess $v$ is the value of accepting at every given wage. Below we show the sequence of value functions after 5 iterations.
# set parameters and constant objects
c = 25
β = 0.99
num_plots = 6
# Operator
T(v) = max.(w/(1 - β), c + β * E*v) # (5) broadcasts over the w, fixes the v
# alternatively, T(v) = [max(wval/(1 - β), c + β * E*v) for wval in w]
# fill in matrix of vs
vs = zeros(n + 1, 6) # data to fill
vs[:, 1] .= w / (1-β) # initial guess of "accept all"
# manually applying operator
for col in 2:num_plots
v_last = vs[:, col - 1]
vs[:, col] .= T(v_last) # apply operator
end
plot(vs)
Now we implement the iteration untile the deviation is quite low.
# iteration function for solving the DP problem (direct implementation)
function compute_reservation_wage_direct(params; v_iv = collect(w ./(1-β)), max_iter = 500,
tol = 1e-6)
@unpack c, β, w = params
# create a closure for the T operator
T(v) = max.(w/(1 - β), c + β * E*v) # (5) fixing the parameter values
v = copy(v_iv) #(important!!) start at initial value. copy to prevent v_iv modification
v_next = similar(v)
i = 0
error = Inf
while i < max_iter && error > tol
v_next .= T(v) # (4)
error = norm(v_next - v)
i += 1
v .= v_next # copy contents into v. Also could have used v[:] = v_next
end
# now compute the reservation wage
return (1 - β) * (c + β * E*v) # (2)
end
compute_reservation_wage_direct (generic function with 1 method)
# Compute the reservation wage
mcm = @with_kw (c=25.0, β=0.99, w=w) # named tuples
compute_reservation_wage_direct(mcm()) # call with default parameters
47.316499757360674
We can also use the fixedpoint
algorithm in the NLsolve
package.
# Example of using fixedpoint
p = 1.0
β1 = 0.9
f(v) = p .+ β1 * v # broadcast the +
sol = fixedpoint(f, [0.8, 0.1],iterations = 100, ftol = 1e-6, m = 6)
println("Fixed point = $(sol.zero), and |f(x) - x| = $(norm(f(sol.zero) - sol.zero)) in " *
"$(sol.iterations) iterations")
Fixed point = [9.999999999999986, 9.999999999999986], and |f(x) - x| = 2.5121479338940403e-15 in 3 iterations
# using the fixedpoint function in the NLsolve package
function compute_reservation_wage(params; v_iv = 1.1*collect(w ./(1-β)), iterations = 500,
ftol = 1e-6, m = 6)
@unpack c, β, w = params
T(v) = max.(w/(1 - β), c + β * E*v) # (5) fixing the parameter values
v_star = fixedpoint(T, v_iv, iterations = iterations, ftol = ftol, m = m).zero # (5)
return (1 - β) * (c + β * E*v_star) # (3)
end
compute_reservation_wage (generic function with 1 method)
# Compute the reservation wage
mcm = @with_kw (c=25.0, β=0.99, w=w) # named tuples
compute_reservation_wage(mcm() )# call with default parameters
47.31649976654636
Note: The convergence is very senstive to the initial values.
#starting at different values can cause the iterations not to converge
#gmethod(5) not working, method 1-4 can work
compute_reservation_wage(mcm(); v_iv = collect(w./(1-β)), m=5)# call with default parameters
LAPACKException(5)
Stacktrace:
[1] chklapackerror at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.5\LinearAlgebra\src\lapack.jl:38 [inlined]
[2] trtrs!(::Char, ::Char, ::Char, ::SubArray{Float64,2,Array{Float64,2},Tuple{UnitRange{Int64},UnitRange{Int64}},false}, ::SubArray{Float64,1,Array{Float64,1},Tuple{UnitRange{Int64}},true}) at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.5\LinearAlgebra\src\lapack.jl:3396
[3] ldiv! at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.5\LinearAlgebra\src\triangular.jl:767 [inlined]
[4] anderson_(::OnceDifferentiable{Array{Float64,1},Array{Float64,2},Array{Float64,1}}, ::Array{Float64,1}, ::Float64, ::Float64, ::Int64, ::Bool, ::Bool, ::Bool, ::Int64, ::Int64, ::Int64, ::NLsolve.AndersonCache{Array{Float64,1},Array{Float64,1},Array{Array{Float64,1},1},Array{Float64,1},Array{Float64,2},Array{Float64,2}}) at C:\Users\Zhiyuan Chen\.julia\packages\NLsolve\gJL1I\src\solvers\anderson.jl:160
[5] anderson(::OnceDifferentiable{Array{Float64,1},Array{Float64,2},Array{Float64,1}}, ::Array{Float64,1}, ::Float64, ::Float64, ::Int64, ::Bool, ::Bool, ::Bool, ::Int64, ::Int64, ::Int64, ::NLsolve.AndersonCache{Array{Float64,1},Array{Float64,1},Array{Array{Float64,1},1},Array{Float64,1},Array{Float64,2},Array{Float64,2}}) at C:\Users\Zhiyuan Chen\.julia\packages\NLsolve\gJL1I\src\solvers\anderson.jl:203
[6] anderson(::OnceDifferentiable{Array{Float64,1},Array{Float64,2},Array{Float64,1}}, ::Array{Float64,1}, ::Float64, ::Float64, ::Int64, ::Bool, ::Bool, ::Bool, ::Int64, ::Int64, ::Int64, ::Int64) at C:\Users\Zhiyuan Chen\.julia\packages\NLsolve\gJL1I\src\solvers\anderson.jl:188
[7] nlsolve(::OnceDifferentiable{Array{Float64,1},Array{Float64,2},Array{Float64,1}}, ::Array{Float64,1}; method::Symbol, xtol::Float64, ftol::Float64, iterations::Int64, store_trace::Bool, show_trace::Bool, extended_trace::Bool, linesearch::Static, linsolve::NLsolve.var"#27#29", factor::Float64, autoscale::Bool, m::Int64, beta::Int64, aa_start::Int64, droptol::Int64) at C:\Users\Zhiyuan Chen\.julia\packages\NLsolve\gJL1I\src\nlsolve\nlsolve.jl:30
[8] fixedpoint(::var"#T#46"{Float64,Float64,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}, ::Array{Float64,1}; method::Symbol, xtol::Float64, ftol::Float64, iterations::Int64, store_trace::Bool, show_trace::Bool, extended_trace::Bool, linesearch::Static, factor::Float64, autoscale::Bool, m::Int64, beta::Int64, droptol::Int64, autodiff::Symbol, inplace::Bool) at C:\Users\Zhiyuan Chen\.julia\packages\NLsolve\gJL1I\src\nlsolve\fixedpoint.jl:34
[9] compute_reservation_wage(::NamedTuple{(:c, :β, :w),Tuple{Float64,Float64,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}}; v_iv::Array{Float64,1}, iterations::Int64, ftol::Float64, m::Int64) at .\In[80]:7
[10] top-level scope at In[88]:3
[11] include_string(::Function, ::Module, ::String, ::String) at .\loading.jl:1091
#starting at different values can cause the iterations not to converge
compute_reservation_wage(mcm(); v_iv = 0.9*collect(w./(1-β)))# call with default parameters
47.316499766546144
We want to know how the reservation wage differs for different parameters $(\beta, c)$
Intuition:
grid_size = 50
R = rand(grid_size, grid_size)
c_vals = range(10.0, 30.0, length = grid_size)
β_vals = range(0.9, 0.99, length = grid_size)
for (i, c) in enumerate(c_vals)
for (j, β) in enumerate(β_vals)
R[i, j] = compute_reservation_wage(mcm(c=c, β=β); v_iv = collect(w ./(1-β)), m =1) # change from defaults
end
end
contour(c_vals, β_vals, R',
title = "Reservation Wage",
xlabel = "c",
ylabel = "β",
fill = true)
surface(c_vals, β_vals, R,
title = "Reservation Wage",
xlabel = "c",
ylabel = "β",
zlabel = "w̄")
Let $ \psi $ denote the value of not accepting a job in this period but then behaving optimally in all subsequent periods.
That is,
$$ \psi = c + \beta \sum_{i=1}^n V(w_i) p_i \tag{6} $$
where $ V $ is the value function.
By the Bellman equation, we then have
$$ V(w_i) = \max \left\{ \frac{w_i}{1 - \beta}, \, \psi \right\} $$Substituting this last equation into (6) gives
$$ \psi = c + \beta \sum_{i=1}^n \max \left\{ \frac{w_i}{1 - \beta}, \psi \right\} p_i \tag{7} $$
Which we could also write as $ \psi = T_{\psi}(\psi) $ for the appropriate operator.
This is a nonlinear equation that we can solve for $ \psi $.
One solution method for this kind of nonlinear equation is iterative.
That is,
Step 1: pick an initial guess $ \psi $.
Step 2: compute the update $ \psi' $ via
$$ \psi' = c + \beta \sum_{i=1}^n \max \left\{ \frac{w_i}{1 - \beta}, \psi \right\} p_i \tag{8} $$
Step 3: calculate the deviation $ |\psi - \psi'| $.
Step 4: if the deviation is larger than some fixed tolerance, set $ \psi = \psi' $ and go to step 2, else continue.
Step 5: return $ \psi $.
Once again, one can use the Banach contraction mapping theorem to show that this process always converges.
The big difference here, however, is that we’re iterating on a single number, rather than an $ n $-vector.
Here’s an implementation using two methods:
# fixed point iteration
function compute_reservation_wage_ψ(c, β; ψ_iv = E * w ./ (1 - β), max_iter = 500,
tol = 1e-6)
T_ψ(ψ) = [c + β * E*max.((w ./ (1 - β)), ψ[1])] # (7)
# using vectors since fixedpoint doesn't support scalar
ψ_star = fixedpoint(T_ψ, [ψ_iv]).zero[1]
return (1 - β) * ψ_star # (2)
end
compute_reservation_wage_ψ(c, β)
47.316499766546116
# root finding
function compute_reservation_wage_ψ2(c, β; ψ_iv = E * w ./ (1 - β), max_iter = 500,
tol = 1e-6)
root_ψ(ψ) = c + β * E*max.((w ./ (1 - β)), ψ) - ψ # (7)
ψ_star = find_zero(root_ψ, ψ_iv)
return (1 - β) * ψ_star # (2)
end
compute_reservation_wage_ψ2(c, β)
47.316499766546144
We are also interested in the duration of unemployment in the economy. Now we compute the average duration of unemployment for different $c$s, which in reality can be affected by the government.
Simulation Method:
function stopping_time(w̄; seed = 1212)
Random.seed!(seed)
t = 1
# @assert length(w)-1 in support(dist) && ̄w <=w[end]
end
stopping_time (generic function with 1 method)
function compute_stopping_time(w̄; seed=1234)
Random.seed!(seed)
stopping_time = 0
t = 1
# make sure the constraint is sometimes binding
@assert length(w) - 1 ∈ support(dist) && w̄≤ w[end]
while true
# Generate a wage draw
w_val = w[rand(dist)] # the wage dist set up earlier
if w_val ≥ w̄
stopping_time = t
break
else
t += 1
end
end
return stopping_time
end
compute_stopping_time (generic function with 1 method)
compute_mean_stopping_time(w̄, num_reps=10000) = mean(i ->
compute_stopping_time(w̄,
seed = i), 1:num_reps)
c_vals = range(10, 40, length = 25)
stop_times = similar(c_vals)
beta = 0.99
for (i, c) in enumerate(c_vals)
w̄ = compute_reservation_wage_ψ(c, beta)
stop_times[i] = compute_mean_stopping_time(w̄)
end
p1 = plot(c_vals, stop_times, label = "mean unemployment duration: β=0.99",
xlabel = "c", ylabel = "months", lw = 3);
beta = 0.9
for (i, c) in enumerate(c_vals)
w̄ = compute_reservation_wage_ψ(c, beta)
stop_times[i] = compute_mean_stopping_time(w̄)
end
p2 = plot(c_vals, stop_times,label = "mean unemployment duration:β=0.9",
xlabel = "c", ylabel = "months", lw = 3);
plot(p1, p2, layout=(1,2))
Increasing the unemployment compensation will lead to the average length of unemployment duration to grow. The agent will wait for a sufficiently good wage draw before getting to work.
As the agent becomes less patient, the average duration of unemployment decreases and will be more sensitive to the unemployment compensation.
function compute_std_stopping_time(w̄, num_reps=10000)
stopping_time_vec = [compute_stopping_time(w̄, seed = i) for i in 1:num_reps]
return var_stopping_time = std(stopping_time_vec)
end
c_vals = range(10, 40, length = 25)
stop_times = similar(c_vals)
beta = 0.99
for (i, c) in enumerate(c_vals)
w̄ = compute_reservation_wage_ψ(c, beta)
stop_times[i] = compute_std_stopping_time(w̄)
end
p1 = plot(c_vals, stop_times, label = "std unemployment duration: β=0.99",
xlabel = "c", ylabel = "months", lw = 3);
beta = 0.9
for (i, c) in enumerate(c_vals)
w̄ = compute_reservation_wage_ψ(c, beta)
stop_times[i] = compute_std_stopping_time(w̄)
end
p2 = plot(c_vals, stop_times,label = "std unemployment duration:β=0.9",
xlabel = "c", ylabel = "months", lw = 3);
plot(p1, p2, layout=(1,2))
stop_times = similar(c_vals)
beta = 0.99
w̄ = compute_reservation_wage_ψ(10, beta)
stopping_time_vec = [compute_stopping_time(w̄, seed = i) for i in 1:10000]
p1 = histogram(stopping_time_vec,label = "c= 10")
w̄ = compute_reservation_wage_ψ(20, beta)
stopping_time_vec = [compute_stopping_time(w̄, seed = i) for i in 1:10000]
p2 = histogram(stopping_time_vec,label = "c= 20")
w̄ = compute_reservation_wage_ψ(40, beta)
stopping_time_vec = [compute_stopping_time(w̄, seed = i) for i in 1:10000]
p3 = histogram(stopping_time_vec,label = "c= 40")
plot(p1,p2, p3, layout=(3, 1))