Undergraduate Computational Macro – Deterministic Dynamics and Introduction to Growth Models

Summary of Equations

Exogenous $z_{t}$ sequence. e.g., $z_{t + 1} / z_{t} = 1 + g_{z}$ given some initial $z_{0}$
Population growth $N_{t + 1} / N_{t} = 1 + g_{N}$ given some initial $N_{0}$

$k_{t + 1} = \frac{1}{1 + g_{N}} [(1 - δ) k_{t} + s z_{t} f (k_{t})], given k_{0}$
- Output per capita $y_{t} = z_{t} f (k_{t})$
- Consumption per capita $c_{t} = (1 - s) y_{t} = (1 - s) z_{t} f (k_{t})$
Wages $w_{t} = (1 - α) z_{t} k_{t}^{α}$ and rental rate of capital $r_{t} = α z_{t} k_{t}^{α - 1}$
Steady state capital $\bar{k} = {(\frac{s \bar{z}}{g_{N} + δ})}^{\frac{1}{1 - α}}$ if $g_{z} = 0$ and $z_{0} = \bar{z}$

45 Degree Diagram

With a fixed point $k_{t + 1} = h (k_{t})$ note that a fixed point is when $\bar{k} = h (\bar{k})$
We can plot the dynamics of the sequence comparing the functions to the 45 degree line where that occurs
This diagram will help us interpret stability and convergence

Iteration

First, lets write a general function to iterate a (univariate) map

function iterate_map(f, x0, T)
    x = zeros(T + 1)
    x[1] = x0
    for t in 2:(T + 1)
        x[t] = f(x[t - 1])
    end
    return x
end

iterate_map (generic function with 1 method)

Plotting the Dynamics

function plot45(f, xmin, xmax, x0, T; num_points = 100, label = L"h(k)",
                xlabel = "k", size = (600, 500))
    # Plot the function and the 45 degree line
    x_grid = range(xmin, xmax, num_points)
    plt = plot(x_grid, f.(x_grid); xlim = (xmin, xmax), ylim = (xmin, xmax),
               linecolor = :black, lw = 2, label, size)
    plot!(x_grid, x_grid; linecolor = :blue, lw = 2, label = nothing)

    # Iterate map and add ticks
    x = iterate_map(f, x0, T)
    if !isnothing(xlabel) && T > 1
      xticks!(x, [L"%$(xlabel)_{%$i}" for i in 0:T])
      yticks!(x, [L"%$(xlabel)_{%$i}" for i in 0:T])
    end

    # Plot arrows and dashes
    for i in 1:T
        plot!([x[i], x[i]], [x[i], x[i + 1]], arrow = :closed, linecolor = :black,
              alpha = 0.5, label = nothing)
        plot!([x[i], x[i + 1]], [x[i + 1], x[i + 1]], arrow = :closed,
              linecolor = :black, alpha = 0.5, label = nothing)
        plot!([x[i + 1], x[i + 1]], [0, x[i + 1]], linestyle = :dash,
              linecolor = :black, alpha = 0.5, label = nothing)
    end
    plot!([x[1], x[1]], [0, x[1]], linestyle = :dash, linecolor = :black,
          alpha = 0.5, label = nothing)
end

plot45 (generic function with 1 method)

Fixed Points

# Implementation of our capital dynamics
h(k; p) = (1 / (1 + p.g_N)) * (
   p.s * p.z_bar * k^p.alpha
   + (1 - p.delta) * k)
k_bar(p) = (p.s * p.z_bar /
      (p.g_N + p.delta))^(1/(1-p.alpha))
p = (;z_bar = 2.0, s = 0.3, alpha = 0.3,
     delta = 0.4, g_N = 0.0)
@show k_bar(p)
@show h(k_bar(p); p)
@show h(0.0;p);

k_bar(p) = 1.7846741842265788
h(k_bar(p); p) = 1.7846741842265788
h(0.0; p) = 0.0

45 Degree Diagram for Solow

k_min = 0.0
k_max = 3.0
k_0 = 0.25
T = 6
plot45(k -> h(k; p), k_min, k_max, k_0,
       T;label = L"k_{t+1} = h(k_t)")

Transition Dynamics

k_vals = iterate_map(k -> h(k; p), k_0, T)
plot(0:T, k_vals; label =[L"k_t" nothing],
     title=L"Dynamics from $k_0 = %$k_0$",
     seriestype = [:scatter, :line],
     xlabel = "t", size=(600, 400))

Rental Rate of Capital

Why does it decrease?

r(k;p) = p.alpha * p.z_bar * k^(p.alpha - 1)
w(k;p) = (1 - p.alpha) * p.z_bar * k^(p.alpha)

@show k_0
plot(0:T, r.(k_vals; p); label = L"r_t",
     title="Real Rental Rate of Capital",
     xlabel = "t", size=(600, 400))
hline!([r(k_bar(p);p)];linestyle=:dash,
       label=L"r(\bar{k})", alpha=0.5)

k_0 = 0.25

Wages

Why does it increase?

@show k_0
plot(0:T, w.(k_vals; p); label = L"w_t",
     title="Real Wages",
      xlabel = "t", size=(600, 400))
hline!([w(k_bar(p);p)];linestyle=:dash,
      label=L"w(\bar{k})", alpha=0.5)

k_0 = 0.25

Above the Steady State?

k_0 = 2.5
plot45(k -> h(k; p), k_min, k_max, k_0,
       T;label = L"k_{t+1} = h(k_t)")

Transition Dynamics

k_vals = iterate_map(k -> h(k; p), k_0, T)
plot(0:T, k_vals; label =[L"k_t" nothing],
     title=L"Dynamics from $k_0 = %$k_0$",
     seriestype = [:scatter, :line],
     xlabel = "t", size=(600, 400))

Rental Rate of Capital

Why does it increase?

@show k_0
plot(0:T, r.(k_vals; p); label = L"r_t",
     title="Real Rental Rate of Capital",
     xlabel = "t", size=(600, 400))
hline!([r(k_bar(p);p)];linestyle=:dash,
       label=L"r(\bar{k})", alpha=0.5)

k_0 = 2.5

Wages

Why does it decrease?

@show k_0
plot(0:T, w.(k_vals; p); label = L"w_t",
     title="Real Wages",
      xlabel = "t", size=(600, 400))
hline!([w(k_bar(p);p)];linestyle=:dash,
      label=L"w(\bar{k})", alpha=0.5)

k_0 = 2.5

At Zero Capital

k_0 = 0.0
@show h(k_0; p)
plot45(k -> h(k; p), k_min, k_max, k_0,
       T;label = L"k_{t+1} = h(k_t)")

h(k_0; p) = 0.0

Perturb Zero Capital

k_0 = 0.0001
@show h(k_0; p)
plot45(k -> h(k; p), k_min, k_max, k_0,
       T;label = L"k_{t+1} = h(k_t)")

h(k_0; p) = 0.03791744066881159

Transition Dynamics

k_vals = iterate_map(k -> h(k; p), k_0, T)
plot(0:T, k_vals; label =[L"k_t" nothing],
     title=L"Dynamics from $k_0 = %$k_0$",
     seriestype = [:scatter, :line],
     xlabel = "t", size=(600, 400))

Marginal Product of Capital

Recall that the Marginal Product of Capital (MPK) is $z_{t} \frac{\partial F (K_{t}, N_{t})}{\partial K_{t}} = z_{t} f^{'} (K_{t} / N_{t}) = α z_{t} k_{t}^{α - 1}$
How does the MPK change as the economy grows?

Visualizing the MPK

Strictly positive, monotonically decreasing, asymptote at $k = 0$

# No g_N here, that is in dynamics
MPK(k; p) = p.alpha * p.z_bar * k^(p.alpha - 1)
k_vals = range(0.0, 10.0; length=100)
plot(k_vals, MPK.(k_vals; p);
     label = L"\alpha\bar{z}k^{\alpha-1}",
     xlabel = "k",
     title="Marginal Product of Capital",
     size=(600, 500), ylim=(0.0, 2.0))

Visualizing the MPK Transition

Converges to a point where capital can’t accumulate without smaller $s$

k_0 = 0.25
T = 10
k_vals = iterate_map(k -> h(k; p), k_0, T)
plot(0:T, MPK.(k_vals; p);
     label = L"\alpha\bar{z}k_t^{\alpha-1}",
     xlabel = "t",
     title="MPK Transition Dynamics",
     size=(600, 500), ylim=(0.0, 2.0))

Gradients of $h (k)$

$\nabla h (k; p) = \frac{1}{1 + g_{N}} (α s z_{t} k^{1 - α} + 1 - δ)$

Recall fixed points are $k^{*} = h (k^{*})$
Contraction mappings are a global property where points get closer
Locally, we can ask the same question. Gradients help us understand whether points expand or contract

Plotting $\nabla h (k)$

h_k(k; p) = (1 / (1 + p.g_N)) * (
   p.alpha * p.s * p.z_bar * k^(p.alpha-1)
   + 1 - p.delta)
k_vals = range(0.0, 10.0; length=100)
plot(k_vals, h_k.(k_vals; p);
     label = L"\nabla h(k)",
     xlabel = "k",
     title=L"Gradient of $h(k)$",
     size=(600, 500), ylim=(0.0, 2.0))
vline!([k_bar(p)];label=L"h(\bar{k})=\bar{k}",
       linestyle=:dash)
vline!([0.0];label=L"h(0)=0", linestyle=:dash)
hline!([1.0];linestyle=:dot, label=nothing,
       alpha=0.5, lw=1)

Local and Global Stability

Consider a fixed point $\bar{k} = h (\bar{k})$
$\bar{k}$ is local stable if there exists an $ϵ > 0$
- For all $| k - \bar{k} | < ϵ$ , $lim_{T \to \infty} h^{T} (k) = \bar{k}$
i.e., starting close to the steady state it converges to steady state
- Global stability if $ϵ = \infty$
- Given some continuity assumptions, can show that $\nabla h (\bar{k}) < 1$ is sufficient
Solow: $k = 0$ blows up, so this is only locally stable for $ϵ = \bar{k}$

Non-convexity in Production

Let $f (k) = - \frac{1}{3} k^{3} + k^{2} + \frac{1}{3} k, h (k) = s f (k) + (1 - δ) k$

function h_multi(k;s=0.95, delta=0.99)
  return s*(-k^3/3 + k^2 + k/3) + (1 - delta)*k
end
k_min = 0.0
k_max = 2.2
k_vals = range(k_min, k_max; length=100)
plot(k_vals, h_multi.(k_vals);
     label = L"h(k)",
     xlabel = "k",
     title="Multiple Crossings",
     size=(600, 500))
plot!(k_vals, k_vals; label=nothing,
      style=:dash)

Fixed Points

There are now 3 fixed points
In this case we can eyeball the graph for initial conditions

h_multi_vec(k_vec) = [h_multi(k_vec[1])] # pack/unpack as vector
k_star_1 = fixedpoint(h_multi_vec, [0.0]).zero[1]
k_star_2 = fixedpoint(h_multi_vec, [1.0]).zero[1]
k_star_3 = fixedpoint(h_multi_vec, [2.0]).zero[1]
@show k_star_1, k_star_2, k_star_3;

(k_star_1, k_star_2, k_star_3) = (0.0, 1.1483123394919963, 1.8516876604860064)

Detour into Roots of Equations

A root or zero of a function $\hat{h} (\cdot)$ is a point $k^{*}$ such that $\hat{h} (k^{*}) = 0$
Fixed points $k^{*} = h (k^{*})$ is a root of the equation $\hat{h} (k) \equiv h (k) - k$
Alternative algorithms for univariate solvers can bracket a solution

sol = nlsolve(k -> h_multi_vec(k) .- k, [1.0]) # vectorized, NLsolve
@show sol.zero
@show find_zero(k -> h_multi(k) - k, (1.0, 1.5)); # bracketed, from Roots

sol.zero = [1.1483123394996155]
find_zero((k->begin
            #= In[23]:3 =#
            h_multi(k) - k
        end), (1.0, 1.5)) = 1.1483123395307762

Steady States and Gradients

There are now 3 fixed points, lets look at the gradients
Remember that $\nabla h (\bar{k}) < 1$ is sufficient for stability

h_multi_k(k;s=0.95, delta=0.99) = s*(-k^2 + 2 * k + 1/3) + (1 - delta)
@show k_star_1, h_multi_k(k_star_1)
@show k_star_2, h_multi_k(k_star_2)
@show k_star_3, h_multi_k(k_star_3);

(k_star_1, h_multi_k(k_star_1)) = (0.0, 0.32666666666666666)
(k_star_2, h_multi_k(k_star_2)) = (1.1483123394919963, 1.2557699441233567)
(k_star_3, h_multi_k(k_star_3)) = (1.8516876604860064, 0.5875633891937458)

In Increasing Returns Region

k_0 = 0.8
T = 5
plot45(h_multi, k_min, k_max, k_0,
       T;label = L"k_{t+1} = h(k_t)")

In Decreasing Returns Region

k_0 = 1.5
T = 5
plot45(h_multi, k_min, k_max, k_0,
       T;label = L"k_{t+1} = h(k_t)")

In Decreasing Returns Region

k_0 = 2.0
T = 3
plot45(h_multi, k_min, k_max, k_0,
       T;label = L"k_{t+1} = h(k_t)")

Visualizing the MPK

function MPK_multi(k;s=0.95, delta=0.99)
    return s*(-k^2 + 2 * k + 1/3)
end
k_vals = range(0.0, 2.0; length=100)
plot(k_vals, MPK_multi.(k_vals);
     label = "MPK",
     xlabel = "k",
     title="Marginal Product of Capital",
     size=(600, 500), ylim=(0.0, 2.0))