text: 
  import numpy as np import scipy.stats as stats from interpolation import interp from numba import njit, prange import matplotlib.pyplot as plt %matplotlib inline from math import gamma
  Model¶
  - $
  $$ x_{t+1} = (1 - b_{t+1}) g(x_t, \phi_t) + b_{t+1} \max \{ g(x_t, \phi_t), u_{t+1}\} \tag{1} $$
  Agent’s objective: maximize expected discounted sum of wages via controls $ \{s_t\} $ and $ \{\phi_t\} $.
  Taking the expectation of $ v(x_{t+1}) $ and using (1), the Bellman equation for this problem can be written as
  $$ v(x) = \max_{s + \phi \leq 1} \left\{ x (1 - s - \phi) + \beta (1 - \pi(s)) v[g(x, \phi)] + \beta \pi(s) \int v[g(x, \phi) \vee u] f(du) \right\} \tag{2} $$
  Here nonnegativity of $ s $ and $ \phi $ is understood, while $ a \vee b := \max\{a, b\} $.
  Parameterization¶
  In the implementation below, we will focus on the parameterization$$ g(x, \phi) = A (x \phi)^{\alpha}, \quad \pi(s) = \sqrt s \quad \text{and} \quad f = \text{Beta}(2, 2) $$
  with default parameter values
  - $ A = 1.4 $
  - $ \alpha = 0.6 $
  - $ \beta = 0.96 $
  The $ \text.
  class JVWorker: r""" A Jovanovic-type model of employment with on-the-job search. """ def __init__(self, A=1.4, α=0.6, β=0.96, # Discount factor π=np.sqrt, # Search effort function a=2, # Parameter of f b=2, # Parameter of f grid_size=50, mc_size=100, ɛ=1e-4): self.A, self.α, self.β, self.π = A, α, β, π self.mc_size, self.ɛ = mc_size, ɛ self.g = njit(lambda x, ϕ: A * (x * ϕ)**α) # Transition function self.f_rvs = np.random.beta(a, b, mc_size) # Max of grid is the max of a large quantile value for f and the # fixed point y = g(y, 1) ɛ = 1e-4 grid_max = max(A**(1 / (1 - α)), stats.beta(a, b).ppf(1 - ɛ)) # Human capital self.x_grid = np.linspace(ɛ, grid_max, grid_size)
  The function
  operator_factory takes an instance of this class and returns a
  jitted version of the Bellman operator
  T, ie.
  where
  $$ w(s, \phi) := x (1 - s - \phi) + \beta (1 - \pi(s)) v[g(x, \phi)] + \beta \pi(s) \int v[g(x, \phi) \vee u] f(du) \tag{3} $$
  When we represent $ v $, it will be with a NumPy array
  v giving values on grid
  x_grid.
  But to evaluate the right-hand side of (3), we need a function, so
  we replace the arrays
  v and
  x_grid with a function
  v_func that gives linear
  interpolation of
  v on
  x_grid.
  Inside the
  for loop, for each
  x in the grid over the state space, we
  set up the function $ w(z) = w(s, \phi) $ defined in (3).
  The function is maximized over all feasible $ (s, \phi) $ pairs.
  Another function,
  get_greedy returns the optimal policies of
  s and $ \phi $
  given a value function.
  def operator_factory(jv, parallel_flag=True): """ Returns a jitted version of the Bellman operator T jv is an instance of JVWorker """ π, β = jv.π, jv.β x_grid, ɛ, mc_size = jv.x_grid, jv.ɛ, jv.mc_size f_rvs, g = jv.f_rvs, jv.g @njit def objective(z, x, v): s, ϕ = z v_func = lambda x: interp(x_grid, v, x) integral = 0 for m in range(mc_size): u = f_rvs[m] integral += v_func(max(g(x, ϕ), u)) integral = integral / mc_size q = π(s) * integral + (1 - π(s)) * v_func(g(x, ϕ)) return x * (1 - ϕ - s) + β * q @njit(parallel=parallel_flag) def T(v): """ The Bellman operator """ v_new = np.empty_like(v) for i in prange v_new[i] = max_val return v_new @njit def get_greedy(v): """ Computes the v-greedy policy of a given function v """ s_policy, ϕ_policy = np.empty_like(v), np.empty_like(v) for i in range max_s, max_ϕ = s, ϕ s_policy[i], ϕ_policy[i] = max_s, max_ϕ return s_policy, ϕ_policy return T, get_greedy
  To solve the model, we will write a function that uses the Bellman operator and iterates to find a fixed point.
  def solve_model(jv, use_parallel=True, tol=1e-4, max_iter=1000, verbose=True, print_skip=25): """ Solves the model by value function iteration * jv is an instance of JVWorker """ T, _ = operator_factory(jv, parallel_flag=use_parallel) # Set up loop v = jv.x_grid * 0.5 #
  jv = JVWorker() T, get_greedy = operator_factory(jv) v_star = solve_model(jv) s_star, ϕ_star = get_greedy(v_star) plots = [s_star, ϕ_star, v_star] titles = ["ϕ policy", "s policy", "value function"] fig, axes = plt.subplots(3, 1, figsize=(12, 12)) for ax, plot, title in zip(axes, plots, titles): ax.plot(jv.x_grid, plot) ax.set(title=title) ax.grid() axes[-1].set_xlabel("x") plt.show()
  Error at iteration 25 is 0.15110768048491963. Error at iteration 50 is 0.05445871191314744. Error at iteration 75 is 0.019626740977814805. Error at iteration 100 is 0.007073413011029928. Error at iteration 125 is 0.0025492348261568765. Error at iteration 150 is 0.0009187358618483188. Error at iteration 175 is 0.00033110938827363157. Error at iteration 200 is 0.00011933073645664649. Converged in 205 iterations.
  The horizontal axis is the state $ x $, while the vertical axis gives $ s(x) $ and $ \phi(x) $.
  Overall, the policies match well with our predictions from above
  - Worker switches from one investment strategy to the other depending on relative return.
  - For low values of $ x $, the best option is to search for a new job.
  - Once $ x $ is larger, worker does better by investing in human capital specific to the current position., setting
  jv = JVWorker(grid_size=25, mc_size=50) plot_grid_max, plot_grid_size = 1.2, 100 plot_grid = np.linspace(0, plot_grid_max, plot_grid_size) fig, ax = plt.subplots() ax.set_xlim(0, plot_grid_max) ax.set_ylim(0, plot_grid_max)?
  jv = JVWorker(grid_size=25, mc_size=50) π, g, f_rvs, x_grid = jv.π, jv.g, jv.f_rvs, jv.x_grid T, get_greedy = operator_factory(jv) v_star = solve_model(jv, verbose=False) s_policy, ϕ_policy = get_greedy(v_star) # Turn the policy function arrays into actual functions s = lambda y: interp(x_grid, s_policy, y) ϕ = lambda y: interp(x_grid, ϕ_policy, y) def h(x, b, u): return (1 - b) * g(x, ϕ(x)) + b * max(g(x, ϕ(x)), u) plot_grid_max, plot_grid_size = 1.2, 100 plot_grid = np.linspace(0, plot_grid_max, plot_grid_size) fig, ax = plt.subplots(figsize=(8, 8)) ticks = (0.25, 0.5, 0.75, 1.0) ax.set(xticks=ticks, yticks=ticks, xlim=(0, plot_grid_max), ylim=(0, plot_grid_max), xlabel='$x_t$', ylabel='$x_{t+1}$') ax.plot(plot_grid, plot_grid, 'k--', alpha=0.6) # 45 degree line for x in plot_grid: for i in range(jv.mc_size): b = 1 if np.random.uniform(0, 1) < π(s(x)) else 0 u = f_rvs[i] y = h(x, b, u) ax.plot(x, y, 'go', alpha=0.25) plt.show()
  Looking at the dynamics, we can see that
  - If $ x_t $ is below about 0.2 the dynamics are random, but $ x_{t+1} > x_t $ is very likely.
  - As $ x_t $ increases the dynamics become deterministic, and $ x_t $ converges to a steady state value close to 1.
  Referring back to the figure here we see that $ x_t \approx 1 $ means that $ s_t = s(x_t) \approx 0 $ and $ \phi_t = \phi(x_t) \approx 0.6 $.
  jv = JVWorker() def xbar(ϕ): A, α = jv.A, jv.α return (A * ϕ**α)**(1 / (1 - α)) ϕ_grid = np.linspace(0, 1, 100) fig, ax = plt.subplots(figsize=(9, 7)) ax.set(xlabel='$\phi$') ax.plot(ϕ_grid, [xbar(ϕ) * (1 - ϕ) for ϕ in ϕ_grid], label='$w^*(\phi)$') ax.legend() plt.show()
  Observe that the maximizer is around 0.6.
  This is similar to the long-run value for $ \phi $ obtained in exercise 1.
  Hence the behavior of the infinitely patent worker is similar to that of the worker with $ \beta = 0.96 $.
  This seems reasonable and helps us confirm that our dynamic programming solutions are probably correct.