DiscreteDP Example: Job Search¶

Daisuke Oyama

Faculty of Economics, University of Tokyo

We study an optimal stopping problem, in the context of job search as discussed in http://quant-econ.net/py/lake_model.html.

In [1]:
%matplotlib inline

In [2]:
from __future__ import division, print_function
import numpy as np
import scipy.stats
import scipy.optimize
import scipy.sparse
from numba import jit
import matplotlib.pyplot as plt
from quantecon.markov import DiscreteDP


Optimal solution¶

We skip the description of the model, just writing down the Bellman equation: \begin{aligned} U &= u(c) + \beta \left[(1 - \gamma) U + \gamma E[V_s]\right], \\ V_s &= \max\left\{U, u(w_s) + \beta \left[(1 - \alpha) V_s + \alpha U\right] \right\}. \end{aligned} For this class of problem, we can characterize the solution analytically.

The optimal policy $\sigma^*$ is monotone; it is characterized by a threshold $s^*$, for which $\sigma^*(s) = 1$ if and only if $s \geq s^*$, where actions $0$ and $1$ represent "reject" and "accept", respectively. The threshold is defined as follows: Let \begin{aligned} g(s) &= u(w_s) - u(c), \\ h(s) &= \frac{\beta \gamma}{1 - \beta (1 - \alpha)} \sum_{s' \geq s} p_s u(w_s). \end{aligned} It is easy to see that $g$ is increasing and $h$ is decreasing. Then the threshold $s^*$ is such that $s \geq s^*$ if and only if $g(s) > h(s)$.

Given $s^*$, the optimal values can be computed as follows: \begin{aligned} U &= \frac{\{1 - (1 - \alpha) \beta\} u(c) + \beta \gamma \sum_{s \geq s^*} p_s u(w_s)} {(1 - \beta) \left[\{1 - (1 - \alpha) \beta\} + \beta \gamma \sum_{s \geq s^*} p_s\right]}, \\ V_s &= \begin{cases} U & \text{if s < s^*} \\ \dfrac{u(w_s) + \alpha \beta U}{1 - (1 - \alpha) \beta} & \text{if s \geq s^*}. \end{cases} \end{aligned}

The optimal policy defines a Markov chain over $\{\text{unemployed}, \text{employed}\}$. Its stationary distribution is $\pi = \left(\frac{\alpha}{\alpha + \lambda}, \frac{\lambda}{\alpha + \lambda}\right)$, where $\lambda = \gamma \sum_{s \geq s^*} p(w_s)$; note that the flow from unemployed to employed is $\lambda$, while the flow from employed to unemployed is $\alpha$.

The expected value at the stationary distribution is $$\pi_0 U + \pi_1 \frac{\sum_{s \geq s^*} p_s V_s}{\sum_{s \geq s^*} p_s}.$$

The following implements the job search problem with the analytical solution above:

In [3]:
class JobSearchModel(object):
"""
Job search model.

Parameters
----------
w : array_like(float, ndim=1)
Array containing wage levels. Must be ordered in ascending order.

pdf : array_like(float, ndim=1)
Wage distribution.

beta : scalar(float)
Discount factor

alaph :scalar(float)
Firing probability.

gamma : scalar(float)
Wage offer arrival probability.

rho: scalar(float)
Degree of (constant) relative risk aversion.

"""
def __init__(self, w, pdf, beta, alpha=0, gamma=1, rho=0):
w = np.asarray(w)
self.pdf = np.asarray(pdf)
self.beta, self.alpha, self.gamma, self.rho = beta, alpha, gamma, rho
self.u_w = self.u(w)

def u(self, y):
"""
y must be array_like.

"""
rho = self.rho
small_number = -9999999
y = np.asarray(y, dtype=float)
nonpositive = (y <= 0)
if rho == 1:
util = np.log(y)
else:
util = (y**(1 - rho) - 1)/(1 - rho)
util[nonpositive] = small_number
return util

def solve(self, c, *args, **kwargs):
"""
Solve directly s_star and U and V_s.

"""
S = len(self.u_w)

a0 = 1 - (1 - self.alpha) * self.beta
a1 = self.beta * self.gamma
coeff = a1 / a0
u_c = self.u(np.array([c]))[0]
s_star = _bisect(self.u_w, self.pdf, u_c, coeff)

C = np.zeros(S, dtype=int)
C[s_star:] = 1

U = a0 * u_c + a1 * self.u_w[s_star:].dot(self.pdf[s_star:])
U /= a0 + a1 * self.pdf[s_star:].sum()
U /= 1 - self.beta

V = np.empty(S)
V[:s_star] = U
V[s_star:] = (self.u_w[s_star:] + self.alpha * self.beta * U) / a0

return V, U, C

def stationary_distribution(self, C):
lamb = self.pdf.dot(C) * self.gamma
pi = np.array([self.alpha, lamb])
pi /= pi.sum()
return pi

@jit(nopython=True)
def _bisect(u_w, pdf, u_c, coeff):
lo = -1
hi = len(u_w)
while(lo < hi-1):
m = (lo + hi) // 2
lhs = u_w[m] - u_c
rhs = 0
for i in range(m+1, len(u_w)):
rhs += (u_w[i] - u_w[m]) * pdf[i]
rhs *= coeff
if lhs > rhs:
hi = m
else:
lo = m
return hi


For comparison, let us also consider the implementation with the DiscreteDP class:

In [4]:
class JobSearchModelDiscreteDP(JobSearchModel):
"""
Job search model with DiscreteDP.

"""
def __init__(self, w, pdf, beta, alpha=0, gamma=1, rho=0):
super(JobSearchModelDiscreteDP, self).__init__(w, pdf, beta, alpha, gamma, rho)

# Number of states
# s = 0, ..., len(w)-1: wage w[s] offered, s = len(w): no offer
num_states = len(w) + 1

# Number of actions: 0: reject, 1: accept
num_actions = 2

L = num_states*num_actions - 1
s_indices, a_indices = np.empty(L), np.empty(L)
s_indices[-1], a_indices[-1] = len(w), 0
s_indices[:-1] = np.repeat(np.arange(len(w)), num_actions)
a_indices[:-1] = np.tile(np.arange(num_actions), len(w))

R0 = np.zeros(L)
R0[[num_actions*i+1 for i in range(len(w))]] = self.u_w

Q = scipy.sparse.lil_matrix((L, num_states))
it = np.nditer((s_indices, a_indices))
for (s, a) in it:
i = it.iterindex
if a == 0:
Q[i, -1] = 1 - self.gamma
Q[i, :len(w)] = self.pdf*self.gamma
else:  # if a == 1
Q[i, s], Q[i, -1] = 1 - self.alpha, self.alpha

self.ddp = DiscreteDP(R0, Q, beta, s_indices, a_indices)

self.num_iter = None

def solve(self, c, *args, **kwargs):
n, m = self.ddp.num_states, self.ddp.num_actions
self.ddp.R[[m*i for i in range(n)]] = self.u(np.array([c]))[0]
res = self.ddp.solve(*args, **kwargs)
V = res.v[:-1]  # Values of jobs
U = res.v[-1]  # Value of unemployed
C = res.sigma[:-1]
self.num_iter = res.num_iter

return V, U, C


The following paramter values are from lakemodel_example.py.

In [5]:
w = np.linspace(0, 175, 201)  # wage grid

# compute probability of each wage level
logw_dist = scipy.stats.norm(np.log(20.),1)
cdf = logw_dist.cdf(np.log(w))
pdf = cdf[1:]-cdf[:-1]
pdf /= pdf.sum()
w = (w[1:] + w[:-1])/2

In [6]:
gamma = 1
alpha = 0.013  # Monthly
alpha_q = (1-(1-alpha)**3)  # Quarterly
beta = 0.99
rho = 2  # risk-aversion

In [7]:
js = JobSearchModel(w, pdf, beta, alpha_q, gamma, rho)

In [8]:
js_ddp = JobSearchModelDiscreteDP(w, pdf, beta, alpha_q, gamma, rho)


Let us check that the results coincide:

In [9]:
cs = np.linspace(1, 75, 25)
bools = []
for c in cs:
V, U, C = js.solve(c=c)
V1, U1, C1 = js_ddp.solve(c=c)
bools.append(np.allclose(V, V1))
bools.append(np.allclose(U, U1))
bools.append(np.array_equal(C, C1))
print(all(bools))

True


Take a look at the optimal solution for $c = 40$ for example:

In [10]:
c = 40
V, U, C = js.solve(c=c)

In [11]:
s_star = len(w) - C.sum()
print(r"Optimal policy: Accept if and only if w >= {0}".format(w[s_star]))

Optimal policy: Accept if and only if w >= 65.1875

In [12]:
fig, ax = plt.subplots(figsize=(8,5))
ax.plot(w, V, label=r'$V$')
ax.plot((w[0], w[-1]), (U, U), 'r--', label=r'$U$')
ax.set_xlabel('Wage')
ax.set_ylabel('Value')
ax.set_title('Optimal value function')
plt.legend(loc=2)
plt.show()


Performance comparison:

In [13]:
c = 40
%timeit js.solve(c=c)

10000 loops, best of 3: 58.7 Âµs per loop

In [14]:
%timeit js_ddp.solve(c=c)

100 loops, best of 3: 8.04 ms per loop


Optimal unemployment insurance policy¶

We compute the optimal level of unemployment insurance as in the lecture, mimicking lakemodel_example.py.

In [15]:
class UnemploymentInsurancePolicy(object):
def __init__(self, w, pdf, beta, alpha=0, gamma=1, rho=0):
self.w, self.pdf, self.beta, self.alpha, self.gamma, self.rho = \
w, pdf, beta, alpha, gamma, rho

def solve_job_search_model(self, c, T):
js = JobSearchModel(self.w-T, self.pdf, self.beta,
self.alpha, self.gamma, self.rho)
V, U, C = js.solve(c=c-T)
pi = js.stationary_distribution(C)

return V, U, C, pi

def implement(self, c):

def budget_balance(T):
_, _, _, pi = self.solve_job_search_model(c, T)
return T - pi[0]*c

# Budget balancing tax given c
T = scipy.optimize.brentq(budget_balance, 0, c)

V, U, C, pi = self.solve_job_search_model(c, T)

EV = (C*V).dot(self.pdf)/(C.dot(self.pdf))
W = pi[0] * U + pi[1] * EV

return T, W, pi

In [16]:
uip = UnemploymentInsurancePolicy(w, pdf, beta, alpha_q, gamma, rho)

In [17]:
grid_size = 501  #25
cvec = np.linspace(5, 135, grid_size)
Ts, Ws = np.empty(grid_size), np.empty(grid_size)
pis = np.empty((grid_size, 2))
for i, c in enumerate(cvec):
T, W, pi = uip.implement(c=c)
Ts[i], Ws[i], pis[i] = T, W, pi
i_max = Ws.argmax()
print('Optimal unemployment benefit:', cvec[i_max])

Optimal unemployment benefit: 67.4

In [18]:
def plot(ax, y_vec, title):
ax.plot(cvec, y_vec)
ax.set_xlabel(r"$c$")
ax.vlines(cvec[i_max], ax.get_ylim()[0], y_vec[i_max], "k", "-.")
ax.set_title(title)

fig, axes = plt.subplots(2, 2, figsize=(10, 6))
plot(axes[0, 0], Ws, "Welfare")
plot(axes[0, 1], Ts, "Taxes")
plot(axes[1, 0], pis[:, 1], "Employment Rate")
plot(axes[1, 1], pis[:, 0], "Unemployment Rate")

plt.tight_layout()
plt.show()

In [19]:
fig, ax = plt.subplots()
plot(ax, cvec-Ts, "Net Compensation")
plt.show()

In [ ]: