34. Cass-Koopmans Model#

34.1. Overview#

This lecture and Cass-Koopmans Competitive Equilibrium describe a model that Tjalling Koopmans [Koopmans, 1965] and David Cass [Cass, 1965] used to analyze optimal growth.

The model extends the model of Robert Solow described in an earlier lecture.

It does so by making saving rate be a decision, instead of a hard-wired constant.

(Solow assumed a constant saving rate determined outside the model.)

We describe two versions of the model, a planning problem in this lecture, and a competitive equilibrium in this lecture Cass-Koopmans Competitive Equilibrium.

Together, the two lectures illustrate what is, in fact, a more general connection between a planned economy and a decentralized economy organized as a competitive equilibrium.

This lecture is devoted to the planned economy version.

In the planned economy, there are

  • no prices

  • no budget constraints

Instead there is a dictator that tells people

  • what to produce

  • what to invest in physical capital

  • who is to consume what and when

The lecture uses important ideas including

  • A min-max problem for solving a planning problem.

  • A shooting algorithm for solving difference equations subject to initial and terminal conditions.

  • A turnpike property of optimal paths for long but finite-horizon economies.

  • A stable manifold and a phase plane

In addition to what’s in Anaconda, this lecture will need the following libraries:

!pip install quantecon
Hide code cell output
Collecting quantecon
  Downloading quantecon-0.8.1-py3-none-any.whl.metadata (5.2 kB)
Requirement already satisfied: numba>=0.49.0 in /home/runner/miniconda3/envs/quantecon/lib/python3.12/site-packages (from quantecon) (0.60.0)
Requirement already satisfied: numpy>=1.17.0 in /home/runner/miniconda3/envs/quantecon/lib/python3.12/site-packages (from quantecon) (1.26.4)
Requirement already satisfied: requests in /home/runner/miniconda3/envs/quantecon/lib/python3.12/site-packages (from quantecon) (2.32.3)
Requirement already satisfied: scipy>=1.5.0 in /home/runner/miniconda3/envs/quantecon/lib/python3.12/site-packages (from quantecon) (1.13.1)
Requirement already satisfied: sympy in /home/runner/miniconda3/envs/quantecon/lib/python3.12/site-packages (from quantecon) (1.14.0)
Requirement already satisfied: llvmlite<0.44,>=0.43.0dev0 in /home/runner/miniconda3/envs/quantecon/lib/python3.12/site-packages (from numba>=0.49.0->quantecon) (0.43.0)
Requirement already satisfied: charset-normalizer<4,>=2 in /home/runner/miniconda3/envs/quantecon/lib/python3.12/site-packages (from requests->quantecon) (3.3.2)
Requirement already satisfied: idna<4,>=2.5 in /home/runner/miniconda3/envs/quantecon/lib/python3.12/site-packages (from requests->quantecon) (3.7)
Requirement already satisfied: urllib3<3,>=1.21.1 in /home/runner/miniconda3/envs/quantecon/lib/python3.12/site-packages (from requests->quantecon) (2.2.3)
Requirement already satisfied: certifi>=2017.4.17 in /home/runner/miniconda3/envs/quantecon/lib/python3.12/site-packages (from requests->quantecon) (2024.8.30)
Requirement already satisfied: mpmath<1.4,>=1.1.0 in /home/runner/miniconda3/envs/quantecon/lib/python3.12/site-packages (from sympy->quantecon) (1.3.0)
Downloading quantecon-0.8.1-py3-none-any.whl (322 kB)
Installing collected packages: quantecon
Successfully installed quantecon-0.8.1

Let’s start with some standard imports:

import matplotlib.pyplot as plt
from numba import jit, float64
from numba.experimental import jitclass
import numpy as np
from quantecon.optimize import brentq

34.2. The Model#

Time is discrete and takes values t=0,1,,T where T is finite.

(We’ll eventually study a limiting case in which T=+)

A single good can either be consumed or invested in physical capital.

The consumption good is not durable and depreciates completely if not consumed immediately.

The capital good is durable but depreciates.

We let Ct be the total consumption of a nondurable consumption good at time t.

Let Kt be the stock of physical capital at time t.

Let C = {C0,,CT} and K = {K0,,KT+1}.

34.2.1. Digression: Aggregation Theory#

We use a concept of a representative consumer to be thought of as follows.

There is a unit mass of identical consumers indexed by ω[0,1].

Consumption of consumer ω is c(ω).

Aggregate consumption is

C=01c(ω)dω

Consider a welfare problem that chooses an allocation {c(ω)} across consumers to maximize

01u(c(ω))dω

where u() is a concave utility function with u>0,u<0 and maximization is subject to

(34.1)#C=01c(ω)dω.

Form a Lagrangian L=01u(c(ω))dω+λ[C01c(ω)dω].

Differentiate under the integral signs with respect to each ω to obtain the first-order necessary conditions

u(c(ω))=λ.

These conditions imply that c(ω) equals a constant c that is independent of ω.

To find c, use feasibility constraint (34.1) to conclude that

c(ω)=c=C.

This line of argument indicates the special aggregation theory that lies beneath outcomes in which a representative consumer consumes amount C.

It appears often in aggregate economics.

We shall use this aggregation theory here and also in this lecture Cass-Koopmans Competitive Equilibrium.

34.2.1.1. An Economy#

A representative household is endowed with one unit of labor at each t and likes the consumption good at each t.

The representative household inelastically supplies a single unit of labor Nt at each t, so that Nt=1 for all t{0,1,,T}.

The representative household has preferences over consumption bundles ordered by the utility functional:

(34.2)#U(C)=t=0TβtCt1γ1γ

where β(0,1) is a discount factor and γ>0 governs the curvature of the one-period utility function.

Larger γ’s imply more curvature.

Note that

(34.3)#u(Ct)=Ct1γ1γ

satisfies u>0,u<0.

u>0 asserts that the consumer prefers more to less.

u<0 asserts that marginal utility declines with increases in Ct.

We assume that K0>0 is an exogenous initial capital stock.

There is an economy-wide production function

(34.4)#F(Kt,Nt)=AKtαNt1α

with 0<α<1, A>0.

A feasible allocation C,K satisfies

(34.5)#Ct+Kt+1F(Kt,Nt)+(1δ)Ktfor all t{0,1,,T}

where δ(0,1) is a depreciation rate of capital.

34.3. Planning Problem#

A planner chooses an allocation {C,K} to maximize (34.2) subject to (34.5).

Let μ={μ0,,μT} be a sequence of nonnegative Lagrange multipliers.

To find an optimal allocation, form a Lagrangian

(34.6)#L(C,K,μ)=t=0Tβt{u(Ct)+μt(F(Kt,1)+(1δ)KtCtKt+1)}

and pose the following min-max problem:

(34.7)#minμmaxC,KL(C,K,μ)
  • Extremization means maximization with respect to C,K and minimization with respect to μ.

  • Our problem satisfies conditions that assure that second-order conditions are satisfied at an allocation that satisfies the first-order necessary conditions that we are about to compute.

Before computing first-order conditions, we present some handy formulas.

34.3.1. Useful Properties of Linearly Homogeneous Production Function#

The following technicalities will help us.

Notice that

F(Kt,Nt)=AKtαNt1α=NtA(KtNt)α

Define the output per-capita production function

F(Kt,Nt)Ntf(KtNt)=A(KtNt)α

whose argument is capital per-capita.

It is useful to recall the following calculations for the marginal product of capital

(34.8)#F(Kt,Nt)Kt=Ntf(KtNt)Kt=Ntf(KtNt)1Nt(Chain rule)=f(KtNt)|Nt=1=f(Kt)

and the marginal product of labor

F(Kt,Nt)Nt=Ntf(KtNt)Nt(Product rule)=f(KtNt)+Ntf(KtNt)KtNt2(Chain rule)=f(KtNt)KtNtf(KtNt)|Nt=1=f(Kt)f(Kt)Kt

(Here we are using that Nt=1 for all t, so that Kt=KtNt.)

34.3.2. First-order necessary conditions#

We now compute first-order necessary conditions for extremization of Lagrangian (34.6):

(34.9)#Ct:u(Ct)μt=0for allt=0,1,,T
(34.10)#Kt:βμt[(1δ)+f(Kt)]μt1=0for all t=1,2,,T
(34.11)#μt:F(Kt,1)+(1δ)KtCtKt+1=0for all t=0,1,,T
(34.12)#KT+1:μT0, 0 if KT+1=0; =0 if KT+1>0

In computing (34.10) we recognize that Kt appears in both the time t and time t1 feasibility constraints (34.5).

Restrictions (34.12) come from differentiating with respect to KT+1 and applying the following Karush-Kuhn-Tucker condition (KKT) (see Karush-Kuhn-Tucker conditions):

(34.13)#μTKT+1=0

Combining (34.9) and (34.10) gives

βu(Ct)[(1δ)+f(Kt)]u(Ct1)=0 for all t=1,2,,T+1

which can be rearranged to become

(34.14)#βu(Ct+1)[(1δ)+f(Kt+1)]=u(Ct) for all t=0,1,,T

Applying the inverse marginal utility of consumption function on both sides of the above equation gives

Ct+1=u1((βu(Ct)[f(Kt+1)+(1δ)])1)

which for our utility function (34.3) becomes the consumption Euler equation

(34.15)#Ct+1=(βCtγ[f(Kt+1)+(1δ)])1/γ

which we can combine with the feasibility constraint (34.5) to get

(34.16)#Ct+1=Ct(β[f(F(Kt,1)+(1δ)KtCt)+(1δ)])1/γKt+1=F(Kt,1)+(1δ)KtCt.

This is a pair of non-linear first-order difference equations that map Ct,Kt into Ct+1,Kt+1 and that an optimal sequence C,K must satisfy.

It must also satisfy the initial condition that K0 is given and KT+1=0.

Below we define a jitclass that stores parameters and functions that define our economy.

planning_data = [
    ('γ', float64),    # Coefficient of relative risk aversion
    ('β', float64),    # Discount factor
    ('δ', float64),    # Depreciation rate on capital
    ('α', float64),    # Return to capital per capita
    ('A', float64)     # Technology
]
@jitclass(planning_data)
class PlanningProblem():

    def __init__(self, γ=2, β=0.95, δ=0.02, α=0.33, A=1):

        self.γ, self.β = γ, β
        self.δ, self.α, self.A = δ, α, A

    def u(self, c):
        '''
        Utility function
        ASIDE: If you have a utility function that is hard to solve by hand
        you can use automatic or symbolic differentiation
        See https://github.com/HIPS/autograd
        '''
        γ = self.γ

        return c ** (1 - γ) / (1 - γ) if γ!= 1 else np.log(c)

    def u_prime(self, c):
        'Derivative of utility'
        γ = self.γ

        return c ** (-γ)

    def u_prime_inv(self, c):
        'Inverse of derivative of utility'
        γ = self.γ

        return c ** (-1 / γ)

    def f(self, k):
        'Production function'
        α, A = self.α, self.A

        return A * k ** α

    def f_prime(self, k):
        'Derivative of production function'
        α, A = self.α, self.A

        return α * A * k ** (α - 1)

    def f_prime_inv(self, k):
        'Inverse of derivative of production function'
        α, A = self.α, self.A

        return (k / (A * α)) ** (1 / (α - 1))

    def next_k_c(self, k, c):
        ''''
        Given the current capital Kt and an arbitrary feasible
        consumption choice Ct, computes Kt+1 by state transition law
        and optimal Ct+1 by Euler equation.
        '''
        β, δ = self.β, self.δ
        u_prime, u_prime_inv = self.u_prime, self.u_prime_inv
        f, f_prime = self.f, self.f_prime

        k_next = f(k) + (1 - δ) * k - c
        c_next = u_prime_inv(u_prime(c) / (β * (f_prime(k_next) + (1 - δ))))

        return k_next, c_next

We can construct an economy with the Python code:

pp = PlanningProblem()

34.4. Shooting Algorithm#

We use shooting to compute an optimal allocation C,K and an associated Lagrange multiplier sequence μ.

First-order necessary conditions (34.9), (34.10), and (34.11) for the planning problem form a system of difference equations with two boundary conditions:

  • K0 is a given initial condition for capital

  • KT+1=0 is a terminal condition for capital that we deduced from the first-order necessary condition for KT+1 the KKT condition (34.13)

We have no initial condition for the Lagrange multiplier μ0.

If we did, our job would be easy:

  • Given μ0 and k0, we could compute c0 from equation (34.9) and then k1 from equation (34.11) and μ1 from equation (34.10).

  • We could continue in this way to compute the remaining elements of C,K,μ.

However, we woujld not be assured that the Kuhn-Tucker condition (34.13) would be satisfied.

Furthermore, we don’t have an initial condition for μ0.

So this won’t work.

Indeed, part of our task is to compute the optimal value of μ0.

To compute μ0 and the other objects we want, a simple modification of the above procedure will work.

It is called the shooting algorithm.

It is an instance of a guess and verify algorithm that consists of the following steps:

  • Guess an initial Lagrange multiplier μ0.

  • Apply the simple algorithm described above.

  • Compute KT+1 and check whether it equals zero.

  • If KT+1=0, we have solved the problem.

  • If KT+1>0, lower μ0 and try again.

  • If KT+1<0, raise μ0 and try again.

The following Python code implements the shooting algorithm for the planning problem.

(Actually, we modified the preceding algorithm slightly by starting with a guess for c0 instead of μ0 in the following code.)

@jit
def shooting(pp, c0, k0, T=10):
    '''
    Given the initial condition of capital k0 and an initial guess
    of consumption c0, computes the whole paths of c and k
    using the state transition law and Euler equation for T periods.
    '''
    if c0 > pp.f(k0) + (1 - pp.δ) * k0:
        print("initial consumption is not feasible")

        return None

    # initialize vectors of c and k
    c_vec = np.empty(T+1)
    k_vec = np.empty(T+2)

    c_vec[0] = c0
    k_vec[0] = k0

    for t in range(T):
        k_vec[t+1], c_vec[t+1] = pp.next_k_c(k_vec[t], c_vec[t])

    k_vec[T+1] = pp.f(k_vec[T]) + (1 - pp.δ) * k_vec[T] - c_vec[T]

    return c_vec, k_vec

We’ll start with an incorrect guess.

paths = shooting(pp, 0.2, 0.3, T=10)
fig, axs = plt.subplots(1, 2, figsize=(14, 5))

colors = ['blue', 'red']
titles = ['Consumption', 'Capital']
ylabels = ['$c_t$', '$k_t$']

T = paths[0].size - 1
for i in range(2):
    axs[i].plot(paths[i], c=colors[i])
    axs[i].set(xlabel='t', ylabel=ylabels[i], title=titles[i])

axs[1].scatter(T+1, 0, s=80)
axs[1].axvline(T+1, color='k', ls='--', lw=1)

plt.show()
_images/bbb28ec13b2ea7ea268ee53281aceab9352f04dac6a05f6878f04472363e4243.png

Evidently, our initial guess for μ0 is too high, so initial consumption too low.

We know this because we miss our KT+1=0 target on the high side.

Now we automate things with a search-for-a-good μ0 algorithm that stops when we hit the target Kt+1=0.

We use a bisection method.

We make an initial guess for C0 (we can eliminate μ0 because C0 is an exact function of μ0).

We know that the lowest C0 can ever be is 0 and that the largest it can be is initial output f(K0).

Guess C0 and shoot forward to T+1.

If KT+1>0, we take it to be our new lower bound on C0.

If KT+1<0, we take it to be our new upper bound.

Make a new guess for C0 that is halfway between our new upper and lower bounds.

Shoot forward again, iterating on these steps until we converge.

When KT+1 gets close enough to 0 (i.e., within an error tolerance bounds), we stop.

@jit
def bisection(pp, c0, k0, T=10, tol=1e-4, max_iter=500, k_ter=0, verbose=True):

    # initial boundaries for guess c0
    c0_upper = pp.f(k0) + (1 - pp.δ) * k0  
    c0_lower = 0

    i = 0
    while True:
        c_vec, k_vec = shooting(pp, c0, k0, T)
        error = k_vec[-1] - k_ter

        # check if the terminal condition is satisfied
        if np.abs(error) < tol:
            if verbose:
                print('Converged successfully on iteration ', i+1)
            return c_vec, k_vec

        i += 1
        if i == max_iter:
            if verbose:
                print('Convergence failed.')
            return c_vec, k_vec

        # if iteration continues, updates boundaries and guess of c0
        if error > 0:
            c0_lower = c0
        else:
            c0_upper = c0

        c0 = (c0_lower + c0_upper) / 2
def plot_paths(pp, c0, k0, T_arr, k_ter=0, k_ss=None, axs=None):

    if axs is None:
        fix, axs = plt.subplots(1, 3, figsize=(16, 4))
    ylabels = ['$c_t$', '$k_t$', r'$\mu_t$']
    titles = ['Consumption', 'Capital', 'Lagrange Multiplier']

    c_paths = []
    k_paths = []
    for T in T_arr:
        c_vec, k_vec = bisection(pp, c0, k0, T, k_ter=k_ter, verbose=False)
        c_paths.append(c_vec)
        k_paths.append(k_vec)

        μ_vec = pp.u_prime(c_vec)
        paths = [c_vec, k_vec, μ_vec]

        for i in range(3):
            axs[i].plot(paths[i])
            axs[i].set(xlabel='t', ylabel=ylabels[i], title=titles[i])

        # Plot steady state value of capital
        if k_ss is not None:
            axs[1].axhline(k_ss, c='k', ls='--', lw=1)

        axs[1].axvline(T+1, c='k', ls='--', lw=1)
        axs[1].scatter(T+1, paths[1][-1], s=80)

    return c_paths, k_paths

Now we can solve the model and plot the paths of consumption, capital, and Lagrange multiplier.

plot_paths(pp, 0.3, 0.3, [10]);
_images/a2c74a96ae84a1dcce83454e9c843be7c5bf5c848d23f471f56104d9c2115bbb.png

34.5. Setting Initial Capital to Steady State Capital#

When T+, the optimal allocation converges to steady state values of Ct and Kt.

It is instructive to set K0 equal to the limT+Kt, which we’ll call steady state capital.

In a steady state Kt+1=Kt=K¯ for all very large t.

Evalauating feasibility constraint (34.5) at K¯ gives

(34.17)#f(K¯)δK¯=C¯

Substituting Kt=K¯ and Ct=C¯ for all t into (34.14) gives

1=βu(C¯)u(C¯)[f(K¯)+(1δ)]

Defining β=11+ρ, and cancelling gives

1+ρ=1[f(K¯)+(1δ)]

Simplifying gives

f(K¯)=ρ+δ

and

K¯=f1(ρ+δ)

For production function (34.4), this becomes

αK¯α1=ρ+δ

As an example, after setting α=.33, ρ=1/β1=1/(19/20)1=20/1919/19=1/19, δ=1/50, we get

K¯=(33100150+119)671009.57583

Let’s verify this with Python and then use this steady state K¯ as our initial capital stock K0.

ρ = 1 / pp.β - 1
k_ss = pp.f_prime_inv(ρ+pp.δ)

print(f'steady state for capital is: {k_ss}')
steady state for capital is: 9.57583816331462

Now we plot

plot_paths(pp, 0.3, k_ss, [150], k_ss=k_ss);
_images/d25ade3e3e9413f35cb09bf2bec4f737fad4421486ba15984005c4cefb4cfb98.png

Evidently, with a large value of T, Kt stays near K0 until t approaches T closely.

Let’s see what the planner does when we set K0 below K¯.

plot_paths(pp, 0.3, k_ss/3, [150], k_ss=k_ss);
_images/da8defb5d09b607c4a39bce126fd6ee2e6511b8b93a106f63e7c073077598f53.png

Notice how the planner pushes capital toward the steady state, stays near there for a while, then pushes Kt toward the terminal value KT+1=0 when t closely approaches T.

The following graphs compare optimal outcomes as we vary T.

plot_paths(pp, 0.3, k_ss/3, [150, 75, 50, 25], k_ss=k_ss);
_images/897d6f9779120c366dcdc9c0a14df1b67f4c1d2a4fecbf3b62432d57e763d2e8.png

34.6. A Turnpike Property#

The following calculation indicates that when T is very large, the optimal capital stock stays close to its steady state value most of the time.

plot_paths(pp, 0.3, k_ss/3, [250, 150, 50, 25], k_ss=k_ss);
_images/87d3eb0dc3b0e75d1721c2fd9eea50069d5a27b793466db00c4594f7878cbc87.png

In the above graphs, different colors are associated with different horizons T.

Notice that as the horizon increases, the planner keeps Kt closer to the steady state value K¯ for longer.

This pattern reflects a turnpike property of the steady state.

A rule of thumb for the planner is

  • from K0, push Kt toward the steady state and stay close to the steady state until time approaches T.

The planner accomplishes this by adjusting the saving rate f(Kt)Ctf(Kt) over time.

Exercise 34.1

The turnpike property is independent of the initial condition K0 provided that T is sufficiently large.

Expand the plot_paths function so that it plots trajectories for multiple initial points using k0s = [k_ss*2, k_ss*3, k_ss/3].

Let’s calculate and plot the saving rate.

@jit
def saving_rate(pp, c_path, k_path):
    'Given paths of c and k, computes the path of saving rate.'
    production = pp.f(k_path[:-1])

    return (production - c_path) / production
def plot_saving_rate(pp, c0, k0, T_arr, k_ter=0, k_ss=None, s_ss=None):

    fix, axs = plt.subplots(2, 2, figsize=(12, 9))

    c_paths, k_paths = plot_paths(pp, c0, k0, T_arr, k_ter=k_ter, k_ss=k_ss, axs=axs.flatten())

    for i, T in enumerate(T_arr):
        s_path = saving_rate(pp, c_paths[i], k_paths[i])
        axs[1, 1].plot(s_path)

    axs[1, 1].set(xlabel='t', ylabel='$s_t$', title='Saving rate')

    if s_ss is not None:
        axs[1, 1].hlines(s_ss, 0, np.max(T_arr), linestyle='--')
plot_saving_rate(pp, 0.3, k_ss/3, [250, 150, 75, 50], k_ss=k_ss)
_images/3119fcc18d98a535163e48c47c263d746f70dc0bde210b4e8a089e7d75e5fe65.png

34.7. A Limiting Infinite Horizon Economy#

We want to set T=+.

The appropriate thing to do is to replace terminal condition (34.12) with

limT+βTu(CT)KT+1=0,

a condition that will be satisfied by a path that converges to an optimal steady state.

We can approximate the optimal path by starting from an arbitrary initial K0 and shooting towards the optimal steady state K at a large but finite T+1.

In the following code, we do this for a large T and plot consumption, capital, and the saving rate.

We know that in the steady state that the saving rate is constant and that s¯=f(K¯)C¯f(K¯).

From (34.17) the steady state saving rate equals

s¯=δK¯f(K¯)

The steady state saving rate S¯=s¯f(K¯) is the amount required to offset capital depreciation each period.

We first study optimal capital paths that start below the steady state.

# steady state of saving rate
s_ss = pp.δ * k_ss / pp.f(k_ss)

plot_saving_rate(pp, 0.3, k_ss/3, [130], k_ter=k_ss, k_ss=k_ss, s_ss=s_ss)
_images/60e320700e05982742b2a64aeed91f2c72c19f3eeb31b2e6440a79b54ec871b5.png

Since K0<K¯, f(K0)>ρ+δ.

The planner chooses a positive saving rate that is higher than the steady state saving rate.

Note that f(K)<0, so as K rises, f(K) declines.

The planner slowly lowers the saving rate until reaching a steady state in which f(K)=ρ+δ.

34.8. Stable Manifold and Phase Diagram#

We now describe a classic diagram that describes an optimal (Kt+1,Ct) path.

The diagram has K on the ordinate axis and C on the coordinate axis.

Given an arbitrary and fixed K, a fixed point C of the consumption Euler equation (34.15) satisfies

C=C(β[f(f(K)+(1δ)KC)+(1δ)])1/γ

which implies

(34.18)#C=f(K)+(1δ)Kf1(1β(1δ))C~(K)

A positive fixed point C=C~(K) exists only if f(K)+(1δ)Kf1(1β(1δ))>0

@jit
def C_tilde(K, pp):

    return pp.f(K) + (1 - pp.δ) * K - pp.f_prime_inv(1 / pp.β - 1 + pp.δ)

Next note that given a time-invariant arbitrary C, a fixed point K of the feasibility condition (34.5) solves the following equation

K=f(K)+(1δK)C.

A fixed point of the above equation is described by a function

(34.19)#K=K~(C)
@jit
def K_diff(K, C, pp):
    return pp.f(K) - pp.δ * K - C

@jit
def K_tilde(C, pp):

    res = brentq(K_diff, 1e-6, 100, args=(C, pp))

    return res.root

A steady state (Ks,Cs) is a pair (K,C) that satisfies both equations (34.18) and (34.19).

It is thus the intersection of the two curves C~ and K~ that we’ll plot in Figure Fig. 34.1 below.

We can compute Ks by solving the equation Ks=K~(C~(Ks))

@jit
def K_tilde_diff(K, pp):

    K_out = K_tilde(C_tilde(K, pp), pp)

    return K - K_out
res = brentq(K_tilde_diff, 8, 10, args=(pp,))

Ks = res.root
Cs = C_tilde(Ks, pp)

Ks, Cs
(9.575838163314447, 1.9160839808123402)

We can use the shooting algorithm to compute trajectories that approach (Ks,Cs).

For a given K, let’s compute C and K for a large T , e.g., =200.

We compute C0 by the bisection algorithm that assures that KT=Ks.

Let’s compute two trajectories towards (Ks,Cs) that start from different sides of Ks: K¯0=1e3<Ks<K¯1=15.

c_vec1, k_vec1 = bisection(pp, 5, 15, T=200, k_ter=Ks)
c_vec2, k_vec2 = bisection(pp, 1e-3, 1e-3, T=200, k_ter=Ks)
Converged successfully on iteration  46
Converged successfully on iteration  51

The following code generates Figure Fig. 34.1, which is patterned on a graph that appears on page 411 of [Intriligator, 2002].

Figure Fig. 34.1 is a classic “phase plane” with “state” variable K on the ordinate axis and “co-state” variable C on the coordinate axis.

Figure Fig. 34.1 plots three curves:

  • the blue line graphs C=C~(K) of fixed points described by equation (34.18).

  • the red line graphs K=K~(C) of fixed points described by equation (34.19)

  • the green line graphs the stable traced out by paths that converge to the steady state starting from an arbitrary K0 at time 0.

    • for a given K0, the shooting algorithm sets C0 to the coordinate on the green line in order to initiate a path that converges to the optimal steady state

    • the arrows on the green line show the direction in which dynamics (34.16) push successive (Kt+1,Ct) pairs.

In addition to the three curves, Figure Fig. 34.1 plots arrows that point where the dynamics (34.16) drive the system when, for a given K0, C0 is not on the stable manifold depicted in the green line.

  • If C0 is set below the green line for a given K0, too much capital is accumulated

  • If C0 is set above the green line for a given K0, too little capital is accumulated

Hide code cell source
fig, ax = plt.subplots(figsize=(7, 5))

K_range = np.arange(1e-1, 15, 0.1)
C_range = np.arange(1e-1, 2.3, 0.1)

# C tilde
ax.plot(K_range, [C_tilde(Ks, pp) for Ks in K_range], color='b')
ax.text(11.8, 4, r'$C=\tilde{C}(K)$', color='b')

# K tilde
ax.plot([K_tilde(Cs, pp) for Cs in C_range], C_range, color='r')
ax.text(2, 1.5, r'$K=\tilde{K}(C)$', color='r')

# stable branch
ax.plot(k_vec1[:-1], c_vec1, color='g')
ax.plot(k_vec2[:-1], c_vec2, color='g')
ax.quiver(k_vec1[5], c_vec1[5],
          k_vec1[6]-k_vec1[5], c_vec1[6]-c_vec1[5],
          color='g')
ax.quiver(k_vec2[5], c_vec2[5],
          k_vec2[6]-k_vec2[5], c_vec2[6]-c_vec2[5],
          color='g')
ax.text(12, 2.5, r'stable branch', color='g')

# (Ks, Cs)
ax.scatter(Ks, Cs)
ax.text(Ks-1.2, Cs+0.2, '$(K_s, C_s)$')

# arrows
K_range = np.linspace(1e-3, 15, 20)
C_range = np.linspace(1e-3, 7.5, 20)
K_mesh, C_mesh = np.meshgrid(K_range, C_range)

next_K, next_C = pp.next_k_c(K_mesh, C_mesh)
ax.quiver(K_range, C_range, next_K-K_mesh, next_C-C_mesh)

# infeasible consumption area
ax.text(0.5, 5, "infeasible\n consumption")

ax.set_ylim([0, 7.5])
ax.set_xlim([0, 15])

ax.set_xlabel('$K$')
ax.set_ylabel('$C$')

plt.show()
_images/b96c320263573c2a1794e59d680a75b1921fc3cc6ddeb198f543747b9ff9a37d.png

Fig. 34.1 Stable Manifold and Phase Plane#

34.9. Concluding Remarks#

In Cass-Koopmans Competitive Equilibrium, we study a decentralized version of an economy with exactly the same technology and preference structure as deployed here.

In that lecture, we replace the planner of this lecture with Adam Smith’s invisible hand.

In place of quantity choices made by the planner, there are market prices that are set by a deus ex machina from outside the model, a so-called invisible hand.

Equilibrium market prices must reconcile distinct decisions that are made independently by a representative household and a representative firm.

The relationship between a command economy like the one studied in this lecture and a market economy like that studied in Cass-Koopmans Competitive Equilibrium is a foundational topic in general equilibrium theory and welfare economics.

34.9.1. Exercise#

Exercise 34.2

  • Plot the optimal consumption, capital, and saving paths when the initial capital level begins at 1.5 times the steady state level as we shoot towards the steady state at T=130.

  • Why does the saving rate respond as it does?