Kydland & Prescott (1982): Real Business Cycles

David Evans

2026-05-28

Introduction

From Micro to Macro

Where We’ve Been

• McCall (1970): an unemployed worker searches for a job
  • solve for a reservation wage
• Rust (1987): a bus superintendent decides when to replace engines
  • solve for replacement probabilities

• Both models feature a single agent making decisions in a micro environment
• Today we scale up: a representative agent in a macroeconomy

Kydland & Prescott (1982)

• “Time to Build and Aggregate Fluctuations”
  • Econometrica and 10000+ citations
• Won the 2004 Nobel Prize in Economics

• Key question: Can a competitive equilibrium model with technology shocks explain business cycles?

The RBC Revolution

Before KP (1982), business cycles were thought to require market failures (sticky prices, coordination failures). KP showed that a perfectly competitive economy hit by productivity shocks could generate realistic fluctuations.

The Big Idea

• Business cycles as the optimal response of a competitive economy to random productivity shocks
• No market failures, no government intervention needed to explain fluctuations

• The core model: a stochastic neoclassical growth model
  • Households save and consume optimally
  • Firms produce using capital and labor
  • Productivity fluctuates randomly

Connection to Previous Lectures

The RBC model combines: dynamic programming (McCall, Rust), Markov chains (productivity shocks), and linear algebra (matrix operations for expectations).

A New Computational Challenge

• In McCall and Rust, the state space was discrete (wage offers, mileage bins)
• We stored the value function as a vector and looked up values directly

• In the RBC model, capital \(k\) is continuous — we can’t just index into a vector!

• We’ll introduce linear interpolation
  • a tool for evaluating functions at arbitrary points, not just grid points
• This unlocks a whole class of models with continuous state spaces

Roadmap

What We’ll Cover

1. The Model — households, firms, competitive equilibrium, and the social planner’s problem
2. Calibration — matching long-run moments and the productivity process
3. Linear Interpolation — a new tool for continuous state spaces
4. Solving the Model — value function iteration with interpolation
5. Simulation — generating artificial business cycles with DataFrames
6. Linearized Solution — differentiating the Euler equation for a fast, analytical policy rule

The Stochastic Growth Model

The Household

The Household’s Problem

• A representative household lives forever and maximizes:

\[ \max_{\{c_t, k_{t+1}\}} \; \mathbb{E}_0 \sum_{t=0}^{\infty} \beta^t \log(c_t) \]

• \(c_t\): consumption in period \(t\)
• \(\beta \in (0,1)\): discount factor
• \(\log(c_t)\): log utility
  • more consumption is better, but with diminishing returns

Income and Ownership

• The household owns capital \(k_t\) and supplies one unit of labor each period (inelastically)
• It earns two types of income:

1. Rental income: \(r_t k_t\) from renting capital to firms
2. Wage income: \(w_t\) from supplying labor

• The household takes prices \((r_t, w_t)\) as given: it is a price taker in competitive markets

The Budget Constraint

• Each period, the household faces the budget constraint:

\[ c_t + k_{t+1} = w_t + r_t \, k_t + (1 - \delta) \, k_t \]

• Left side: uses of funds
  • consumption today + capital carried into tomorrow
• Right side: sources of funds
  • wages + capital income + undepreciated capital

• The household chooses \(c_t\) and \(k_{t+1}\) to maximize lifetime utility subject to this constraint

Simplification

We fix labor supply at \(\ell_t = 1\) (inelastic). The full KP (1982) model includes labor-leisure choice and time-to-build for investment, which generate additional dynamics. The core insights carry through without these.

The Firm

The Firm’s Problem

• A representative firm produces output using capital and labor:

\[ y_t = z_t \, k_t^\alpha \, \ell_t^{1-\alpha} \]

• \(z_t\): total factor productivity (TFP): the random shock driving business cycles
• \(\alpha \in (0,1)\): capital’s share of output (typically \(\approx 1/3\))
• Constant returns to scale (CRS): doubling both inputs doubles output

Profit Maximization

• The firm rents capital at rate \(r_t\) and hires labor at wage \(w_t\):

\[ \pi_t = z_t \, k_t^\alpha \, \ell_t^{1-\alpha} - r_t \, k_t - w_t \, \ell_t \]

• First-order conditions (set marginal product = factor price):

\[ r_t = \alpha \, z_t \, k_t^{\alpha - 1} \, \ell_t^{1-\alpha} \qquad\qquad w_t = (1-\alpha) \, z_t \, k_t^{\alpha} \, \ell_t^{-\alpha} \]

Zero Profits

With CRS, paying each factor its marginal product exactly exhausts output: \(r_t k_t + w_t \ell_t = y_t\). The firm earns zero profits in equilibrium

Factor Prices in Equilibrium

• Since the household supplies \(\ell_t = 1\) unit of labor, the factor prices simplify to:

\[ r_t = \alpha \, z_t \, k_t^{\alpha - 1} \qquad\qquad w_t = (1-\alpha) \, z_t \, k_t^{\alpha} \]

• The rental rate falls with capital (diminishing returns)
• The wage rises with capital (more capital per worker \(\rightarrow\) higher productivity)

• Both move with \(z_t\): when productivity is high, all factor prices rise

Technology Shocks

The Driving Force

• Productivity follows an AR(1) process in logs:

\[ \log z_{t+1} = \rho \, \log z_t + \sigma_\varepsilon \, \varepsilon_{t+1}, \qquad \varepsilon_{t+1} \sim N(0, 1) \]

• \(\rho \in (0,1)\): persistence of technology shocks (typically 0.9–0.95)
• \(\sigma_\varepsilon\): volatility of technology innovations

• High \(\rho\) means a good shock today predicts good productivity tomorrow: persistence drives business cycle dynamics

Connecting to Markov Chains

We’ll discretize this continuous AR(1) process into a finite-state Markov chain. This lets us use our Markov chain tools!

Discretizing the AR(1): Rouwenhorst Method

• We approximate the continuous AR(1) with an \(N\)-state Markov chain
• The Rouwenhorst (1995) method:
  1. Create an evenly spaced grid for \(\log z\) spanning \(\pm \sqrt{N-1} \cdot \sigma_z\)
  2. Build the transition matrix recursively

Why Rouwenhorst?

Rouwenhorst is simpler to implement than Tauchen (1986) (an alternative method) and performs better for highly persistent processes (\(\rho\) close to 1): exactly the case in RBC models.

The rouwenhorst Function

function rouwenhorst(N, ρ, σ_ε)
    σ_z = σ_ε / sqrt(1 - ρ^2)                     # unconditional std dev
    p = (1 + ρ) / 2                                 # persistence parameter
    P = [p 1-p; 1-p p]                             # base case: N = 2
    for n in 3:N                                    # recursively expand
        z = zeros(n - 1)
        P = p     * [P z; z' 0] + (1-p) * [z P; 0 z'] +
            (1-p) * [z' 0; P z] + p     * [0 z'; z P]
        P[2:end-1, :] ./= 2                        # normalize interior rows
    end
    z_grid = collect(LinRange(-sqrt(N-1) * σ_z,     # evenly spaced grid
                               sqrt(N-1) * σ_z, N))
    return z_grid, P
end

rouwenhorst (generic function with 1 method)

Testing the Rouwenhorst Approximation

log_z_grid, P_z = rouwenhorst(7, 0.90, 0.02)
z_grid = exp.(log_z_grid)                           # convert to levels
round.(z_grid, digits=4)

7-element Vector{Float64}:
 0.8937
 0.9278
 0.9632
 1.0
 1.0382
 1.0778
 1.1189

round.(P_z, digits=3)

7×7 Matrix{Float64}:
 0.735  0.232  0.031  0.002  0.0    0.0    0.0
 0.039  0.745  0.195  0.02   0.001  0.0    0.0
 0.002  0.078  0.751  0.156  0.012  0.0    0.0
 0.0    0.006  0.117  0.753  0.117  0.006  0.0
 0.0    0.0    0.012  0.156  0.751  0.078  0.002
 0.0    0.0    0.001  0.02   0.195  0.745  0.039
 0.0    0.0    0.0    0.002  0.031  0.232  0.735

• High diagonal values: the economy tends to stay near its current productivity level

Visualizing the Transition Matrix

using Plots
heatmap(P_z, xlabel="Next State", ylabel="Current State",
        title="Transition Matrix for Productivity", color=:viridis, yflip=true)

Competitive Equilibrium

Definition

Competitive Equilibrium

A competitive equilibrium is a sequence of prices \(\{r_t, w_t\}_{t=0}^{\infty}\) and allocations \(\{c_t, k_{t+1}\}_{t=0}^{\infty}\) such that:

1. Household optimizes: given prices, the household maximizes expected lifetime utility subject to its budget constraint
2. Firm optimizes: given prices, the firm’s factor demands satisfy \(r_t = \text{MPK}\) and \(w_t = \text{MPL}\)
3. Markets clear: the labor market clears (\(\ell_t = 1\)) and the goods market clears (\(c_t + k_{t+1} - (1-\delta)k_t = y_t\))

From Budget Constraint to Resource Constraint

• Substituting the zero-profit factor prices into the household’s budget constraint:

\[ c_t + k_{t+1} = \underbrace{(1-\alpha) \, z_t \, k_t^\alpha}_{w_t} + \underbrace{\alpha \, z_t \, k_t^{\alpha-1}}_{r_t} \cdot k_t + (1-\delta) \, k_t \]

\[ c_t + k_{t+1} = z_t \, k_t^\alpha + (1-\delta) \, k_t \]

• The budget constraint collapses to the aggregate resource constraint!
• Wages + rental income = total output (zero profits under CRS)

Why This Matters

• In the competitive equilibrium, we would need to track:
  • Household optimization (Euler equation)
  • Firm optimization (factor prices)
  • Market clearing conditions

• That’s a lot of moving parts! Is there a shortcut?

Yes! The First Welfare Theorem!

The First Welfare Theorem

From Equilibrium to Planning

• The First Welfare Theorem states: under standard conditions (complete markets, no externalities, price-taking behavior), a competitive equilibrium allocation is Pareto optimal

• A Pareto optimal allocation is exactly what a benevolent social planner would choose

The Shortcut

The CE allocation = the planner’s solution. So instead of tracking prices and market clearing, we can solve the much simpler planner’s problem and recover the competitive equilibrium allocations directly.

• To find equilibrium prices, we can always back them out: \(r_t = \alpha z_t k_t^{\alpha-1}\) and \(w_t = (1-\alpha) z_t k_t^\alpha\)

The Social Planner’s Problem

The Sequential Problem

• Given initial conditions \((k_0, z_0)\), the planner chooses entire sequences of consumption and capital to maximize expected lifetime utility:

\[ \max_{\{c_t,\, k_{t+1}\}_{t=0}^{\infty}} \; \mathbb{E}_0 \sum_{t=0}^{\infty} \beta^t \log(c_t) \]

subject to, for every \(t \geq 0\), \[ c_t + k_{t+1} = z_t \, k_t^\alpha + (1-\delta) \, k_t, \qquad c_t \geq 0, \quad k_{t+1} \geq 0 \]

• Productivity follows the finite-state Markov chain from Rouwenhorst:
  • let \(\bar z_1, \ldots, \bar z_S\) denote the \(S\) possible productivity values
  • let \(s_t \in \{1, \ldots, S\}\) be the index, so \(z_t = \bar z_{s_t}\)
  • transition matrix \(P_{s, s'} = \Pr(s_{t+1} = s' \mid s_t = s)\)

• This is an infinite-dimensional optimization problem
  • one choice per period, forever, across every possible realization of \(\{s_t\}\)
• Bellman equations help us break this giant problem into something manageable!

Splitting: Today vs. the Future

• Pull the \(t = 0\) term out of the sum and factor \(\beta\) from the rest:

\[ \mathbb{E}_0 \sum_{t=0}^{\infty} \beta^t \log(c_t) = \log(c_0) + \beta \, \mathbb{E}_0 \sum_{t=1}^{\infty} \beta^{t-1} \log(c_t) \]

• The planner’s problem becomes a choice today plus a continuation problem starting tomorrow:

\[ \max_{c_0,\, k_1} \left\{ \log(c_0) + \beta \, \mathbb{E}_0 \left[\; \max_{\{c_t,\, k_{t+1}\}_{t=1}^\infty} \sum_{t=1}^\infty \beta^{t-1} \log(c_t) \;\right] \right\} \]

• The inner problem has the same form as the original, just starting from \((k_1, s_1)\) instead of \((k_0, s_0)\)
  • the Markov property is key: tomorrow’s distribution depends only on \(s_1\), not the whole history \(\{s_0, s_1\}\)
• Define \(V(k, s)\) as the value of the planner’s problem starting with capital \(k\) when productivity is \(\bar z_s\)
  • then the inner max equals \(V(k_1, s_1)\)

The Bellman Equation

• Substituting \(V(k_1, s_1)\) for the continuation problem and dropping time-0 subscripts gives a recursive problem in just two variables:
• State variables: \((k, s)\) — current capital and productivity index
• Control variable: \(k'\) (which pins down \(c = \bar z_s \, k^\alpha + (1-\delta)k - k'\))

\[ V(k, s) = \max_{k'} \left\{ \log(c) + \beta \, \mathbb{E}\left[V(k', s') \mid s\right] \right\} \]

where \(c = \bar z_s \, k^\alpha + (1-\delta) k - k' > 0\)

• The conditional expectation is a sum over the Markov chain: \[ \mathbb{E}\left[V(k', s') \mid s\right] = \sum_{s' = 1}^{S} P_{s, s'} \, V(k', s') \]

The Deterministic Steady State

What Happens Without Shocks?

• Before solving the full stochastic model, consider: what if productivity were constant at \(z = 1\)?
• In the long run, capital settles at a steady state \(k_{ss}\) where \(k_{t+1} = k_t = k_{ss}\)

• At the steady state, the household consumes all output net of depreciation:
  • Investment exactly replaces worn-out capital: \(i_{ss} = \delta \, k_{ss}\)
  • Consumption: \(c_{ss} = k_{ss}^\alpha - \delta \, k_{ss}\)

• Understanding the steady state helps us in three ways:
  1. It tells us where capital should be on average
  2. It guides our grid construction (center the grid around \(k_{ss}\))
  3. It provides a starting point for simulation

Deriving \(k_{ss}\): A Perturbation Argument

• Suppose the planner is at the steady state and considers saving one extra unit of capital today

• Cost: consuming one less unit today reduces utility by \(\frac{1}{c_{ss}}\)

• Benefit: that extra unit of capital produces \(\alpha \, k_{ss}^{\alpha - 1}\) extra output tomorrow, plus \((1-\delta)\) undepreciated capital comes back
  • valued at \(\frac{1}{c_{ss}}\) and discounted by \(\beta\)

\[ \text{Benefit} = \beta \cdot \frac{\alpha \, k_{ss}^{\alpha - 1} + 1 - \delta}{c_{ss}} \]

The Steady State Condition

• At the optimum, the planner cannot improve by saving more or less: cost = benefit

\[ \frac{1}{c_{ss}} = \beta \cdot \frac{\alpha \, k_{ss}^{\alpha - 1} + 1 - \delta}{c_{ss}} \]

• The \(\frac{1}{c_{ss}}\) terms cancel, leaving:

\[ 1 = \beta \left( \alpha \, k_{ss}^{\alpha - 1} + 1 - \delta \right) \]

• Solving for \(k_{ss}\):

\[ \boxed{k_{ss} = \left(\frac{\alpha \beta}{1 - \beta(1-\delta)}\right)^{\frac{1}{1-\alpha}}} \]

Interpreting the Steady State

• The term \(\alpha \, k_{ss}^{\alpha - 1}\) is the marginal product of capital
  • the steady-state rental rate \(r_{ss}\)

• Our condition says \(\beta(r_{ss} + 1 - \delta) = 1\), or:

\[ r_{ss} = \frac{1}{\beta} - (1-\delta) \]

• The more patient the household (\(\beta\) closer to 1), the lower the equilibrium interest rate
  • patient agents save more, driving down returns

The Logic of the Steady State

The planner keeps accumulating capital until the net return (\(r_{ss} - \delta\)) falls to match the household’s rate of impatience (\(1/\beta - 1\)). At that point, saving one more unit costs exactly as much as it benefits

Setup in Julia

Loading Packages

using Plots
using LinearAlgebra
using Statistics
using Distributions
using Random
using Optim
using Interpolations
using DataFrames

The RBCModel Struct

• A mutable struct holds the model parameters and the objects we’ll compute from them:

@kwdef mutable struct RBCModel
    β::Float64 = 0.96                           # discount factor
    α::Float64 = 0.33                           # capital share
    δ::Float64 = 0.08                           # depreciation rate
    ρ::Float64 = 0.90                           # persistence of TFP
    σ_ε::Float64 = 0.02                         # std dev of TFP innovation
    N_z::Int = 7                                # number of productivity states
    N_k::Int = 200                              # number of capital grid points


    z_grid::Vector{Float64} = Float64[]         # filled by setup!
    P_z::Matrix{Float64} = zeros(0,0)           # filled by setup!
    k_grid::Vector{Float64} = Float64[]         # filled by setup!
end

RBCModel

• The economic parameters (\(\beta, \alpha, \delta, \rho, \sigma_\varepsilon\)) have sensible defaults
• The grids and transition matrix start empty
  • we fill them with a setup! function

The setup! Function

• Compute everything the model needs from the primitives:

function setup!(m::RBCModel)
    (; α, β, δ, ρ, σ_ε, N_z, N_k) = m
    # Discretize productivity
    log_z, P = rouwenhorst(N_z, ρ, σ_ε)
    m.z_grid = exp.(log_z)
    m.P_z = P
    # Capital grid centered on steady state
    k_ss = (α * β / (1 - β*(1-δ)))^(1/(1-α))
    m.k_grid = collect(LinRange(0.2k_ss, 1.8k_ss, N_k))
    return m
end

setup! (generic function with 1 method)

• setup! mutates the struct (note the ! convention)
  • it fills in z_grid, P_z, and k_grid
• Calling setup! again after changing a parameter rebuilds everything consistently

Creating and Inspecting the Model

m = setup!(RBCModel())                               # create with defaults and set up

RBCModel(0.96, 0.33, 0.08, 0.9, 0.02, 7, 200, [0.8936953824472262, 0.927811339397771, 0.9632296400120643, 1.0, 1.0381740329206182, 1.0778053226306605, 1.118949498498781], [0.7350918906249998 0.23213428125000016 … 1.781250000000008e-6 1.5625000000000085e-8; 0.03868904687500002 0.74527321875 … 2.82187500000001e-5 2.968750000000013e-7; … ; 2.968750000000013e-7 2.8218750000000093e-5 … 0.74527321875 0.03868904687500002; 1.5625000000000085e-8 1.781250000000008e-6 … 0.23213428125000016 0.7350918906249998], [0.8867671427443429, 0.9224160731059244, 0.9580650034675061, 0.9937139338290878, 1.0293628641906694, 1.065011794552251, 1.1006607249138325, 1.1363096552754142, 1.1719585856369958, 1.2076075159985775  …  7.660063911444851, 7.695712841806433, 7.731361772168014, 7.767010702529595, 7.802659632891177, 7.838308563252759, 7.873957493614341, 7.909606423975922, 7.945255354337504, 7.980904284699085])

println("z_grid: ", round.(m.z_grid, digits=4))
println("k_grid: [$(round(m.k_grid[1], digits=3)), ..., $(round(m.k_grid[end], digits=3))]")
println("N_k = $(m.N_k) grid points")

z_grid: [0.8937, 0.9278, 0.9632, 1.0, 1.0382, 1.0778, 1.1189]
k_grid: [0.887, ..., 7.981]
N_k = 200 grid points

The Steady State

• A helper to compute steady-state values from the model parameters:

function steady_state(m::RBCModel)
    (; α, β, δ) = m
    k_ss = (α * β / (1 - β*(1-δ)))^(1/(1-α))
    y_ss = k_ss^α
    c_ss = y_ss - δ * k_ss
    i_ss = δ * k_ss
    r_ss = α * k_ss^(α-1)
    return (k=k_ss, y=y_ss, c=c_ss, i=i_ss, r=r_ss)
end

steady_state (generic function with 1 method)

ss = steady_state(m)
println("k_ss = $(round(ss.k, digits=3)),  y_ss = $(round(ss.y, digits=3)),  c_ss = $(round(ss.c, digits=3)),  r_ss = $(round(ss.r, digits=3))")

k_ss = 4.434,  y_ss = 1.635,  c_ss = 1.28,  r_ss = 0.122

Linear Interpolation

Why We Need Interpolation

The Problem

• In McCall and Rust, the state space was discrete
  • we stored \(V\) as a vector and looked up values by index
• In the RBC model, capital \(k\) is continuous

• We store \(V\) on a grid: \(V(k_1, s), V(k_2, s), \ldots, V(k_N, s)\)
• But when we optimize, the optimal \(k'\) will not land on a grid point!

• We need to evaluate \(V(k', s)\) for any \(k'\), not just grid points
• Solution: interpolation
  • use known grid values to approximate \(V\) at off-grid points

Visual Intuition

• We know \(V\) at the blue dots
  • but we need \(V(2.7)\)!

How Linear Interpolation Works

The Idea

• Given two adjacent grid points \(k_i\) and \(k_{i+1}\) with known values \(V_i\) and \(V_{i+1}\):

\[ V(k') \approx V_i + \frac{k' - k_i}{k_{i+1} - k_i} \cdot (V_{i+1} - V_i) \]

• Using the interpolation weight \(w = \frac{k' - k_i}{k_{i+1} - k_i}\):

\[ V(k') \approx (1 - w) \cdot V_i + w \cdot V_{i+1} \]

Intuition

Linear interpolation draws a straight line between adjacent known points and reads off the value at \(k'\). The weight \(w\) measures how far \(k'\) is from \(k_i\) toward \(k_{i+1}\).

Step by Step

1. Find the interval: locate \(i\) such that \(k_i \leq k' \leq k_{i+1}\)
2. Compute the weight: \(w = \frac{k' - k_i}{k_{i+1} - k_i} \in [0, 1]\)
3. Interpolate: \(\hat{V} = (1-w) V_i + w V_{i+1}\)

• If \(w = 0\): \(k'\) is exactly at \(k_i\) → \(\hat{V} = V_i\)
• If \(w = 1\): \(k'\) is exactly at \(k_{i+1}\) → \(\hat{V} = V_{i+1}\)
• If \(w = 0.5\): \(k'\) is midway → \(\hat{V} = \frac{V_i + V_{i+1}}{2}\)

Interpolations.jl

• Julia’s Interpolations.jl package does all of this for us with clean call syntax:

# Known grid points
k_demo = [1.0, 2.0, 3.0, 4.0, 5.0]
V_demo = [0.0, 0.7, 1.2, 1.5, 1.6]

V_itp = linear_interpolation(k_demo, V_demo)       # build interpolation object

5-element extrapolate(interpolate((::Vector{Float64},), ::Vector{Float64}, Gridded(Linear())), Throw()) with element type Float64:
 0.0
 0.7
 1.2
 1.5
 1.6

• Now V_itp behaves like a function — just call it with any \(k'\):

V_itp(2.0)         # exact grid point → 0.7

0.7

V_itp(2.5)         # midpoint → (0.7 + 1.2)/2 = 0.95

0.95

V_itp(2.7)         # arbitrary point → 0.3×0.7 + 0.7×1.2 = 1.05

1.05

Call Syntax

linear_interpolation(xs, ys) returns a callable object. Writing V_itp(k') is just like calling a function — clean, readable, and exactly what we need inside an optimizer.

Visualizing the Interpolation

k_fine = LinRange(1.0, 5.0, 200)

plot(k_fine, V_itp.(k_fine), linewidth=2, color=:steelblue,
     label="V_itp(k)", legend=:bottomright, size=(700, 400))
scatter!(k_demo, V_demo, markersize=8, color=:red, label="Grid points")
xlabel!("k")
ylabel!("V(k)")
title!("Linear Interpolation in Action")

• The interpolation connects the dots with straight line segments
• Note the broadcasting: V_itp.(k_fine) evaluates the interpolant at every point in the vector

Accuracy vs. Grid Size

• With more grid points, the piecewise-linear approximation better captures curvature:

Rule of Thumb

For smooth economic functions, 100–500 grid points typically give excellent accuracy with linear interpolation.

Solving the RBC Model

The Planning Problem

Recall: The Bellman Equation

\[ V(k, s) = \max_{k'} \left\{ \log(c) + \beta \, \mathbb{E}\big[V(k', s') \mid s\big] \right\}, \qquad c = \bar z_s \, k^\alpha + (1-\delta) k - k' \]

• An infinite-horizon, stochastic, recursive problem — it looks intimidating
• But at its heart is a simple one-period trade-off: consume now vs. save for later

The plan

Boil this down to a problem you can solve in a few lines of Julia, then iterate it to recover \(V\).

A Simple Problem

• Strip away the complexity: a planner has resources \(R\) to split between consumption \(c\) and savings \(k'\):

\[ c + k' = R \]

• It values consumption today, \(\log(c)\), and future capital, \(\beta\, v^f(k')\):

\[ \max_{k'} \;\; \log(\underbrace{R - k'}_{=\,c}) \;+\; \beta\, v^f(k') \]

• \(v^f(k')\) is the (known) value of having capital \(k'\) into the future
  • we’ll have to find this eventually, but for now we’ll take it as given
• Assume that planner must choose \(k'\) inside the range of the grid

This is ALL you need

Learn to solve this two-variable problem we can solve the entire infinite-horizon Bellman equation.

Step 1: Solve It in Julia

Numerical Optimization with Optim.jl

• The problem rarely has a closed form
• so we maximize \(\log(R - k') + \beta\, v^f(k')\) with Optim.jl:

obj(kp) = -(log(resources - kp) + β * vf(kp))     # negative — we minimize
res = optimize(obj, k_min, k_max)

• resources is \(R\); vf is the future value \(v^f\)
• We minimize the negative payoff (Optim minimizes by convention)

Example: \(v^f(k') = 0\)

• The planner places no value on the future, assuming \(R=5\)

#inputs
resources = 5.
vf_zero(kp) = 0.0

#solve for max
obj(kp) = -(log(resources - kp) + m.β * vf_zero(kp))
k_max = min(m.k_grid[end], resources - 1e-8)
k_min = m.k_grid[1]
res = optimize(obj, k_min, k_max, Brent())

(kp = Optim.minimizer(res), c = resources - Optim.minimizer(res))

(kp = 0.8867671586152314, c = 4.113232841384769)

• Compare to

k_min

0.8867671427443429

• With no future value, the planner sets \(k'\) at the lower bound and consumes nearly everything

Example: \(v^f(k') = \log(k')\)

• Now give capital a future value of \(\log(k')\):

#inputs
vf_log(kp) = log(kp)
resources = 5.

#solve for max
obj(kp) = -(log(resources - kp) + m.β * vf_log(kp))   # same code, new future value
k_max = min(m.k_grid[end], resources - 1e-8)
k_min = m.k_grid[1]
res = optimize(obj, k_min, k_max)

(kp = Optim.minimizer(res), c = resources - Optim.minimizer(res))

(kp = 2.448979593459382, c = 2.551020406540618)

• The planner now saves since there is a future value that rewards capital
  • The optimum is interior (away from the bounds)

Step 2: Wrap It in a Function

optimal_policy

• Both examples ran the same code: we can wrap it in a function
  • taking resources \(R\) and the future value vf as inputs

function optimal_policy(m, R, vf)
    (; β, k_grid) = m
    obj(kp) = -(log(R - kp) + β * vf(kp))

    k_max = min(k_grid[end], R - 1e-8) #avoid negative consumption
    k_min = k_grid[1]
    res = optimize(obj, k_min, k_max)

    return (V  = -Optim.minimum(res),
            kp = Optim.minimizer(res),
            c  = R - Optim.minimizer(res))
end

optimal_policy (generic function with 1 method)

• Same body as the examples, now taking resources \(R\) and the future value vf directly
• Returns a named tuple (V, kp, c)
  • everything we need: values and choices

Step 3: Iterate the Bellman Equation

The Bellman Equation is the Simple Problem

• The Bellman equation at state (k, s) \[ V(k, s) = \max_{k'} \left\{ \log(c) + \beta \, \mathbb{E}\big[V(k', s') \mid s\big] \right\}, \qquad c = \bar z_s \, k^\alpha + (1-\delta) k - k' \] is exactly our simple problem with two substitutions:

\[ R = \underbrace{\bar z_s \, k^\alpha + (1-\delta)k}_{\text{resources at }(k,\,s)}, \qquad v^f(k') = \underbrace{\mathbb{E}\big[V(k', s') \mid s\big]}_{\text{expected continuation value}} \]

• So optimal_policy already solves the Bellman equation at one state if we can supply \(v^f\)

• But \(v^f\) depends on \(V\), which is what we’re solving for!

• We need to find the fixed point: \(V\) that we put into the bellman equation is the V we get out.

The iteration loop

Guess \(V\) \(\rightarrow\) build \(v^f = \mathbb{E}[V \mid s]\) \(\rightarrow\) run optimal_policy at every \((k,s)\) \(\rightarrow\) get \(V_{\text{new}}\) \(\rightarrow\) repeat until convergence.

From \(V\) to \(v^f\)

• optimal_policy needs \(v^f(\cdot)\) as a callable function of \(k'\)

• But \(V\) is stored on a grid: \(V[i, s'] = V(k_i, s')\)

• Compute the expected continuation value on the grid (rows = \(k\), columns = current state \(s\)):

\[ \mathbb{E}V[i', s] \;=\; \sum_{s'=1}^{S} P_{s, s'} \, V[i', s'] \]

• Interpolate each column \(\mathbb{E}V[:, s]\) into the callable \(v^f(\cdot, s)\)

compute_vf: The Future Value for State \(s\)

• compute_vf(m, V, s) builds the future value \(v^f(\cdot)\) for a single current state \(s\):

function compute_vf(m, V, s)
    (; k_grid, P_z) = m
    Nk = length(k_grid)
    S = length(P_z[s, :])
    EV_s = zeros(Nk)                        # E[V(k', s') | s] on the grid
    for i′ in 1:Nk
        EV_s[i′] = sum(P_z[s, s′] * V[i′, s′] for s′ in 1:S)
    end
    #construct the interpolated callable function
    return linear_interpolation(k_grid, EV_s)
end

compute_vf (generic function with 1 method)

• Interpolating it returns the callable \(v^f(k')\)
  • exactly what optimal_policy expects as input

bellman_operator: Pseudocode

function bellman_operator(m::RBCModel, V)
    # For each state s:
    #   vf = compute_vf(m, V, s)            # expected continuation value
    #   For each capital k_i:
    #     R = resources at (k_i, s)
    #     optimal_policy(m, R, vf)          # optimize
    # Return V_new, kp_policy, c_policy
end

bellman_operator: Full Function

function bellman_operator(m::RBCModel, V)
    (; N_k, N_z, α, δ, z_grid, k_grid) = m
    V_new     = similar(V)
    kp_policy = similar(V)
    c_policy  = similar(V)
    for s in 1:N_z
        vf = compute_vf(m, V, s)                     # E[V(k', s')|s] as a function
        z  = z_grid[s]
        for ik in 1:N_k
            k = k_grid[ik]
            R = z * k^α + (1 - δ) * k                # resources at (k, s)
            result = optimal_policy(m, R, vf)
            V_new[ik, s]     = result.V
            kp_policy[ik, s] = result.kp
            c_policy[ik, s]  = result.c
        end
    end
    return V_new, kp_policy, c_policy
end

bellman_operator (generic function with 1 method)

• For each state \(s\), compute_vf builds \(v^f\) once, then we reuse it across all \(k\)
• For each \(k\), we form resources \(R\) and hand \((R, v^f)\) to optimal_policy

Iterating to Convergence

The VFI Algorithm

1. Start with an initial guess \(V^0(k, s) = 0\)
2. Apply bellman_operator to get \(V^{n+1}\) and policies
3. Repeat until \(\|V^{n+1} - V^n\|_\infty < \text{tol}\)

Convergence

The Bellman operator is a contraction mapping with modulus \(\beta\). By the Banach fixed-point theorem, iterating converges to the unique fixed point \(V^*\) at a geometric rate.

solveRBC: Pseudocode

function solveRBC(m::RBCModel; tol=1e-6, maxiter=1000, verbose=true)
    # Initial guess V = 0
    # Iterate Bellman operator until convergence
    # Build interpolated policy functions
    # Return solution (grids + interpolated functions)
end

solveRBC: Full Function

function solveRBC(m::RBCModel; tol=1e-6, maxiter=1000, verbose=true)
    (; k_grid, N_k, N_z) = m
    V = zeros(N_k, N_z)                              # initial guess: V = 0
    V_new, kp_policy, c_policy = bellman_operator(m, V)
    dist = 1.0; iter = 0
    while dist > tol && iter < maxiter
        V_new, kp_policy, c_policy = bellman_operator(m, V)
        dist = maximum(abs.(V_new - V))
        V = V_new
        iter += 1
    end
    if dist < tol && verbose
        println("Converged in $iter iterations (dist = $(round(dist,digits=8)))")
    end
    kp_itp = [linear_interpolation(k_grid, kp_policy[:, s]) for s in 1:N_z]
    c_itp  = [linear_interpolation(k_grid, c_policy[:, s])  for s in 1:N_z]
    return (V=V, kp=kp_policy, c=c_policy, kp_itp=kp_itp, c_itp=c_itp)
end

solveRBC (generic function with 1 method)

Running the Solver

m = setup!(RBCModel())
sol = solveRBC(m)

Converged in 294 iterations (dist = 9.6e-7)

(V = [0.5450767096792496 0.9524145298017854 … 2.6015946443534506 3.019041640967933; 0.6288520459016571 1.0350077981943682 … 2.679696895724538 3.096067287082053; … ; 7.561499679633558 7.877213571940097 … 9.168728441903232 9.498811277963839; 7.582033994564345 7.897513180113208 … 9.188095953883254 9.517944455040283], kp = [1.1539249754885443 1.1719585978272085 … 1.2641278395360356 1.2854963280075271; 1.1875771745434114 1.2076075203881291 … 1.2998968026169335 1.321739180753339; … ; 7.358439432055147 7.402519069664838 … 7.589068507785396 7.644303445005629; 7.388480155611816 7.432574853973526 … 7.6244149648272 7.674711618118721], c = [0.5208484138806513 0.5356042861903008 … 0.5875971126193718 0.6057730735793494; 0.5312381980315675 0.5444266125860389 … 0.5981867053624921 0.6164064753278229; … ; 1.7222221900487993 1.7457497838398366 … 1.8564418373720857 1.8827418787174723; 1.7275968224640046 1.7512093081018767 … 1.8570501398142119 1.8884090088945928], kp_itp = Interpolations.Extrapolation{Float64, 1, Interpolations.GriddedInterpolation{Float64, 1, Vector{Float64}, Gridded{Linear{Throw{OnGrid}}}, Tuple{Vector{Float64}}}, Gridded{Linear{Throw{OnGrid}}}, Throw{Nothing}}[200-element extrapolate(interpolate((::Vector{Float64},), ::Vector{Float64}, Gridded(Linear())), Throw()) with element type Float64:
 1.1539249754885443
 1.1875771745434114
 1.2210546041585544
 1.254394521373463
 1.2877132936648057
 1.3209209761025729
 1.3539903647998963
 1.38687416806246
 1.4215010942520319
 1.4571500386302887
 ⋮
 7.143981576389601
 7.1740246277605815
 7.20407646735
 7.2341289457713644
 7.2679255897914095
 7.299585013582776
 7.32875276669357
 7.358439432055147
 7.388480155611816, 200-element extrapolate(interpolate((::Vector{Float64},), ::Vector{Float64}, Gridded(Linear())), Throw()) with element type Float64:
 1.1719585978272085
 1.2076075203881291
 1.2432564423329713
 1.2789053815738296
 1.314554293340328
 1.3485214664181069
 1.3818308248663593
 1.4149994120024603
 1.448032677319324
 1.4809721818651598
 ⋮
 7.186995182999644
 7.217099560710081
 7.2471965899186115
 7.277251820006923
 7.307196745588469
 7.3392235511660875
 7.372599333029496
 7.402519069664838
 7.432574853973526, 200-element extrapolate(interpolate((::Vector{Float64},), ::Vector{Float64}, Gridded(Linear())), Throw()) with element type Float64:
 1.1966980794855495
 1.2312146967315636
 1.265516107828779
 1.2996393305904002
 1.3335755409516568
 1.3673011903049008
 1.4008665058087282
 1.4342978264030621
 1.467710287221131
 1.501013499357431
 ⋮
 7.232276638015115
 7.262706359280592
 7.292906507588428
 7.323010212386788
 7.353145228765818
 7.383271854058047
 7.412869916999862
 7.446170343382582
 7.478401314881104, 200-element extrapolate(interpolate((::Vector{Float64},), ::Vector{Float64}, Gridded(Linear())), Throw()) with element type Float64:
 1.214438507753562
 1.249329182933468
 1.284060793261448
 1.318617342007283
 1.3529948319138827
 1.387173814553634
 1.4215010995059496
 1.4571500446129377
 1.4927989689595673
 1.5284478884385593
 ⋮
 7.27483353929229
 7.30486847092937
 7.33922350881025
 7.3706292900679875
 7.400670341696727
 7.430708359222595
 7.460874916235503
 7.4907740914615495
 7.520998813424199, 200-element extrapolate(interpolate((::Vector{Float64},), ::Vector{Float64}, Gridded(Linear())), Throw()) with element type Float64:
 1.243256439108657
 1.2789053611810666
 1.314326504384777
 1.349168110882513
 1.383838030878106
 1.4183457271039948
 1.452683227649133
 1.4868561944267642
 1.5208988390855722
 1.5548233755828587
 ⋮
 7.32453710656391
 7.354351547505431
 7.384474358280346
 7.414727775174974
 7.4461704150363035
 7.481012310954662
 7.511367300478453
 7.541802961686753
 7.572105934333229, 200-element extrapolate(interpolate((::Vector{Float64},), ::Vector{Float64}, Gridded(Linear())), Throw()) with element type Float64:
 1.2641278395360356
 1.2998968026169335
 1.335421971809969
 1.3707277829524782
 1.4058371175140867
 1.4407628925214644
 1.475511036917388
 1.5100448785023561
 1.5444170894210265
 1.5786451316609817
 ⋮
 7.3748724256684515
 7.407212675621696
 7.436888151196912
 7.467277755805652
 7.497793765074748
 7.52829386891891
 7.558739445264765
 7.589068507785396
 7.6244149648272, 200-element extrapolate(interpolate((::Vector{Float64},), ::Vector{Float64}, Gridded(Linear())), Throw()) with element type Float64:
 1.2854963280075271
 1.321739180753339
 1.35770883878131
 1.3934440944236473
 1.4289641715571866
 1.4642853158249232
 1.4994149733272448
 1.5344172925630561
 1.5692669026857295
 1.6039408120276473
 ⋮
 7.4258341984430976
 7.456311363010271
 7.486785276776135
 7.517468307491838
 7.5529894065043655
 7.583250020060552
 7.613775363219099
 7.644303445005629
 7.674711618118721], c_itp = Interpolations.Extrapolation{Float64, 1, Interpolations.GriddedInterpolation{Float64, 1, Vector{Float64}, Gridded{Linear{Throw{OnGrid}}}, Tuple{Vector{Float64}}}, Gridded{Linear{Throw{OnGrid}}}, Throw{Nothing}}[200-element extrapolate(interpolate((::Vector{Float64},), ::Vector{Float64}, Gridded(Linear())), Throw()) with element type Float64:
 0.5208484138806513
 0.5312381980315675
 0.5415152206899543
 0.5516598791682972
 0.5615717885070037
 0.5713554145071493
 0.5810511979532422
 0.590718534532372
 0.598439816122051
 0.6049463549252909
 ⋮
 1.6885484496770946
 1.6939772091502023
 1.699388929773538
 1.704791824077538
 1.706442427670308
 1.710222187997756
 1.7164856163865876
 1.7222221900487993
 1.7275968224640046, 200-element extrapolate(interpolate((::Vector{Float64},), ::Vector{Float64}, Gridded(Linear())), Throw()) with element type Float64:
 0.5356042861903008
 0.5444266125860389
 0.5529504324503285
 0.5611940561608639
 0.5691741211222301
 0.578587413740062
 0.5884237521585731
 0.5981786591113816
 0.6078582090857514
 0.6174314381582562
 ⋮
 1.7124338513133859
 1.7179033921560416
 1.7233717159093986
 1.7288733387461717
 1.7344768308088456
 1.7379900715653633
 1.7401460279246512
 1.7457497838398366
 1.7512093081018767, 200-element extrapolate(interpolate((::Vector{Float64},), ::Vector{Float64}, Gridded(Linear())), Throw()) with element type Float64:
 0.5449060064284283
 0.5553062906948543
 0.5656118787408753
 0.5758047807078135
 0.5859110467024122
 0.5959698760044441
 0.605945308254044
 0.6158240502591259
 0.6255029307834394
 0.6350833403922604
 ⋮
 1.7366052094874025
 1.7418554123470686
 1.7473262955900095
 1.7528847978249376
 1.758403231189643
 1.763921364764399
 1.769959435424009
 1.7722865821999845
 1.7756746874648446, 200-element extrapolate(interpolate((::Vector{Float64},), ::Vector{Float64}, Gridded(Linear())), Throw()) with element type Float64:
 0.5625062700187375
 0.5729951612623883
 0.5833213847655176
 0.5935206915866327
 0.6036149626122915
 0.6136398943915031
 0.6232634887648316
 0.631325931646812
 0.6391613245596512
 0.6467810732481227
 ⋮
 1.7661524153768404
 1.7719074600971494
 1.7733331665338659
 1.777698968208818
 1.7834204079221534
 1.7891358590957944
 1.7947138162556646
 1.800550267973044
 1.8060523522168337, 200-element extrapolate(interpolate((::Vector{Float64},), ::Vector{Float64}, Gridded(Linear())), Throw()) with element type Float64:
 0.5703781272564015
 0.5805890983042137
 0.5906938337977019
 0.6010646000048081
 0.6113121130414281
 0.6214437785528406
 0.6314829452817026
 0.6414380120912255
 0.6512876620845687
 0.6610315389518908
 ⋮
 1.7913054598325227
 1.7973952486181766
 1.8031670831064508
 1.8087988003630944
 1.813231855999745
 1.8142562885061464
 1.819758331049984
 1.8251704754185765
 1.8307061508889015, 200-element extrapolate(interpolate((::Vector{Float64},), ::Vector{Float64}, Gridded(Linear())), Throw()) with element type Float64:
 0.5875971126193718
 0.5981867053624921
 0.6086733268928708
 0.6190538326933466
 0.6293246163653212
 0.639490273193849
 0.6495608380634386
 0.6595875807065921
 0.6695312159751745
 0.679386592183288
 ⋮
 1.818684331215687
 1.8223669258712993
 1.828704344563314
 1.834317760032386
 1.8397949718763495
 1.8452783644930637
 1.850806633371512
 1.8564418373720857
 1.8570501398142119, 200-element extrapolate(interpolate((::Vector{Float64},), ::Vector{Float64}, Gridded(Linear())), Throw()) with element type Float64:
 0.6057730735793494
 0.6164064753278229
 0.6269530692711671
 0.6373961671729069
 0.6477365560344299
 0.6579761710400291
 0.668124139320271
 0.678131467155956
 0.6880374234258815
 0.6978777800332117
 ⋮
 1.8484034129948244
 1.8540722359265898
 1.8597339794733694
 1.8651763549432578
 1.865770489085567
 1.8716150128028435
 1.8771847872545626
 1.8827418787174723
 1.8884090088945928])

Patience

This takes a moment because each Bellman iteration solves \(N_k \times N_z\) optimization problems. More efficient methods exist (policy iteration, endogenous grid method), but VFI with interpolation is the clearest to understand.

Calibration

Why Calibrate?

• We need values for five parameters: \((\alpha, \beta, \delta, \rho, \sigma_\varepsilon)\)
• Rather than estimating them econometrically, the RBC literature calibrates: pick values so that the model’s steady state matches long-run averages in the data

Calibration vs. Estimation

Calibration uses different data moments than the ones we’re trying to explain. We match long-run averages (growth accounting, interest rates) and then ask: can the model reproduce business cycle fluctuations?

Matching Long-Run Moments

Parameter	Target	Data Source	Value
\(\alpha\)	Capital’s share of income	National accounts: \(rK/Y \approx 1/3\)	0.33
\(\beta\)	Real interest rate	\(r = 1/\beta - 1\approx 4\%\)/yr	0.96
\(\delta\)	Investment–capital ratio	\(I/K \approx 8\%\)/yr	0.08

• \(\alpha = 0.33\): from NIPA data, labor’s share of income \(\approx 2/3\), so capital’s share \(\approx 1/3\)
• \(\beta = 0.96\): steady-state Euler equation \(r_{ss} = 1/\beta - 1\); with \(r \approx 4\%\) we get \(\beta \approx 0.96\)
• \(\delta = 0.08\): annual depreciation of the capital stock is roughly 8%

The Productivity Process

• KP model the Solow residual (total factor productivity) as an AR(1):

\[ \log z_t = \rho \, \log z_{t-1} + \varepsilon_t, \qquad \varepsilon_t \sim N(0, \sigma_\varepsilon^2) \]

• The Solow residual is computed from data:

\[ \log z_t = \log Y_t - \alpha \log K_t - (1-\alpha) \log L_t \]

• Estimate the AR(1) by OLS regression of \(\log z_t\) on \(\log z_{t-1}\):
  • \(\hat{\rho} \approx 0.90\) — productivity shocks are very persistent
  • \(\hat{\sigma}_\varepsilon \approx 0.02\) — innovations are small (≈ 2% of output)

Our Calibration

println("Calibrated Parameters:")
println("  α  = $(m.α)   (capital share)")
println("  β  = $(m.β)   (discount factor → r ≈ $(round(1/m.β - 1 + m.δ, digits=3)))")
println("  δ  = $(m.δ)   (depreciation rate)")
println("  ρ  = $(m.ρ)   (TFP persistence)")
println("  σ_ε = $(m.σ_ε)  (TFP innovation std dev)")

Calibrated Parameters:
  α  = 0.33   (capital share)
  β  = 0.96   (discount factor → r ≈ 0.122)
  δ  = 0.08   (depreciation rate)
  ρ  = 0.9   (TFP persistence)
  σ_ε = 0.02  (TFP innovation std dev)

println("Implied Steady State:")
println("  K/Y  = $(round(ss.k / ss.y, digits=2))    (data: ≈ 2.8)")
println("  I/Y  = $(round(ss.i / ss.y, digits=2))    (data: ≈ 0.25-0.30)")
println("  C/Y  = $(round(ss.c / ss.y, digits=2))    (data: ≈ 0.65-0.70)")

Implied Steady State:
  K/Y  = 2.71    (data: ≈ 2.8)
  I/Y  = 0.22    (data: ≈ 0.25-0.30)
  C/Y  = 0.78    (data: ≈ 0.65-0.70)

The Calibration Strategy

1. Match \(\alpha\) from factor shares in national income accounts
2. Match \(\beta\) and \(\delta\) from long-run interest rates and investment rates
3. Match \(\rho\) and \(\sigma_\varepsilon\) from the estimated Solow residual process
4. Then ask: does the model reproduce business cycle statistics it was not calibrated to match?

Results

The Value Function

plot(legend=:bottomright, size=(700, 400))
for s in [1, 4, 7]
    plot!(m.k_grid, sol.V[:, s], linewidth=2,
          label="z = $(round(m.z_grid[s], digits=3))")
end
xlabel!("Capital k")
ylabel!("V(k, s)")
title!("Value Function")
vline!([ss.k], color=:gray, linestyle=:dash, label="k_ss")

• Higher productivity \(z\) → higher value (more output at every \(k\))
• Value is concave in \(k\) (diminishing returns to capital)

The Savings Policy \(k'(k, s)\)

plot(legend=:topleft, size=(700, 400))
for s in [1, 4, 7]
    plot!(m.k_grid, sol.kp[:, s], linewidth=2,
          label="z = $(round(m.z_grid[s], digits=3))")
end
plot!(m.k_grid, m.k_grid, color=:gray, linestyle=:dash, label="45° line")
xlabel!("Current Capital k")
ylabel!("Next Period Capital k'")
title!("Savings Policy Function")
vline!([ss.k], color=:gray, linestyle=:dot, label="")

• Above the 45° line: capital is growing (\(k' > k\))
• Below the 45° line: capital is shrinking (\(k' < k\))
• The crossing point is the steady state for each \(z\)

The Consumption Policy \(c(k, s)\)

plot(legend=:topleft, size=(700, 400))
for s in [1, 4, 7]
    plot!(m.k_grid, sol.c[:, s], linewidth=2,
          label="z = $(round(m.z_grid[s], digits=3))")
end
xlabel!("Capital k")
ylabel!("Consumption c(k, s)")
title!("Consumption Policy Function")
vline!([ss.k], color=:gray, linestyle=:dash, label="k_ss")

• Consumption is increasing in both \(k\) and \(z\)
• High productivity → more output → consume more and save more

The Investment Policy

i_policy = sol.kp .- (1 - m.δ) .* m.k_grid         # i = k' - (1-δ)k

plot(legend=:topleft, size=(700, 400))
for s in [1, 4, 7]
    plot!(m.k_grid, i_policy[:, s], linewidth=2,
          label="z = $(round(m.z_grid[s], digits=3))")
end
xlabel!("Capital k")
ylabel!("Investment i(k, s)")
title!("Investment Policy Function")
hline!([m.δ * ss.k], color=:gray, linestyle=:dash, label="Steady state i")

• Investment responds strongly to productivity
  • when \(z\) is high, the agent invests heavily
• This foreshadows a key business cycle fact: investment is volatile

Simulation

Simulating the Economy

The Simulation Function

function simulateRBC(m::RBCModel, sol, T; seed=42)
    # Initialize at steady state
    # For each period: 
        # use policies to get k', c, i
        # draw next period's productivity from the Markov chain
        # update state k_{t+1} = k'
    # Return a DataFrame with columns t, k, z, y, c, i
end

The Simulation Function

function simulateRBC(m::RBCModel, sol, T; seed=42)
    Random.seed!(seed)
    (; α, δ, z_grid, P_z, N_z) = m
    k_sim = zeros(T); z_sim = zeros(T); s_sim = zeros(Int, T)
    y_sim = zeros(T); c_sim = zeros(T); i_sim = zeros(T)
    k_sim[1] = steady_state(m).k                    # start at steady state
    s_sim[1] = (N_z + 1) ÷ 2                        # start at median z
    for t in 1:T
        s = s_sim[t]
        z = z_grid[s]
        kp = sol.kp_itp[s](k_sim[t])               # interpolated policy!
        z_sim[t] = z
        y_sim[t] = z * k_sim[t]^α
        c_sim[t] = y_sim[t] + (1-δ)*k_sim[t] - kp
        i_sim[t] = kp - (1-δ)*k_sim[t]
        if t < T
            k_sim[t+1] = kp
            s_sim[t+1] = rand(Categorical(P_z[s, :]))
        end
    end
    return DataFrame(t=1:T, k=k_sim, z=z_sim, y=y_sim, c=c_sim, i=i_sim)
end

simulateRBC (generic function with 1 method)

Using Interpolated Policies

• Since solveRBC returns interpolated functions, simulation just calls them:

kp = sol.kp_itp[s](k_sim[t])         # reads like math: k'(k_t, s_t)

• The current capital \(k_t\) may not be on our grid (because last period’s optimal \(k'\) was found by continuous optimization)
• The interpolated policy function evaluates \(k'(k_t, s_t)\) at any \(k_t\)

Interpolation Everywhere

Linear interpolation is used in two places: (1) inside VFI to evaluate the continuation value, and (2) in simulation to evaluate the policy function. Both are essential for handling the continuous state space.

Running the Simulation

sim = simulateRBC(m, sol, 200)
first(sim, 6)

6×6 DataFrame

Row	t	k	z	y	c	i
	Int64	Float64	Float64	Float64	Float64	Float64
1	1	4.43384	1.0	1.6347	1.27999	0.354707
2	2	4.43384	1.0	1.6347	1.27999	0.354707
3	3	4.43384	1.0	1.6347	1.27999	0.354707
4	4	4.43384	1.0	1.6347	1.27999	0.354707
5	5	4.43384	1.0	1.6347	1.27999	0.354707
6	6	4.43384	1.0	1.6347	1.27999	0.354707

• Each row is one period; columns are the macro aggregates
• DataFrames give us a clean, tabular view of the simulation and powerful tools for computing statistics

Output, Consumption, and Investment

p1 = plot(sim.y, linewidth=2, label="", title="Output", size=(800, 150))
p2 = plot(sim.c, linewidth=2, label="", title="Consumption", size=(800, 150))
p3 = plot(sim.i, linewidth=2, label="", title="Investment", size=(800, 150))
plot(p1, p2, p3, layout=(3, 1), size=(800, 500), xlabel="Period")

• Consumption is smoother than output
  • households smooth consumption over time
• Investment is more volatile than output
  • the main adjustment margin

Productivity and Capital

p1 = plot(sim.z, linewidth=2, color=:steelblue, label="Productivity z",
          title="Technology Shock", size=(800, 200))
p2 = plot(sim.k, linewidth=2, color=:darkgreen, label="Capital k",
          title="Capital Stock", size=(800, 200))
hline!(p2, [ss.k], color=:gray, linestyle=:dash, label="k_ss")
plot(p1, p2, layout=(2,1), size=(800, 450))

• Capital responds sluggishly to productivity shocks
  • it takes time to build (or deplete) capital
• This is the propagation mechanism that turns short-lived shocks into persistent cycles

Business Cycle Statistics

Computing Moments

• Simulate a long series and construct log-deviation :

sim_long = simulateRBC(m, sol, 10_000)

sim_long.y_dev = log.(sim_long.y) .- mean(log.(sim_long.y))
sim_long.c_dev = log.(sim_long.c) .- mean(log.(sim_long.c))
sim_long.i_dev = log.(sim_long.i) .- mean(log.(sim_long.i))
sim_long.k_dev = log.(sim_long.k) .- mean(log.(sim_long.k))

first(sim_long, 6) #report the first 6 rows

6×10 DataFrame

Row	t	k	z	y	c	i	y_dev	c_dev	i_dev	k_dev
	Int64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64
1	1	4.43384	1.0	1.6347	1.27999	0.354707	-0.00359466	-0.00367429	-0.000783353	-0.00496687
2	2	4.43384	1.0	1.6347	1.27999	0.354707	-0.00359466	-0.00367429	-0.000783354	-0.00496687
3	3	4.43384	1.0	1.6347	1.27999	0.354707	-0.00359466	-0.00367429	-0.000783355	-0.00496687
4	4	4.43384	1.0	1.6347	1.27999	0.354707	-0.00359466	-0.00367429	-0.000783356	-0.00496688
5	5	4.43384	1.0	1.6347	1.27999	0.354707	-0.00359466	-0.00367429	-0.000783357	-0.00496688
6	6	4.43384	1.0	1.6347	1.27999	0.354707	-0.00359466	-0.00367429	-0.000783359	-0.00496688

Standard Deviations

stats = DataFrame(
    Variable   = ["Output", "Consumption", "Investment", "Capital"],
    σ_pct      = [100*std(sim_long.y_dev), 100*std(sim_long.c_dev), 100*std(sim_long.i_dev), 100*std(sim_long.k_dev)],
    σ_relative = [1.0, std(sim_long.c_dev)/std(sim_long.y_dev), std(sim_long.i_dev)/std(sim_long.y_dev), std(sim_long.k_dev)/std(sim_long.y_dev)],
    corr_y     = [1.0, cor(sim_long.y_dev, sim_long.c_dev), cor(sim_long.y_dev, sim_long.i_dev), cor(sim_long.y_dev, sim_long.k_dev)]
)
stats

4×4 DataFrame

Row	Variable	σ_pct	σ_relative	corr_y
	String	Float64	Float64	Float64
1	Output	6.15173	1.0	1.0
2	Consumption	5.29062	0.860021	0.965192
3	Investment	11.1569	1.81362	0.892204
4	Capital	6.4229	1.04408	0.853609

Business Cycle Observations

Key Business Cycle Facts

The model reproduces two central stylized facts:

1. Consumption is less volatile than output (\(\sigma_c / \sigma_y < 1\)) — consumption smoothing
2. Investment is more volatile than output — investment absorbs the shocks
3. All variables are procyclical (positively correlated with output)

The RBC Achievement

With just one source of randomness (technology shocks), the model generates realistic comovement among macro aggregates. This was revolutionary in the early 1980s.

Comparative Statics

Effect of Persistence (\(\rho\))

plot(legend=:topright, size=(700, 400), title="Output: Effect of Shock Persistence")
for ρ_val in [0.5, 0.75, 0.95]
    m_trial = setup!(RBCModel(ρ=ρ_val, N_k=100))
    sol_trial = solveRBC(m_trial, tol=1e-5, verbose=false)
    sim_trial = simulateRBC(m_trial, sol_trial, 150)
    plot!(sim_trial.y, linewidth=2, label="ρ = $ρ_val")
end
xlabel!("Period")
ylabel!("Output")

• Higher persistence → longer, more pronounced business cycles
• Low persistence → output fluctuates rapidly but reverts quickly

Effect of Capital Share (\(\alpha\))

plot(legend=:topleft, size=(700, 400), title="Savings Policy: Effect of α")
for α_val in [0.25, 0.33, 0.40]
    m_trial = setup!(RBCModel(α=α_val, N_k=100))
    sol_trial = solveRBC(m_trial, tol=1e-5, verbose=false)
    s_mid = (m_trial.N_z + 1) ÷ 2
    plot!(m_trial.k_grid, sol_trial.kp[:, s_mid], linewidth=2, label="α = $α_val")
end
xlabel!("Current Capital k")
ylabel!("Next Period Capital k'")

• Higher \(\alpha\) → capital is more productive → agent saves more → higher steady-state capital

Effect of Discount Factor (\(\beta\))

plot(legend=:topleft, size=(700, 400), title="Savings Policy: Effect of β")
for β_val in [0.90, 0.96, 0.99]
    m_trial = setup!(RBCModel(β=β_val, N_k=100))
    sol_trial = solveRBC(m_trial, tol=1e-5, verbose=false)
    s_mid = (m_trial.N_z + 1) ÷ 2
    plot!(m_trial.k_grid, sol_trial.kp[:, s_mid], linewidth=2, label="β = $β_val")
end
xlabel!("Current Capital k")
ylabel!("Next Period Capital k'")

• More patient agents (\(\beta\) closer to 1) save more → higher steady-state capital
• Impatient agents consume more today, accumulate less capital

Linearized Solution

The Euler Equation

Notation

• Define deviations from steady state:

\[ \hat{k}_t \equiv k_t - k_{ss}, \qquad \hat{c}_t \equiv c_t - c_{ss}, \qquad \hat{z}_t \equiv z_t - 1 \]

• We’ll linearize by taking a first-order Taylor expansion around the steady state, dropping terms like \(\hat{k}_t \cdot \hat{z}_t\) (second-order small)

Why Linearize?

VFI is accurate but slow (\(N_k \times N_z\) optimization problems per iteration). Linearization gives a policy rule that’s nearly exact near the steady state — and almost instant to compute.

Linearizing the Resource Constraint

The Resource Constraint

• Recall the aggregate resource constraint:

\[ c_t + k_{t+1} = z_t \, k_t^\alpha + (1-\delta) \, k_t \]

• Write each variable as steady state + deviation:

\[ (c_{ss} + \hat{c}_t) + (k_{ss} + \hat{k}_{t+1}) = (1 + \hat{z}_t)(k_{ss} + \hat{k}_t)^\alpha + (1-\delta)(k_{ss} + \hat{k}_t) \]

Differentiating the Resource Constraint

• First-order approximation of \((k_{ss} + \hat{k}_t)^\alpha\):

\[ (k_{ss} + \hat{k}_t)^\alpha \approx k_{ss}^\alpha + \alpha \, k_{ss}^{\alpha-1} \, \hat{k}_t \]

• Expanding the right side (dropping \(\hat{z}_t \cdot \hat{k}_t\)):

\[ \text{RHS} \approx k_{ss}^\alpha + \underbrace{\alpha k_{ss}^{\alpha-1}}_{r_{ss}} \hat{k}_t + \underbrace{k_{ss}^\alpha}_{y_{ss}} \hat{z}_t + (1-\delta) k_{ss} + (1-\delta) \hat{k}_t \]

Linearized Resource Constraint

• Subtracting the steady-state equation (\(c_{ss} + k_{ss} = y_{ss} + (1-\delta) k_{ss}\)), the constant terms cancel:

\[ \hat{c}_t + \hat{k}_{t+1} = \underbrace{(r_{ss} + 1 - \delta)}_{= 1/\beta} \hat{k}_t + y_{ss} \, \hat{z}_t \]

\[ \boxed{\hat{c}_t = \frac{1}{\beta} \, \hat{k}_t - \hat{k}_{t+1} + y_{ss} \, \hat{z}_t} \]

• This tells us \(\hat{c}_t\) as a linear function of \(\hat{k}_t\), \(\hat{k}_{t+1}\), and \(\hat{z}_t\)
• Used \(r_{ss} + 1 - \delta = 1/\beta\) (the steady-state Euler equation)

Linearizing the Euler Equation

First-Order Condition

• Taking the FOC of the Bellman equation with respect to \(k'\):

\[ \frac{1}{c_t} = \beta \, \mathbb{E}_t \left[ \frac{\alpha z_{t+1} k_{t+1}^{\alpha-1} + 1 - \delta}{c_{t+1}} \right] \]

• Left side: marginal cost of saving one more unit (foregone consumption today)
• Right side: marginal benefit (discounted expected return, valued in consumption units)

The Euler Equation

At the optimum, the household is indifferent between consuming one more unit today vs. saving it and consuming the return tomorrow.

Differentiating the Euler Equation

• The Euler equation is \(\displaystyle \frac{1}{c_t} = \beta \, \mathbb{E}_t \left[ \frac{R_{t+1}}{c_{t+1}} \right]\) where \(R_{t+1} = \alpha z_{t+1} k_{t+1}^{\alpha-1} + 1 - \delta\)

• Linearize left side: \(\;\displaystyle \frac{1}{c_{ss} + \hat{c}_t} \approx \frac{1}{c_{ss}} - \frac{\hat{c}_t}{c_{ss}^2}\)

• Linearize right side: expand \(R_{t+1}/c_{t+1}\) around steady state

\[ \frac{R_{ss} + \hat{R}_{t+1}}{c_{ss} + \hat{c}_{t+1}} \approx \frac{R_{ss}}{c_{ss}} + \frac{\hat{R}_{t+1}}{c_{ss}} - \frac{R_{ss} \, \hat{c}_{t+1}}{c_{ss}^2} \]

Linearized Return

• The deviation of the gross return is:

\[ \hat{R}_{t+1} = \underbrace{\frac{(\alpha - 1) \, r_{ss}}{k_{ss}}}_{\text{diminishing MPK}} \hat{k}_{t+1} + \underbrace{r_{ss}}_{\text{productivity effect}} \hat{z}_{t+1} \]

• More capital \(\rightarrow\) lower return (diminishing returns, since \(\alpha - 1 < 0\))
• Higher productivity \(\rightarrow\) higher return

Linearized Euler Equation

• Setting LHS = \(\beta \cdot\) RHS and using \(\beta R_{ss} = 1\) to cancel constant terms:

\[ -\frac{\hat{c}_t}{c_{ss}^2} = \frac{\beta}{c_{ss}} \mathbb{E}_t[\hat{R}_{t+1}] - \frac{1}{c_{ss}^2} \mathbb{E}_t[\hat{c}_{t+1}] \]

• Multiply by \(-c_{ss}^2\):

\[ \boxed{\hat{c}_t = \mathbb{E}_t[\hat{c}_{t+1}] - \beta \, c_{ss} \, \mathbb{E}_t[\hat{R}_{t+1}]} \]

• This is an intertemporal substitution equation: when expected returns are high (\(\hat{R}_{t+1} > 0\)), the household saves more today (\(\hat{c}_t < \mathbb{E}_t[\hat{c}_{t+1}]\))

Solving the Linearized System

Consumption Is Linear in \(\hat{k}_t\) and \(\hat{z}_t\)

• Guess a linear policy rule: \(\;\hat{k}_{t+1} = \phi_k \, \hat{k}_t + \phi_z \, \hat{z}_t\)

• Substitute into the linearized resource constraint \(\hat{c}_t = \tfrac{1}{\beta}\hat{k}_t - \hat{k}_{t+1} + y_{ss}\hat{z}_t\):

\[ \hat{c}_t = \underbrace{\left(\tfrac{1}{\beta} - \phi_k\right)}_{\equiv\; \psi_k} \hat{k}_t + \underbrace{(y_{ss} - \phi_z)}_{\equiv\; \psi_z} \hat{z}_t \]

\[ \boxed{\hat{c}_t = \psi_k \, \hat{k}_t + \psi_z \, \hat{z}_t} \]

• Once we know \((\phi_k, \phi_z)\), consumption is pinned down by the resource constraint

Expected Future Consumption and Return

• Using \(\mathbb{E}_t[\hat{z}_{t+1}] = \rho \, \hat{z}_t\) and \(\hat{k}_{t+1} = \phi_k \hat{k}_t + \phi_z \hat{z}_t\):

\[ \mathbb{E}_t[\hat{c}_{t+1}] = \psi_k \underbrace{(\phi_k \hat{k}_t + \phi_z \hat{z}_t)}_{\hat{k}_{t+1}} + \psi_z \underbrace{\rho \, \hat{z}_t}_{\mathbb{E}_t[\hat{z}_{t+1}]} \]

• Define \(\gamma \equiv \tfrac{(\alpha-1)\,r_{ss}}{k_{ss}}\) (the curvature of the MPK). Then:

\[ \mathbb{E}_t[\hat{R}_{t+1}] = \gamma \underbrace{(\phi_k \hat{k}_t + \phi_z \hat{z}_t)}_{\hat{k}_{t+1}} + r_{ss} \underbrace{\rho \, \hat{z}_t}_{\mathbb{E}_t[\hat{z}_{t+1}]} \]

• Both expectations are linear in \((\hat{k}_t, \hat{z}_t)\) — so matching coefficients will pin down \(\phi\)

Finding \(\phi_k\): Match Coefficients on \(\hat{k}_t\)

• Plug into the Euler equation \(\hat{c}_t = \mathbb{E}_t[\hat{c}_{t+1}] - \beta c_{ss}\,\mathbb{E}_t[\hat{R}_{t+1}]\) and collect terms on \(\hat{k}_t\):

\[ \psi_k = \phi_k\bigl(\psi_k - \beta\, c_{ss}\, \gamma\bigr) \]

• Since \(\psi_k = \tfrac{1}{\beta} - \phi_k\), this is one equation in one unknown (\(\phi_k\)):

\[ \boxed{\frac{1}{\beta} - \phi_k = \phi_k\!\left(\frac{1}{\beta} - \phi_k - \beta \, c_{ss}\, \gamma\right)} \]

• This is a quadratic in \(\phi_k\) — pick the stable root (\(|\phi_k| < 1\))

Finding \(\phi_z\): Match Coefficients on \(\hat{z}_t\)

• Once \(\phi_k\) is known (\(\Rightarrow\) \(\psi_k\) is known), collect terms on \(\hat{z}_t\):

\[ \psi_z = \psi_k \phi_z + \psi_z \rho - \beta\, c_{ss}\bigl(\gamma \phi_z + r_{ss}\,\rho\bigr) \]

• Substitute \(\psi_z = y_{ss} - \phi_z\) and collect \(\phi_z\) terms:

\[ \boxed{\phi_z = \frac{y_{ss}(1 - \rho) + \beta\, c_{ss}\, r_{ss}\, \rho}{\psi_k - \beta\, c_{ss}\, \gamma + (1 - \rho)}} \]

• This is a linear equation — just arithmetic once \(\phi_k\) is known

Finding \(\phi_z\): Match Coefficients on \(\hat{z}_t\)

Summary of the Approach

1. Guess \(\hat{k}_{t+1} = \phi_k \hat{k}_t + \phi_z \hat{z}_t\) → consumption \(\hat{c}_t = \psi_k \hat{k}_t + \psi_z \hat{z}_t\) follows from the RC
2. Compute \(\mathbb{E}_t[\hat{c}_{t+1}]\) and \(\mathbb{E}_t[\hat{R}_{t+1}]\) as linear functions of \((\hat{k}_t, \hat{z}_t)\)
3. Match \(\hat{k}_t\) coefficients in the Euler equation → quadratic for \(\phi_k\)
4. Match \(\hat{z}_t\) coefficients → linear equation for \(\phi_z\)

Julia Implementation

Solving for \(\phi_k\)

• Use NonlinearSolve.jl to solve the \(\phi_k\) equation, then compute \(\phi_z\) directly:

using NonlinearSolve

function linearize_rbc(m)
    (; α, β, δ, ρ) = m
    ss_vals = steady_state(m)
    k_ss = ss_vals.k;  y_ss = ss_vals.y;  c_ss = ss_vals.c
    r_ss = α * k_ss^(α-1)
    γ = (α - 1) * r_ss / k_ss                        # curvature of MPK

    # Step 1: solve for ϕ_k  (ψ_k = ψ_k·ϕ_k − β·c_ss·γ·ϕ_k)
    function ϕk_residual(x, p)
        ψ_k = 1/β - x[1]
        return [ψ_k - x[1] * (ψ_k - β * c_ss * γ)]
    end
    prob = NonlinearProblem(ϕk_residual, [0.9])
    ϕ_k = solve(prob, NewtonRaphson()).u[1]

    # Step 2: compute ϕ_z directly
    ψ_k = 1/β - ϕ_k
    ϕ_z = (y_ss*(1 - ρ) + β*c_ss*r_ss*ρ) / (ψ_k - β*c_ss*γ + (1 - ρ))

    return (ϕ_k=ϕ_k, ϕ_z=ϕ_z)
end

linearize_rbc (generic function with 1 method)

The Linearized Policy Rule

lin = linearize_rbc(m)
println("ϕ_k = $(round(lin.ϕ_k, digits=4))")
println("ϕ_z = $(round(lin.ϕ_z, digits=4))")

ϕ_k = 0.8784
ϕ_z = 1.0427

• The linearized capital policy:

\[ k_{t+1} \approx k_{ss} + \phi_k \, (k_t - k_{ss}) + \phi_z \, (z_t - 1) \]

• \(|\phi_k| < 1\): capital mean-reverts toward the steady state
• \(\phi_z > 0\): higher productivity \(\rightarrow\) save more

Comparing to VFI

Policy Functions

s_mid = (m.N_z + 1) ÷ 2
kp_linear = ss.k .+ lin.ϕ_k .* (m.k_grid .- ss.k) .+ lin.ϕ_z .* (m.z_grid[s_mid] - 1.0)

plot(m.k_grid, sol.kp[:, s_mid], linewidth=2, label="VFI",
     legend=:topleft, size=(700, 400))
plot!(m.k_grid, kp_linear, linewidth=2, linestyle=:dash, label="Linearized")
plot!(m.k_grid, m.k_grid, color=:gray, linestyle=:dot, label="45° line")
vline!([ss.k], color=:gray, linestyle=:dashdot, label="k_ss", alpha=0.5)
xlabel!("Current Capital k")
ylabel!("Next Period Capital k'")
title!("Savings Policy: VFI vs. Linearized (median z)")

• Near \(k_{ss}\): the two solutions nearly coincide
• Far from \(k_{ss}\): the linearized solution diverges — it can’t capture curvature

Policy Functions Across Productivity Levels

plot(legend=:topleft, size=(700, 400), title="VFI (solid) vs. Linearized (dashed)")
for s in [1, 4, 7]
    z = m.z_grid[s]
    kp_lin = ss.k .+ lin.ϕ_k .* (m.k_grid .- ss.k) .+ lin.ϕ_z .* (z - 1.0)
    plot!(m.k_grid, sol.kp[:, s], linewidth=2,
          label="z = $(round(z, digits=3))", color=palette(:default)[s ÷ 3 + 1])
    plot!(m.k_grid, kp_lin, linewidth=2, linestyle=:dash,
          label="", color=palette(:default)[s ÷ 3 + 1])
end
plot!(m.k_grid, m.k_grid, color=:gray, linestyle=:dot, label="45°")
xlabel!("Current Capital k")
ylabel!("Next Period Capital k'")

Simulating the Linearized Model

function simulateRBC_linear(m, lin, T; seed=42)
    Random.seed!(seed)
    (; α, δ, z_grid, P_z, N_z) = m
    k_ss = steady_state(m).k
    k_sim = zeros(T);  s_sim = zeros(Int, T)
    y_sim = zeros(T);  c_sim = zeros(T);  i_sim = zeros(T)
    k_sim[1] = k_ss
    s_sim[1] = (N_z + 1) ÷ 2
    for t in 1:T
        z = z_grid[s_sim[t]]
        kp = k_ss + lin.ϕ_k * (k_sim[t] - k_ss) + lin.ϕ_z * (z - 1.0)
        y_sim[t] = z * k_sim[t]^α
        c_sim[t] = y_sim[t] + (1-δ)*k_sim[t] - kp
        i_sim[t] = kp - (1-δ)*k_sim[t]
        if t < T
            k_sim[t+1] = kp
            s_sim[t+1] = rand(Categorical(P_z[s_sim[t], :]))
        end
    end
    return DataFrame(t=1:T, k=k_sim, z=z_grid[s_sim], y=y_sim, c=c_sim, i=i_sim)
end

simulateRBC_linear (generic function with 1 method)

Comparing Simulated Output

sim_vfi = simulateRBC(m, sol, 200)
sim_lin = simulateRBC_linear(m, lin, 200)

plot(sim_vfi.y, linewidth=2, label="VFI", legend=:topright, size=(700, 400))
plot!(sim_lin.y, linewidth=2, linestyle=:dash, label="Linearized")
xlabel!("Period")
ylabel!("Output")
title!("Simulated Output: VFI vs. Linearized")

• The two simulations track each other very closely
• Near the steady state, linearization is an excellent approximation

Comparing Simulated Consumption and Investment

p1 = plot(sim_vfi.c, linewidth=2, label="VFI", title="Consumption", size=(700, 200))
plot!(p1, sim_lin.c, linewidth=2, linestyle=:dash, label="Linearized")
p2 = plot(sim_vfi.i, linewidth=2, label="VFI", title="Investment", size=(700, 200))
plot!(p2, sim_lin.i, linewidth=2, linestyle=:dash, label="Linearized")
plot(p1, p2, layout=(2,1), size=(700, 450))

Linearization: Pros and Cons

Pros: Almost instant to compute, gives closed-form policy rules, easy to analyze theoretically.

Cons: Only accurate near the steady state, misses nonlinearities (e.g., precautionary savings, occasionally binding constraints).

Most of the macro literature uses linearization (or higher-order perturbation). VFI is the gold standard for accuracy but is computationally demanding.

Summary

Key Takeaways

Concepts

• Real Business Cycle theory
• Calibration: matching long-run moments and the Solow residual
• Stochastic neoclassical growth model
• Social planner’s Bellman equation
• Linear interpolation for continuous states
• Rouwenhorst discretization of AR(1) processes
• Business cycle statistics: volatilities, correlations
• Linearization via method of undetermined coefficients and root solving

Julia Tools

• Interpolations.jl: linear_interpolation with call syntax
• Optim.jl: Brent() for 1D optimization
• NonlinearSolve.jl: NonlinearProblem, NewtonRaphson()
• DataFrames.jl: tabular simulation output
• Mutable structs with setup! and (; ...) = m destructuring
• Categorical for Markov simulation
• Broadcasting for element-wise operations