How to read this lecture...

Code should execute sequentially if run in a Jupyter notebook

A Lake Model of Employment and Unemployment


This lecture describes what has come to be called a lake model

The lake model is a basic tool for modeling unemployment

It allows us to analyze

  • flows between unemployment and employment
  • how these flows influence steady state employment and unemployment rates

It is a good model for interpreting monthly labor department reports on gross and net jobs created and jobs destroyed

The “lakes” in the model are the pools of employed and unemployed

The “flows” between the lakes are caused by

  • firing and hiring
  • entry and exit from the labor force

For the first part of this lecture, the parameters governing transitions into and out of unemployment and employment are exogenous

Later, we’ll determine some of these transition rates endogenously using the McCall search model

The only knowledge required for this lecture is

We’ll use some nifty concepts like ergodicity, which provides a fundamental link between cross sectional and long run time series distributions

These concepts will help us build an equilibrium model of ex ante homogeneous workers whose different luck generates variations in their ex post experiences

The Model

The economy is inhabited by a very large number of ex-ante identical workers

The workers live forever, spending their lives moving between unemployment and employment

Their rates of transition between employment and unemployment are governed by the following parameters:

  • \(\lambda\), the job finding rate for currently unemployed workers
  • \(\alpha\), the dismissal rate for currently employed workers
  • \(b\), the entry rate into the labor force
  • \(d\), the exit rate from the labor force

The growth rate of the labor force evidently equals \(g=b-d\)

Aggregate Variables

We want to derive the dynamics of the following aggregates

  • \(E_t\), the total number of employed workers at date \(t\)
  • \(U_t\), the total number of unemployed workers at \(t\)
  • \(N_t\), the number of workers in the labor force at \(t\)

We also want to know the values of the following objects

  • The employment rate \(e_t := E_t/N_t\)
  • The unemployment rate \(u_t := U_t/N_t\)

(Here and below, capital letters represent stocks and lowercase letters represent flows)

Laws of Motion for Stock Variables

We begin by constructing laws of motion for the aggregate variables \(E_t,U_t, N_t\)

Of the mass of workers \(E_t\) who are employed at date \(t\),

  • \((1-d)E_t\) will remain in the labor force
  • of these, \((1-\alpha)(1-d)E_t\) will remain employed

Of the mass of workers \(U_t\) workers who are currently unemployed,

  • \((1-d)U_t\) will remain in the labor force
  • of these, \((1-d) \lambda U_t\) will become employed

Therefore, the number of workers who will be employed at date \(t+1\) will be

\[E_{t+1} = (1-d)(1-\alpha)E_t + (1-d)\lambda U_t\]

A similar analysis implies

\[U_{t+1} = (1-d)\alpha E_t + (1-d)(1-\lambda)U_t + b (E_t+U_t)\]

The value \(b(E_t+U_t)\) is the mass of new workers entering the labor force unemployed

The total stock of workers \(N_t=E_t+U_t\) evolves as

\[N_{t+1} = (1+b-d)N_t = (1+g)N_t\]

Letting \(X_t := \left(\begin{matrix}E_t\\U_t\end{matrix}\right)\), the law of motion for \(X\) is

\[\begin{split}X_{t+1} = A X_t \quad \text{where} \quad A := \begin{pmatrix} (1-d)(1-\alpha) & (1-d)\lambda\\ (1-d)\alpha + b & (1-d)(1-\lambda) + b \end{pmatrix}\end{split}\]

This law tells us how total employment and unemployment evolve over time

Laws of Motion for Rates

Now let’s derive the law of motion for rates

To get these we can divide both sides of \(X_{t+1} = A X_t\) by \(N_{t+1}\) to get

\[\begin{split}\begin{pmatrix} E_{t+1}/N_{t+1} \\ U_{t+1}/N_{t+1} \end{pmatrix} = \frac1{1+g} A \begin{pmatrix} E_{t}/N_{t} \\ U_{t}/N_{t} \end{pmatrix}\end{split}\]


\[\begin{split}x_t := \left(\begin{matrix} e_t\\ u_t \end{matrix}\right) = \left(\begin{matrix} E_t/N_t\\ U_t/N_t \end{matrix}\right)\end{split}\]

we can also write this as

\[x_{t+1} = \hat A x_t \quad \text{where} \quad \hat A := \frac{1}{1 + g} A\]

You can check that \(e_t + u_t = 1\) implies that \(e_{t+1}+u_{t+1} = 1\)

This follows from the fact that the columns of \(\hat A\) sum to 1


Let’s code up these equations

Here’s the code:


@author : Victoria Gregory, John Stachurski

@date: 2016-08-03

Provides a class that simulates the dynamics of unemployment and employment in
the lake model.


type LakeModel
    lambda :: Float64
    alpha :: Float64
    b :: Float64
    d :: Float64
    g :: Float64
    A :: Matrix{Float64}
    A_hat :: Matrix{Float64}

Constructor with default values for `LakeModel`

##### Fields of `LakeModel`

  - lambda : job finding rate
  - alpha : dismissal rate
  - b : entry rate into labor force
  - d : exit rate from labor force
  - g : net entry rate
  - A : updates stock
  - A_hat : updates rate

function LakeModel(;lambda=0.283, alpha=0.013, b=0.0124, d=0.00822)

    g = b - d
    A = [ (1-d) * (1-alpha)  (1-d) * lambda;
          (1-d) * alpha + b (1-lambda) * (1-d) + b]
    A_hat = A ./ (1 + g)

    return LakeModel(lambda, alpha, b, d, g, A, A_hat)

Finds the steady state of the system :math:`x_{t+1} = \hat A x_{t}`

##### Arguments

  - lm : instance of `LakeModel`
  - tol: convergence tolerance

##### Returns

  - x : steady state vector of employment and unemployment rates


function rate_steady_state(lm::LakeModel, tol=1e-6)
    x = 0.5 * ones(2)
    error = tol + 1
    while (error > tol)
        new_x = lm.A_hat * x
        error = maxabs(new_x - x)
        x = new_x
    return x

Simulates the the sequence of Employment and Unemployent stocks

##### Arguments

  - X0 : contains initial values (E0, U0)
  - T : number of periods to simulate

##### Returns

  - X_path : contains sequence of employment and unemployment stocks


function simulate_stock_path(lm::LakeModel, X0::Vector{Float64}, T::Int)
    X_path = Array(Float64, 2, T)
    X = copy(X0)
    for t in 1:T
        X_path[:, t] = X
        X = lm.A * X
    return X_path

Simulates the the sequence of employment and unemployent rates.

##### Arguments

  - X0 : contains initial values (E0, U0)
  - T : number of periods to simulate

##### Returns

  - X_path : contains sequence of employment and unemployment rates

function simulate_rate_path(lm::LakeModel, x0::Vector{Float64}, T::Int)
    x_path = Array(Float64, 2, T)
    x = copy(x0)
    for t in 1:T
        x_path[:, t] = x
        x = lm.A_hat * x
    return x_path
lm = LakeModel()

2×2 Array{Float64,2}:
 0.978887   0.280674
 0.0252931  0.723506
lm = LakeModel(alpha = 2)

2×2 Array{Float64,2}:
 -0.99178  0.280674
  1.99596  0.723506

Aggregate Dynamics

Let’s run a simulation under the default parameters (see above) starting from \(X_0 = (138, 12)\)

using Plots

lm = LakeModel()
N_0 = 150      # Population
e_0 = 0.92     # Initial employment rate
u_0 = 1 - e_0  # Initial unemployment rate
T = 50         # Simulation length

E_0 = e_0 * N_0
U_0 = u_0 * N_0
X_0 = [E_0; U_0]

X_path = simulate_stock_path(lm, X_0, T)

titles = ["Employment" "Unemployment" "Labor Force"]
dates = collect(1:T)

x1 = squeeze(X_path[1, :], 1)
x2 = squeeze(X_path[2, :], 1)
x3 = squeeze(sum(X_path, 1), 1)

plot(dates, Vector[x1, x2, x3], title=titles, layout=(3, 1), legend=:none)

The aggregates \(E_t\) and \(U_t\) don’t converge because their sum \(E_t + U_t\) grows at rate \(g\)

On the other hand, the vector of employment and unemployment rates \(x_t\) can be in a steady state \(\bar x\) if there exists an \(\bar x\) such that

  • \(\bar x = \hat A \bar x\)
  • the components satisfy \(\bar e + \bar u = 1\)

This equation tells us that a steady state level \(\bar x\) is an eigenvector of \(\hat A\) associated with a unit eigenvalue

We also have \(x_t \to \bar x\) as \(t \to \infty\) provided that the remaining eigenvalue of \(\hat A\) has modulus less that 1

This is the case for our default parameters:

lm = LakeModel()

e,f = eigvals(lm.A_hat)

abs(e), abs(f)

Let’s look at the convergence of the unemployment and employment rate to steady state levels (dashed red line)

lm = LakeModel()
e_0 = 0.92     # Initial employment rate
u_0 = 1 - e_0  # Initial unemployment rate
T = 50         # Simulation length

xbar = rate_steady_state(lm)
x_0 = [e_0; u_0]
x_path = simulate_rate_path(lm, x_0, T)

titles = ["Employment rate" "Unemployment rate"]
dates = collect(1:T)

plot(dates, x_path', layout=(2, 1), title=titles, legend=:none)
hline!(xbar', layout=(2, 1), color=:red, linestyle=:dash,
       ylims=[0.9999*minimum(x_path[1,:]) Inf;
       -Inf 1.001*maximum(x_path[2,:])])

Dynamics of an Individual Worker

A individual worker’s employment dynamics are governed by a finite state Markov process

The worker can be in one of two states:

  • \(s_t=0\) means unemployed
  • \(s_t=1\) means employed

Let’s start off under the assumption that \(b = d = 0\)

The associated transition matrix is then

\[\begin{split}P = \left( \begin{matrix} 1 - \lambda & \lambda \\ \alpha & 1 - \alpha \end{matrix} \right)\end{split}\]

Let \(\psi_t\) denote the marginal distribution over employment / unemployment states for the worker at time \(t\)

As usual, we regard it as a row vector

We know from an earlier discussion that \(\psi_t\) follows the law of motion

\[\psi_{t+1} = \psi_t P\]

We also know from the lecture on finite Markov chains that if \(\alpha \in (0, 1)\) and \(\lambda \in (0, 1)\), then \(P\) is stationary and ergodic, with a unique stationary distribution \(\psi^*\)

The unique stationary distribution is defined by

\[\psi^*[0] = \frac{\alpha}{\alpha + \lambda}\]

Not surprisingly, probability mass on the unemployment state increases with the dismissal rate and falls with the job finding rate rate


Let’s look at a typical lifetime of employment-unemployment spells

We want to compute the average amounts of time an infinitely lived worker would spend employed and unemployed


\[\bar s_{u,T} := \frac1{T} \sum_{t=1}^T \mathbb 1\{s_t = 0\}\]


\[\bar s_{e,T} := \frac1{T} \sum_{t=1}^T \mathbb 1\{s_t = 1\}\]

(As usual, \(\mathbb 1\{Q\} = 1\) if statement \(Q\) is true and 0 otherwise)

These are the fraction of time a worker spends unemployed and employed, respectively, up until period \(T\)

As stated above, if \(\alpha \in (0, 1)\) and \(\lambda \in (0, 1)\), then \(P\) is ergodic, and hence we have

\[\lim_{T \to \infty} \bar s_{u, T} = \psi^*[0] \quad \text{and} \quad \lim_{T \to \infty} \bar s_{e, T} = \psi^*[1]\]

with probability one

Inspection tells us that \(P\) is exactly the transpose of \(\hat A\) under the assumption \(b=d=0\)

Thus, the percentages of time that an infinitely lived worker spends employed and unemployed equal the fractions of workers employed and unemployed in the steady state distribution

Convergence rate

How long does it take for time series sample averages to converge to cross sectional averages?

We can use QuantEcon.jl‘s MarkovChain class to investigate this

Let’s plot the path of the sample averages over 5,000 periods

Agent dynamics in the lake model

using QuantEcon

lm = LakeModel(d=0, b=0)
T = 5000     # Simulation length

alpha, lambda = lm.alpha, lm.lambda
P = [(1 - lambda) lambda; alpha (1 - alpha)]

mc = MarkovChain(P, [0; 1])     # 0=unemployed, 1=employed
xbar = rate_steady_state(lm)

s_path = simulate(mc, T; init=2)
s_bar_e = cumsum(s_path) ./ (1:T)
s_bar_u = 1 - s_bar_e
s_bars = [s_bar_u s_bar_e]

titles = ["Proportion of time unemployed" "Proportion of time employed"]
dates = collect(1:T)

plot(dates, s_bars, layout=(2, 1), title=titles, legend=:none)
hline!(reverse(xbar)', layout=(2, 1), color=:red, linestyle=:dash)

The stationary probabilities are given by the dashed red line

In this case it takes much of the sample for these two objects to converge

This is largely due to the high persistence in the Markov chain

Endogenous Job Finding Rate

We now make the hiring rate endogenous

The transition rate from unemployment to employment will be determined by the McCall search model [McC70]

All details relevant to the following discussion can be found in our treatment of that model

Reservation Wage

The most important thing to remember about the model is that optimal decisions are characterized by a reservation wage \(\bar w\)

  • if the wage offer \(w\) in hand is greater than or equal to \(\bar w\), then the worker accepts
  • otherwise, the worker rejects

As we saw in our discussion of the model, the reservation wage depends on the wage offer distribution and the parameters

  • \(\alpha\) (the separation rate)
  • \(\beta\) (the discount factor)
  • \(\gamma\) (the offer arrival rate)
  • \(c\) (unemployment compensation)

Linking the McCall Search Model to the Lake Model

Suppose that all workers inside a Lake Model behave according to the McCall search model

The exogenous probability of leaving employment remains \(\alpha\)

But their optimal decision rules determine the probability \(\lambda\) of leaving unemployment

This is now

(1)\[\lambda = \gamma \mathbb P \{ w_t \geq \bar w\} = \gamma \sum_{w' \geq \bar w} p(w')\]

Fiscal Policy

We can use the McCall search version of the Lake Model to find an optimal level of unemployment insurance

We assume that the government sets unemployment compensation \(c\)

The government imposes a lump sum tax \(\tau\) sufficient to finance total unemployment payments

To attain a balanced budget at a steady state, taxes, the steady state unemployment rate \(u\), and the unemployment compensation rate must satisfy

\[\tau = u c\]

The lump sum tax applies to everyone, including unemployed workers

Thus, the post-tax income of an employed worker with wage \(w\) is \(w - \tau\)

The post-tax income of an unemployed worker is \(c - \tau\).

For each specification \((c, \tau)\) of government policy, we can solve for the worker’s optimal reservation wage

This determines \(\lambda\) via (1) evaluated at post tax wages, which in turn determines a steady state unemployment rate \(u(c, \tau)\)

For a given level of unemployment benefit \(c\), we can solve for a tax that balances the budget in the steady state

\[\tau = u(c, \tau) c\]

To evaluate alternative government tax-unemployment compensation pairs, we require a welfare criterion

We use a steady state welfare criterion

\[W := e \, {\mathbb E} [V \, | \, \text{employed}] + u \, U\]

where the notation \(V\) and \(U\) is as defined in the McCall search model lecture

The wage offer distribution will be a discretized version of the distribution \(\mathcal{N}(\log(20),1)\), as shown in the next figure


We take a period to be a month

We set \(b\) and \(d\) to match monthly birth and death rates, respectively, in the U.S. population

  • \(b = 0.0124\)
  • \(d = 0.00822\)

Following [DFH06], we set \(\alpha\), the hazard rate of leaving employment, to

  • \(\alpha = 0.013\)

Fiscal Policy Code

We will make use of code we wrote in the McCall model lecture, embedded below for convenience

The first piece of code implements value function iteration


Implements iteration on the Bellman equations to solve the McCall growth model


using Distributions

# A default utility function

function u(c, sigma)
    if c > 0
        return (c^(1 - sigma) - 1) / (1 - sigma)
        return -10e6

# default wage vector with probabilities

const n = 60                                   # n possible outcomes for wage
const default_w_vec = linspace(10, 20, n)   # wages between 10 and 20
const a, b = 600, 400                          # shape parameters
const dist = BetaBinomial(n-1, a, b)
const default_p_vec = pdf(dist)

type McCallModel
    alpha::Float64        # Job separation rate
    beta::Float64         # Discount rate
    gamma::Float64        # Job offer rate
    c::Float64            # Unemployment compensation
    sigma::Float64        # Utility parameter
    w_vec::Vector{Float64} # Possible wage values
    p_vec::Vector{Float64} # Probabilities over w_vec

    function McCallModel(alpha=0.2,

        return new(alpha, beta, gamma, c, sigma, w_vec, p_vec)

A function to update the Bellman equations.  Note that V_new is modified in
place (i.e, modified by this function).  The new value of U is returned.

function update_bellman!(mcm, V, V_new, U)
    # Simplify notation
    alpha, beta, sigma, c, gamma = mcm.alpha, mcm.beta, mcm.sigma, mcm.c, mcm.gamma

    for (w_idx, w) in enumerate(mcm.w_vec)
        # w_idx indexes the vector of possible wages
        V_new[w_idx] = u(w, sigma) + beta * ((1 - alpha) * V[w_idx] + alpha * U)

    U_new = u(c, sigma) + beta * (1 - gamma) * U +
                    beta * gamma * sum(max(U, V) .* mcm.p_vec)

    return U_new

function solve_mccall_model(mcm; tol::Float64=1e-5, max_iter::Int=2000)

    V = ones(length(mcm.w_vec))  # Initial guess of V
    V_new = similar(V)           # To store updates to V
    U = 1.0                        # Initial guess of U
    i = 0
    error = tol + 1

    while error > tol && i < max_iter
        U_new = update_bellman!(mcm, V, V_new, U)
        error_1 = maximum(abs(V_new - V))
        error_2 = abs(U_new - U)
        error = max(error_1, error_2)
        V[:] = V_new
        U = U_new
        i += 1

    return V, U

The second piece of code repeated from the McCall model lecture is used to complete the reservation wage

Compute the reservation wage in the McCall model.


Computes the reservation wage of an instance of the McCall model
by finding the smallest w such that V(w) > U.

If V(w) > U for all w, then the reservation wage w_bar is set to
the lowest wage in mcm.w_vec.

If v(w) < U for all w, then w_bar is set to np.inf.

mcm : an instance of McCallModel
return_values : bool (optional, default=false)
    Return the value functions as well 

w_bar : scalar
    The reservation wage
function compute_reservation_wage(mcm; return_values=false)

    V, U = solve_mccall_model(mcm)
    w_idx = searchsortedfirst(V - U, 0)  

    if w_idx == length(V)
        w_bar = Inf
        w_bar = mcm.w_vec[w_idx]

    if return_values == false
        return w_bar
        return w_bar, V, U


Now let’s compute and plot welfare, employment, unemployment, and tax revenue as a function of the unemployment compensation rate

using Distributions
using Roots

# Some global variables that will stay constant
alpha = 0.013
alpha_q = (1-(1-alpha)^3)
b = 0.0124
d = 0.00822
beta = 0.98
gamma = 1.0
sigma = 2.0

# The default wage distribution: a discretized log normal
log_wage_mean, wage_grid_size, max_wage = 20, 200, 170
w_vec = linspace(1e-3, max_wage, wage_grid_size + 1)
logw_dist = Normal(log(log_wage_mean), 1)
cdf_logw = cdf(logw_dist, log(w_vec))
pdf_logw = cdf_logw[2:end] - cdf_logw[1:end-1]
p_vec = pdf_logw ./ sum(pdf_logw)
w_vec = (w_vec[1:end-1] + w_vec[2:end]) / 2

Compute the reservation wage, job finding rate and value functions of the
workers given c and tau.

function compute_optimal_quantities(c::Float64, tau::Float64)
  mcm = McCallModel(alpha_q,
                    c-tau,           # post-tax compensation
                    collect(w_vec-tau),   # post-tax wages

  w_bar, V, U = compute_reservation_wage(mcm, return_values=true)
  lmda = gamma * sum(p_vec[w_vec-tau .> w_bar])

  return w_bar, lmda, V, U

Compute the steady state unemployment rate given c and tau using optimal
quantities from the McCall model and computing corresponding steady state

function compute_steady_state_quantities(c::Float64, tau::Float64)
  w_bar, lmda, V, U = compute_optimal_quantities(c, tau)

  # Compute steady state employment and unemployment rates
  lm = LakeModel(lambda=lmda, alpha=alpha_q, b=b, d=d)
  x = rate_steady_state(lm)
  e_rate, u_rate = x

  # Compute steady state welfare
  w = sum(V .* p_vec .* (w_vec - tau .> w_bar)) / sum(p_vec .* (w_vec - tau .> w_bar))
  welfare = e_rate .* w + u_rate .* U

  return e_rate, u_rate, welfare

Find tax level that will induce a balanced budget.

function find_balanced_budget_tax(c::Float64)
  function steady_state_budget(t::Float64)
    e_rate, u_rate, w = compute_steady_state_quantities(c, t)
    return t - u_rate * c

  tau = fzero(steady_state_budget, 0.0, 0.9 * c)

  return tau

# Levels of unemployment insurance we wish to study
Nc = 60
c_vec = linspace(5.0, 140.0, Nc)

tax_vec = Array(Float64, Nc,)
unempl_vec = Array(Float64, Nc,)
empl_vec = Array(Float64, Nc,)
welfare_vec = Array(Float64, Nc,)

for i = 1:Nc
  t = find_balanced_budget_tax(c_vec[i])
  e_rate, u_rate, welfare = compute_steady_state_quantities(c_vec[i], t)
  tax_vec[i] = t
  unempl_vec[i] = u_rate
  empl_vec[i] = e_rate
  welfare_vec[i] = welfare

titles = ["Unemployment" "Employment" "Tax" "Welfare"]
plot(c_vec, [unempl_vec empl_vec tax_vec welfare_vec],
    color=:blue, lw=2, alpha=0.7, title=titles, legend=:none, layout=(2,2))

The figure that the preceding code listing generates is shown below


Welfare first increases and then decreases as unemployment benefits rise

The level that maximizes steady state welfare is approximately 62


Exercise 1

Consider an economy with initial stock of workers \(N_0 = 100\) at the steady state level of employment in the baseline parameterization

  • \(\alpha = 0.013\)
  • \(\lambda = 0.283\)
  • \(b = 0.0124\)
  • \(d = 0.00822\)

(The values for \(\alpha\) and \(\lambda\) follow [DFH06])

Suppose that in response to new legislation the hiring rate reduces to \(\lambda = 0.2\)

Plot the transition dynamics of the unemployment and employment stocks for 50 periods

Plot the transition dynamics for the rates

How long does the economy take to converge to its new steady state?

What is the new steady state level of employment?

Exercise 2

Consider an economy with initial stock of workers \(N_0 = 100\) at the steady state level of employment in the baseline parameterization

Suppose that for 20 periods the birth rate was temporarily high (\(b = 0.0025\)) and then returned to its original level

Plot the transition dynamics of the unemployment and employment stocks for 50 periods

Plot the transition dynamics for the rates

How long does the economy take to return to its original steady state?


Exercise 1

Construct the class containing the default parameters and record the steady state.

lm = LakeModel()
x0 = rate_steady_state(lm)
println("Initial Steady State: $x0")
Initial Steady State: [0.917332,0.0826681]

Initialize the simulation values

N0 = 100
T = 50

New legislation changes \(\lambda\) to \(0.2\)

lm =LakeModel(lambda = 0.2)
LakeModel(0.2,0.013,0.0124,0.00822,0.00418,[0.978887 0.198356; 0.0252931 0.805824],[0.974812 0.19753; 0.0251879 0.80247])
xbar = rate_steady_state(lm) # new steady state
X_path = simulate_stock_path(lm,x0 * N0, T)
x_path = simulate_rate_path(lm,x0, T)
println("New Steady State: $xbar")
New Steady State: [0.886904,0.113096]

Now plot stocks

titles = ["Employment" "Unemployment" "Labor Force"]
dates = collect(1:T)

x1 = X_path[1, :]
x2 = X_path[2, :]
x3 = squeeze(sum(X_path, 1), 1)

plot(dates, Vector[x1, x2, x3], title=titles, layout=(3, 1), legend=:none)

And how the rates evolve:

titles = ["employment rate" "unemployment rate"]
dates = collect(1:T)

x1 = x_path[1, :]
x2 = x_path[2, :]

p=plot(dates, Vector[x1, x2], title=titles, layout=(2, 1), legend=:none)
hline!((xbar)', layout=(2, 1), color=:red, linestyle=:dash)

We see that it takes 20 periods for the economy to converge to it’s new steady state levels

Exercise 2

This next exercise has the economy expriencing a boom in entrances to the labor market and then later returning to the original levels. For 20 periods the economy has a new entry rate into the labor market

Let’s start off at the baseline parameterization and record the steady state:

lm = LakeModel()
x0 = rate_steady_state(lm)
2-element Array{Float64,1}:

Here are the other parameters:

b_hat = 0.003
T_hat = 20

Let’s increase \(b\) to the new value and simulate for 20 periods

lm = LakeModel(b=b_hat)
X_path1 = simulate_stock_path(lm,x0 * N0, T_hat) # simulate stocks
x_path1 = simulate_rate_path(lm,x0, T_hat) # simulate rates

Now we reset \(b\) to the original value and then, using the state after 20 periods for the new initial conditions, we simulate for the additional 30 periods

lm = LakeModel(b=0.0124)
X_path2 = simulate_stock_path(lm,X_path1[:,end-1],T-T_hat+1) # simulate stocks
x_path2 = simulate_rate_path(lm,x_path1[:,end-1],T-T_hat+1) # simulate rates

Finally we combine these two paths and plot

x_path = hcat(x_path1,x_path2[:,2:end]) # note [1:] to avoid doubling period 20
X_path = hcat(X_path1,X_path2[:,2:end])
titles = ["employment" "unemployment" "labor force"]
dates = collect(1:T)

x1 = X_path[1,:]
x2 = X_path[2,:]
x3 = squeeze(sum(X_path, 1), 1)

plot(dates, Vector[x1, x2, x3], title=titles, layout=(3, 1), legend=:none)

And the rates:

titles = ["employment Rate" "unemployment Rate"]
dates = collect(1:T)

x1 = x_path[1,:]
x2 = x_path[2,:]

p=plot(dates, Vector[x1, x2], title=titles, layout=(2, 1), legend=:none)
hline!((x0)', layout=(2, 1), color=:red, linestyle=:dash)