How to read this lecture...

Code should execute sequentially if run in a Jupyter notebook

Modeling Career Choice

Overview

Next we study a computational problem concerning career and job choices

The model is originally due to Derek Neal [Nea99]

This exposition draws on the presentation in [LS12], section 6.5

Model features

  • career and job within career both chosen to maximize expected discounted wage flow
  • infinite horizon dynamic programming with two states variables

Model

In what follows we distinguish between a career and a job, where

  • a career is understood to be a general field encompassing many possible jobs, and
  • a job is understood to be a position with a particular firm

For workers, wages can be decomposed into the contribution of job and career

  • \(w_t = \theta_t + \epsilon_t\), where

    • \(\theta_t\) is contribution of career at time \(t\)
    • \(\epsilon_t\) is contribution of job at time \(t\)

At the start of time \(t\), a worker has the following options

  • retain a current (career, job) pair \((\theta_t, \epsilon_t)\) — referred to hereafter as “stay put”
  • retain a current career \(\theta_t\) but redraw a job \(\epsilon_t\) — referred to hereafter as “new job”
  • redraw both a career \(\theta_t\) and a job \(\epsilon_t\) — referred to hereafter as “new life”

Draws of \(\theta\) and \(\epsilon\) are independent of each other and past values, with

  • \(\theta_t \sim F\)
  • \(\epsilon_t \sim G\)

Notice that the worker does not have the option to retain a job but redraw a career — starting a new career always requires starting a new job

A young worker aims to maximize the expected sum of discounted wages

(1)\[\mathbb{E} \sum_{t=0}^{\infty} \beta^t w_t\]

subject to the choice restrictions specified above

Let \(V(\theta, \epsilon)\) denote the value function, which is the maximum of (1) over all feasible (career, job) policies, given the initial state \((\theta, \epsilon)\)

The value function obeys

\[V(\theta, \epsilon) = \max\{I, II, III\},\]

where

(2)\[\begin{split}\begin{aligned} & I = \theta + \epsilon + \beta V(\theta, \epsilon) \\ & II = \theta + \int \epsilon' G(d \epsilon') + \beta \int V(\theta, \epsilon') G(d \epsilon') \nonumber \\ & III = \int \theta' F(d \theta') + \int \epsilon' G(d \epsilon') + \beta \int \int V(\theta', \epsilon') G(d \epsilon') F(d \theta') \nonumber \end{aligned}\end{split}\]

Evidently \(I\), \(II\) and \(III\) correspond to “stay put”, “new job” and “new life”, respectively

Parameterization

As in [LS12], section 6.5, we will focus on a discrete version of the model, parameterized as follows:

  • both \(\theta\) and \(\epsilon\) take values in the set linspace(0, B, N) — an even grid of \(N\) points between \(0\) and \(B\) inclusive
  • \(N = 50\)
  • \(B = 5\)
  • \(\beta = 0.95\)

The distributions \(F\) and \(G\) are discrete distributions generating draws from the grid points linspace(0, B, N)

A very useful family of discrete distributions is the Beta-binomial family, with probability mass function

\[p(k \,|\, n, a, b) = {n \choose k} \frac{B(k + a, n - k + b)}{B(a, b)}, \qquad k = 0, \ldots, n\]

Interpretation:

  • draw \(q\) from a Beta distribution with shape parameters \((a, b)\)
  • run \(n\) independent binary trials, each with success probability \(q\)
  • \(p(k \,|\, n, a, b)\) is the probability of \(k\) successes in these \(n\) trials

Nice properties:

  • very flexible class of distributions, including uniform, symmetric unimodal, etc.
  • only three parameters

Here’s a figure showing the effect of different shape parameters when \(n=50\)

../_images/beta-binom.png

The code that generated this figure can be found here

Implementation: career.jl

The code for solving the DP problem described above is found in this file, which is repeated here for convenience

#=
A type to solve the career / job choice model due to Derek Neal.

@author : Spencer Lyon <spencer.lyon@nyu.edu>

@date: 2014-08-05

References
----------

http://quant-econ.net/jl/career.html

[Neal1999] Neal, D. (1999). The Complexity of Job Mobility among Young Men,
Journal of Labor Economics, 17(2), 237-261.
=#

"""
Career/job choice model fo Derek Neal (1999)

##### Fields

- `beta::Real` : Discount factor in (0, 1)
- `N::Int` : Number of possible realizations of both epsilon and theta
- `B::Real` : upper bound for both epsilon and theta
- `theta::AbstractVector` : A grid of values on [0, B]
- `epsilon::AbstractVector` : A grid of values on [0, B]
- `F_probs::AbstractVector` : The pdf of each value associated with of F
- `G_probs::AbstractVector` : The pdf of each value associated with of G
- `F_mean::Real` : The mean of the distribution F
- `G_mean::Real` : The mean of the distribution G

"""
type CareerWorkerProblem
    beta::Real
    N::Int
    B::Real
    theta::AbstractVector
    epsilon::AbstractVector
    F_probs::AbstractVector
    G_probs::AbstractVector
    F_mean::Real
    G_mean::Real
end

"""
Constructor with default values for `CareerWorkerProblem`

##### Arguments

- `beta::Real(0.95)` : Discount factor in (0, 1)
- `B::Real(5.0)` : upper bound for both epsilon and theta
- `N::Real(50)` : Number of possible realizations of both epsilon and theta
- `F_a::Real(1), F_b::Real(1)` : Parameters of the distribution F
- `G_a::Real(1), G_b::Real(1)` : Parameters of the distribution F

##### Notes

There is also a version of this function that accepts keyword arguments for
each parameter
"""
# use key word argument 
function CareerWorkerProblem(;beta::Real=0.95, B::Real=5.0, N::Real=50,
                             F_a::Real=1, F_b::Real=1, G_a::Real=1,
                             G_b::Real=1)
    theta = linspace(0, B, N)
    epsilon = copy(theta)
    F_probs::Vector{Float64} = pdf(BetaBinomial(N-1, F_a, F_b))
    G_probs::Vector{Float64} = pdf(BetaBinomial(N-1, G_a, G_b))
    F_mean = sum(theta .* F_probs)
    G_mean = sum(epsilon .* G_probs)
    CareerWorkerProblem(beta, N, B, theta, epsilon, F_probs, G_probs,
                        F_mean, G_mean)
end


"""
Apply the Bellman operator for a given model and initial value.

##### Arguments

- `cp::CareerWorkerProblem` : Instance of `CareerWorkerProblem`
- `v::Matrix`: Current guess for the value function
- `out::Matrix` : Storage for output
- `;ret_policy::Bool(false)`: Toggles return of value or policy functions

##### Returns

None, `out` is updated in place. If `ret_policy == true` out is filled with the
policy function, otherwise the value function is stored in `out`.

"""
function update_bellman!(cp::CareerWorkerProblem, v::Array, out::Array;
                           ret_policy=false)
    # new life. This is a function of the distribution parameters and is
    # always constant. No need to recompute it in the loop
    v3 = (cp.G_mean + cp.F_mean + cp.beta .*
          cp.F_probs' * v * cp.G_probs)[1]  # don't need 1 element array

    for j=1:cp.N
        for i=1:cp.N
            # stay put
            v1 = cp.theta[i] + cp.epsilon[j] + cp.beta * v[i, j]

            # new job
            v2 = (cp.theta[i] .+ cp.G_mean .+ cp.beta .*
                  v[i, :]'*cp.G_probs)[1]  # don't need a single element array

            if ret_policy
                if v1 > max(v2, v3)
                    action = 1
                elseif v2 > max(v1, v3)
                    action = 2
                else
                    action = 3
                end
                out[i, j] = action
            else
                out[i, j] = max(v1, v2, v3)
            end
        end
    end
end


function update_bellman(cp::CareerWorkerProblem, v::Array; ret_policy=false)
    out = similar(v)
    update_bellman!(cp, v, out, ret_policy=ret_policy)
    return out
end

"""
Extract the greedy policy (policy function) of the model.

##### Arguments

- `cp::CareerWorkerProblem` : Instance of `CareerWorkerProblem`
- `v::Matrix`: Current guess for the value function
- `out::Matrix` : Storage for output

##### Returns

None, `out` is updated in place to hold the policy function

"""
function get_greedy!(cp::CareerWorkerProblem, v::Array, out::Array)
    update_bellman!(cp, v, out, ret_policy=true)
end


function get_greedy(cp::CareerWorkerProblem, v::Array)
    update_bellman(cp, v, ret_policy=true)
end

The code defines

  • a type CareerWorkerProblem that

    • encapsulates all the details of a particular parameterization
    • implement the Bellman operator \(T\)

In this model, \(T\) is defined by \(Tv(\theta, \epsilon) = \max\{I, II, III\}\), where \(I\), \(II\) and \(III\) are as given in (2), replacing \(V\) with \(v\)

The default probability distributions in CareerWorkerProblem correspond to discrete uniform distributions (see the Beta-binomial figure)

In fact all our default settings correspond to the version studied in [LS12], section 6.5.

Hence we can reproduce figures 6.5.1 and 6.5.2 shown there, which exhibit the value function and optimal policy respectively

Here’s the value function

../_images/career_vf.png

Value function with uniform probabilities

The code used to produce this plot was career/career_vf_plot.py

The optimal policy can be represented as follows (see Exercise 3 for code)

../_images/career_opt_pol.png

Interpretation:

  • If both job and career are poor or mediocre, the worker will experiment with new job and new career
  • If career is sufficiently good, the worker will hold it and experiment with new jobs until a sufficiently good one is found
  • If both job and career are good, the worker will stay put

Notice that the worker will always hold on to a sufficiently good career, but not necessarily hold on to even the best paying job

The reason is that high lifetime wages require both variables to be large, and the worker cannot change careers without changing jobs

  • Sometimes a good job must be sacrificed in order to change to a better career

Exercises

Exercise 1

Using the default parameterization in the type CareerWorkerProblem, generate and plot typical sample paths for \(\theta\) and \(\epsilon\) when the worker follows the optimal policy

In particular, modulo randomness, reproduce the following figure (where the horizontal axis represents time)

../_images/career_ex1_fig.png

Hint: To generate the draws from the distributions \(F\) and \(G\), use the type DiscreteRV

Exercise 2

Let’s now consider how long it takes for the worker to settle down to a permanent job, given a starting point of \((\theta, \epsilon) = (0, 0)\)

In other words, we want to study the distribution of the random variable

\[T^* := \text{the first point in time from which the worker's job no longer changes}\]

Evidently, the worker’s job becomes permanent if and only if \((\theta_t, \epsilon_t)\) enters the “stay put” region of \((\theta, \epsilon)\) space

Letting \(S\) denote this region, \(T^*\) can be expressed as the first passage time to \(S\) under the optimal policy:

\[T^* := \inf\{t \geq 0 \,|\, (\theta_t, \epsilon_t) \in S\}\]

Collect 25,000 draws of this random variable and compute the median (which should be about 7)

Repeat the exercise with \(\beta=0.99\) and interpret the change

Exercise 3

As best you can, reproduce the figure showing the optimal policy

Hint: The get_greedy() function returns a representation of the optimal policy where values 1, 2 and 3 correspond to “stay put”, “new job” and “new life” respectively. Use this and contourf from PyPlot.jl to produce the different shadings.

Now set G_a = G_b = 100 and generate a new figure with these parameters. Interpret.

Solutions

Exercise 1

srand(41)  # reproducible results
wp = CareerWorkerProblem()

function solve_wp(wp::CareerWorkerProblem)
    v_init = fill(100.0, wp.N, wp.N)
    func(x) = update_bellman(wp, x)
    v = compute_fixed_point(func, v_init, max_iter=500, verbose=false)
    optimal_policy = get_greedy(wp, v)
    return v, optimal_policy
end

v, optimal_policy = solve_wp(wp)

F = DiscreteRV(wp.F_probs)
G = DiscreteRV(wp.G_probs)

function gen_path(T=20)
    i = j = 1
    theta_ind = Int[]
    epsilon_ind = Int[]

    for t=1:T
        # do nothing if stay put
        if optimal_policy[i, j] == 2  # new job
            j = QuantEcon.draw(G)[1]
        elseif optimal_policy[i, j] == 3  # new life
            i, j = QuantEcon.draw(F)[1], QuantEcon.draw(G)[1]
        end
        push!(theta_ind, i)
        push!(epsilon_ind, j)
    end
    return wp.theta[theta_ind], wp.epsilon[epsilon_ind]
end

n = 2
thetas = []
epsilons = []
for i=1:n
    theta_path, epsilon_path = gen_path()
    push!(thetas, theta_path)
    push!(epsilons, epsilon_path)
end

plot(epsilons, label=["epsilon" ""], layout=grid(2,1), linewidth=2)
plot!(thetas, label=["theta" ""], linewidth=2, ylims=(0, 5))
../_images/career_solutions_ex1_jl.png

Exercise 2

The median for the original parameterization can be computed as follows

function gen_first_passage_time(optimal_policy::Matrix)
    t = 0
    i = j = 1
    while true
        if optimal_policy[i, j] == 1    # Stay put
            return t
        elseif optimal_policy[i, j] == 2  # New job
            j = QuantEcon.draw(G)[1]
        else                            # New life
            i, j = QuantEcon.draw(F)[1], QuantEcon.draw(G)[1]
        end
        t += 1
    end
end


M = 25000
samples = Array{Float64}(M)
for i=1:M
    samples[i] = gen_first_passage_time(optimal_policy)
end
print(median(samples))
7.0

To compute the median with \(\beta=0.99\) instead of the defalut value \(\beta=0.95\), replace wp=CareerWorkerProblem() with wp=CareerWorkerProblem(beta=0.99)

The medians are subject to randomness, but should be about 7 and 11 respectively. Not surprisingly, more patient workers will wait longer to settle down to their final job

wp2 = CareerWorkerProblem(beta=0.99)

v2, optimal_policy2 = solve_wp(wp2)

samples2 = Array{Float64}(M)
for i=1:M
    samples2[i] = gen_first_passage_time(optimal_policy2)
end
print(median(samples2))
14.0

Exercise 3

Here’s the code to reproduce the original figure

function region_plot(wp::CareerWorkerProblem)
    v, optimal_policy = solve_wp(wp)
    region_plot(optimal_policy)
end

function region_plot(optimal_policy::Matrix)
    lvls = [0.5, 1.5, 2.5, 3.5]

    contour(collect(wp.theta), collect(wp.epsilon),
            optimal_policy', fill=true, levels=lvls,
            alpha=0.5, cbar=false)
    xlabel!("theta")
    ylabel!("epsilon")
    annotate!([(1.8, 2.5,text("new life")),
              (4.5, 2.5, text("new job")),
              (4.5, 4.5, text("stay put", :white))])
end
region_plot(optimal_policy)
../_images/career_solutions_ex3_1_jl.png

Now we want to set G_a = G_b = 100 and generate a new figure with these parameters.

To do this replace:

wp = CareerWorkerProblem()

with:

wp = CareerWorkerProblem(G_a=100, G_b=100)

In the new figure, you will see that the region for which the worker will stay put has grown because the distribution for \(\epsilon\) has become more concentrated around the mean, making high-paying jobs less realistic

region_plot(CareerWorkerProblem(G_a=100, G_b=100))
../_images/career_solutions_ex3_2_jl.png