Kydland & Prescott (1982): Real Business Cycles

1 Introduction

1.1 From Micro to Macro

In the McCall search model, an unemployed worker receives random wage offers and decides each period whether to accept or keep searching. In Rust (1987), a bus superintendent observes accumulated mileage and decides whether to replace the engine or keep running it. Both are single-agent models set in a micro environment: one decision-maker, one state variable, one choice.

Kydland & Prescott (1982) scale this up dramatically. The agent is a representative household in a macroeconomy: it owns capital, supplies labor, and makes consumption and savings decisions. Firms produce using capital and labor, and productivity fluctuates randomly over time. The question is whether a perfectly competitive economy hit by productivity shocks can replicate the business cycles we observe in the data.

“Time to Build and Aggregate Fluctuations” was published in Econometrica in 1982 (now with over 10,000 citations) and earned Finn Kydland and Edward Prescott the 2004 Nobel Prize in Economics.

ImportantThe RBC Revolution

Before KP (1982), business cycles were thought to require market failures: sticky prices, coordination failures, or demand shortfalls. KP showed that a perfectly competitive economy hit by productivity shocks could generate realistic fluctuations without any friction. This launched the Real Business Cycle (RBC) research program.

The core model is a stochastic neoclassical growth model: households save and consume optimally, firms produce using capital and labor, and productivity fluctuates randomly. The model combines dynamic programming (as in McCall and Rust), Markov chains (for the productivity process), and linear algebra (for computing expectations).

A key computational challenge distinguishes the RBC model from everything we have done so far. In McCall and Rust, the state space was discrete: we stored the value function as a vector and looked up values directly. In the RBC model, capital \(k\) is continuous. We cannot simply index into a vector to evaluate the value function at an arbitrary capital level. This lecture introduces linear interpolation, the tool that lets us handle continuous state spaces.

1.2 Roadmap

These notes cover six main topics:

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

2 The Stochastic Growth Model

2.1 The Household

A representative household lives forever and maximizes expected discounted lifetime utility:

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

where \(c_t\) is consumption in period \(t\) and \(\beta \in (0,1)\) is the discount factor. The \(\mathbb{E}_0\) denotes the expectation conditional on information available at time 0—the household does not know future productivity but makes the best plan given what it does know.

Why log utility? The function \(\log(c_t)\) has two important properties. First, it is increasing: more consumption is always better. Second, it is concave: the marginal utility of an extra dollar of consumption \(1/c_t\) falls as \(c_t\) rises. A household consuming $100 values an extra dollar much more than a household already consuming $10,000. This concavity is what drives consumption smoothing: the household prefers a steady consumption path to one that swings between feasts and famines.

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

1. Rental income: \(r_t k_t\) — the household rents its capital stock to firms at a competitive rental rate \(r_t\)
2. Wage income: \(w_t\) — the household supplies one unit of labor at the competitive wage

The household takes both prices \((r_t, w_t)\) as given. It is a price taker in competitive factor markets. Each period, it faces the budget constraint:

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

The left side describes uses of funds: the household allocates its resources between consumption today (\(c_t\)) and capital carried forward to tomorrow (\(k_{t+1}\)). The right side describes sources of funds: wage income (\(w_t\)), rental income (\(r_t k_t\)), and undepreciated capital (\((1-\delta) k_t\)). The parameter \(\delta \in (0,1)\) is the depreciation rate: the fraction of the capital stock that wears out each period. If \(\delta = 0.10\), then 10% of the machines break down each year and must be replaced just to keep the capital stock constant.

The household’s problem is to choose sequences of consumption \(\{c_t\}\) and capital \(\{k_{t+1}\}\) to maximize expected lifetime utility subject to this budget constraint in every period. This is a dynamic optimization problem: today’s savings decision affects the capital available tomorrow, which affects future income and future choices. The Bellman equation, which we will derive shortly, captures this recursive structure.

TipSimplification

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 features, and the simplified model is considerably easier to solve.

2.2 The Firm

A representative firm produces output using capital and labor according to a Cobb–Douglas production function:

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

This specification has three components. The variable \(z_t\) is total factor productivity (TFP): the random shock that drives business cycles. When \(z_t\) is high, the same inputs produce more output; when \(z_t\) is low, the economy is less productive. The parameter \(\alpha \in (0,1)\) is capital’s share of output, typically about \(1/3\) based on U.S. national accounts data. The remaining share \((1-\alpha) \approx 2/3\) goes to labor.

The production function exhibits constant returns to scale (CRS): if you double both capital and labor, output exactly doubles. Mathematically, \(f(\lambda k, \lambda \ell) = \lambda f(k, \ell)\) for any \(\lambda > 0\). This ends up being a natural assumption for an aggregate economy. If the U.S. had twice as many workers and twice as many machines, we would expect roughly twice the output.

The firm rents capital at rate \(r_t\) and hires labor at wage \(w_t\), taking both prices as given. It chooses inputs to maximize profits:

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

The firm’s first-order conditions set the marginal product of each factor equal to its price. Taking the derivative of the production function with respect to \(k_t\) and \(\ell_t\):

\[ 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} \]

The rental rate equals the marginal product of capital (MPK) and the wage equals the marginal product of labor (MPL). Each factor is paid exactly what it contributes at the margin.

TipZero Profits

With constant returns to scale, 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. This is Euler’s theorem applied to homogeneous functions. You can verify this directly: \(\alpha z_t k_t^{\alpha} \ell_t^{1-\alpha} + (1-\alpha) z_t k_t^{\alpha} \ell_t^{1-\alpha} = z_t k_t^{\alpha} \ell_t^{1-\alpha} = y_t\).

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} \]

These expressions reveal two important economic forces. First, the rental rate falls with capital: the exponent \(\alpha - 1 < 0\) means that as the economy accumulates more capital, each additional unit is less productive. This is diminishing returns to capital, and it is what prevents the economy from accumulating capital without bound. Second, the wage rises with capital: more capital per worker means each worker is more productive (they have more machines to work with). Both factor prices move proportionally with \(z_t\): when productivity is high, every factor of production is more valuable.

2.3 Technology Shocks

What drives business cycles in this model? The answer is total factor productivity (TFP), \(z_t\). Some periods the economy is highly productive: new technologies diffuse, weather is favorable, supply chains run smoothly. Other periods it is less so. The RBC model takes the position that these productivity fluctuations are the primary source of aggregate fluctuations.

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) \]

The parameter \(\rho \in (0,1)\) governs the persistence of technology shocks. With \(\rho = 0.90\), a 1% productivity innovation today still predicts a \(0.90^4 \approx 0.66\%\) deviation four periods later. This persistence is critical: if shocks were purely temporary (\(\rho = 0\)), they would cause one-period blips in output, not the sustained expansions and contractions we observe in the data. It is the persistence of productivity shocks that turns one-period disturbances into multi-period business cycles.

The parameter \(\sigma_\varepsilon\) controls the volatility of innovations: how large the surprises are in each period. The typical estimate of \(\sigma_\varepsilon \approx 0.02\) means that a one-standard-deviation shock moves productivity by about 2%. This may sound small, but because the shocks are persistent and the economy amplifies them through savings and investment responses, the resulting output fluctuations end up being comparable to what we see in the data.

Why model TFP as an AR(1)? First, it is empirically grounded: regressing the Solow residual on its own lag yields a tight fit with high \(R^2\). Second, the AR(1) is parsimonious: it captures both the volatility and persistence of productivity with just two parameters. Third, the process is stationary when \(|\rho| < 1\), so the model has a well-defined long-run distribution around which business cycles fluctuate.

ImportantConnecting to Markov Chains

The AR(1) process for \(\log z_t\) is continuous: \(z_t\) can take any positive value. To solve the model on a computer, we discretize this process into a finite-state Markov chain with \(N_z\) states and a transition matrix \(P_z\), just like we studied in the Markov chain lectures. This lets us compute the conditional expectation \(\mathbb{E}[V(k', z') \mid z]\) as a finite sum \(\sum_{j=1}^{N_z} P_{ij} V(k', z_j)\) using the Markov chain tools we’ve developed.

2.4 Competitive Equilibrium

We now have all the pieces of the model: a household that saves and consumes, a firm that produces using capital and labor, and a stochastic process for productivity. The next step is to define what it means for all these pieces to fit together consistently.

In a competitive equilibrium, every agent optimizes, every market clears, and no one has an incentive to deviate. Formally:

ImportantDefinition: 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\))

Condition 3 is the goods market clearing condition: every unit of output is either consumed or invested. Investment is \(i_t = k_{t+1} - (1-\delta)k_t\), which is the new capital created beyond what survives depreciation. There is no government, no foreign trade, and no inventory accumulation, so \(y_t = c_t + i_t\) must hold in every period.

Now here is the key simplification. 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 \]

The wage and rental income together equal total output (by Euler’s theorem and zero profits), so this collapses to the aggregate resource constraint:

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

Notice what happened: prices \((r_t, w_t)\) have disappeared entirely. The household’s budget constraint, once we substitute in the competitive factor prices, becomes identical to the economy-wide resource constraint. This means we don’t have to worry about satisfying the household budget constraint if we work directly with the resource constraint for the entire economy. In principle, computing a competitive equilibrium requires tracking household optimization, firm optimization, and market clearing simultaneously. But thanks to the First Welfare Theorem, there is a much simpler approach.

2.5 The First Welfare Theorem

The First Welfare Theorem is one of the most powerful results in economics. It states that under a set of conditions—complete markets, no externalities, and price-taking behavior—a competitive equilibrium allocation is Pareto optimal. An allocation is Pareto optimal if there is no way to make any agent better off without making someone else worse off.

This helps us solve for a competitive equilibrium, because a Pareto optimal allocation is exactly what a benevolent social planner would choose. A benevolent planner maximizes the representative household’s welfare subject only to the aggregate resource constraint. Since the competitive equilibrium is Pareto optimal, and since there is only one household (the representative household), the competitive equilibrium allocation must be identical to the planner’s solution.

This equivalence gives us a powerful computational shortcut. Instead of solving the full competitive equilibrium — which requires finding prices, checking that households optimize given those prices, and verifying that markets clear — we can solve a single optimization problem: the social planner’s problem. The planner faces only the aggregate resource constraint and chooses consumption and savings to maximize welfare. Once we have the planner’s allocation, we can recover the equilibrium prices from the firm’s first-order conditions.

ImportantThe Shortcut

The competitive equilibrium allocation equals the planner’s solution. 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 back them out afterward: \(r_t = \alpha z_t k_t^{\alpha-1}\) and \(w_t = (1-\alpha) z_t k_t^\alpha\).

This approach would not work if there were distortions: taxes, monopoly power, externalities, or incomplete markets would drive a wedge between the competitive equilibrium and the planner’s solution. Much of modern macroeconomics studies exactly these departures. But the frictionless RBC model satisfies all the conditions of the First Welfare Theorem, so the planner’s problem gives us the equilibrium for free.

2.6 The Social Planner’s Problem

The planner’s problem is much easier to work with than the full competitive equilibrium. The planner does not need to think about prices at all—they simply choose how much to save (\(k_{t+1}\)) to maximize the household’s welfare subject to the aggregate resource constraint. The state of the economy is the pair \((k_t, z_t)\): current capital and current productivity. The control is \(k_{t+1}\), the amount of capital carried into the next period. Once \(k_{t+1}\) is chosen, consumption is pinned down by the resource constraint: \(c_t = z_t k_t^\alpha + (1-\delta)k_t - k_{t+1}\).

This is a dynamic programming problem with exactly the same recursive structure as the McCall and Rust models. The key insight is that the planner faces the same trade-off every period: consuming more today is pleasant, but saving more increases the capital stock and raises future output. The optimal plan balances these two forces in every state of the world.

The Bellman equation for the planner’s problem is:

\[ V(k, z) = \max_{k'} \left\{ \log\left(z k^\alpha + (1-\delta) k - k'\right) + \beta \, \mathbb{E}\left[V(k', z') \mid z\right] \right\} \]

subject to \(0 < k' \leq z k^\alpha + (1-\delta) k\) (consumption must be positive). The first term is today’s utility from consumption, where consumption equals total resources minus savings. The second term is the discounted expected future value, where the expectation is over tomorrow’s productivity \(z'\) conditional on today’s \(z\). The planner chooses \(k'\) to make the best trade-off between these two terms.

Notice that the choice variable \(k'\) is continuous: the planner can choose any savings level, not just “accept” or “reject.” This makes the optimization problem richer but also more computationally demanding.

The following table compares the RBC planner’s problem to the models we have solved previously:

	McCall	Rust	RBC
State	Wage offer (discrete)	Mileage (discrete)	Capital \(\times\) Productivity
Choice	Accept/Reject	Replace/Keep	How much to save
State space	Finite	Finite	Continuous \(\times\) Finite
Continuation	\(\sum_j P_{ij} V_j\)	\(\sum_j F_{ij} V_j\)	\(\sum_j P_{ij} V(k', z_j)\)

The key difference is in the last row. In McCall and Rust, the continuation value was a simple weighted sum over a finite vector: we could look up \(V_j\) by indexing into an array. In the RBC model, the continuation value involves \(V(k', z_j)\) where \(k'\) is a continuous variable that will generally not fall on a grid point. We cannot simply index into a vector to evaluate the value function at an arbitrary capital level. This is why we need interpolation: a method for approximating \(V\) between grid points.

2.7 The Deterministic Steady State

Before solving the full stochastic model, it is useful to understand what happens when productivity is 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, investment exactly replaces worn-out capital (\(i_{ss} = \delta k_{ss}\)) and consumption equals output net of depreciation (\(c_{ss} = k_{ss}^\alpha - \delta k_{ss}\)).

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

2.7.1 Deriving the Steady State

Suppose the planner is at the steady state and considers saving one extra unit of capital today. The cost is consuming one less unit today, which reduces utility by \(1/c_{ss}\). The benefit is that the extra unit of capital produces \(\alpha k_{ss}^{\alpha-1}\) extra output tomorrow, plus the \((1-\delta)\) undepreciated capital comes back. All of this is discounted by \(\beta\) and valued at \(1/c_{ss}\):

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

At the optimum, the planner cannot improve by saving more or less, so cost equals benefit:

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

The \(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}}} \]

The term \(\alpha k_{ss}^{\alpha-1}\) is the marginal product of capital: the steady-state rental rate \(r_{ss}\). The condition can be rewritten as \(r_{ss} = 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 until the planner is indifferent between consuming and saving.

TipThe 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.

2.8 Setup in Julia

We now turn to implementing the model in Julia. The setup involves three main tasks: loading packages, discretizing the productivity process, and bundling the model parameters into a reusable struct. This infrastructure will support everything that follows: value function iteration, simulation, and linearization.

2.8.1 Loading Packages

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

2.8.2 Discretizing the AR(1): Rouwenhorst Method

The AR(1) process for \(\log z_t\) is continuous, but our Bellman equation requires computing expectations as finite sums over a Markov chain. We need to approximate the continuous process with an \(N\)-state Markov chain that captures its persistence and volatility.

The Rouwenhorst (1995) method does this in two steps:

1. Grid construction: create an evenly spaced grid for \(\log z\) spanning \(\pm\sqrt{N-1} \cdot \sigma_z\), where \(\sigma_z = \sigma_\varepsilon / \sqrt{1 - \rho^2}\) is the unconditional standard deviation of the process. This ensures the grid covers the same range as the true stationary distribution.
2. Transition matrix: build \(P_z\) recursively, starting from a \(2 \times 2\) base case with parameter \(p = (1+\rho)/2\). Each recursion step adds one state by combining four shifted copies of the previous matrix. The resulting transition matrix captures the AR(1) persistence exactly for any \(N\).

Why Rouwenhorst over the more common Tauchen (1986) method? Rouwenhorst performs much better for highly persistent processes (\(\rho\) close to 1), which is exactly the empirically relevant case for RBC models. Tauchen’s method can severely understate the persistence of the discretized process when \(\rho > 0.9\), leading to business cycles that are too short-lived.

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,
         sqrt(N-1) * σ_z, N))
    return z_grid, P
end

rouwenhorst (generic function with 1 method)

Let’s test the function with a 7-state approximation using the calibrated parameter values (\(\rho = 0.90\), \(\sigma_\varepsilon = 0.02\)). The function returns the grid in log levels, so we exponentiate to get the productivity levels \(z_j\):

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

The seven productivity states range from about 0.91 (a recession) to about 1.10 (a boom), centered at 1.0 (average productivity). These fluctuations of roughly \(\pm 10\%\) around the mean are consistent with the range of Solow residual variation in U.S. data.

The transition matrix shows how likely the economy is to move between productivity states:

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

The high diagonal values (the probability of staying in the same state next period) confirm that the economy tends to remain near its current productivity level. For example, if the economy is in the median state (\(z \approx 1.0\)), it has roughly a 70% chance of staying there and small probabilities of moving to adjacent states. This persistence is exactly what we need: a good shock today predicts good productivity tomorrow, generating the sustained expansions and contractions we observe in the data.

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

2.8.3 The RBCModel Struct

As models grow more complex, passing individual parameters to every function becomes unwieldy and error-prone. Instead, we bundle everything into a struct: a custom data type that holds all the model’s parameters and computed objects in one place. We use @kwdef mutable struct, which gives us two conveniences: @kwdef lets us set default values and construct instances using keyword arguments, and mutable lets us modify fields after creation (which we need for filling in grids).

We use a mutable struct to hold the model parameters and the objects we compute from them:

@kwdef mutable struct RBCModel
    β::Float64 = 0.96           # discount factor
    α::Float64 = 0.33           # capital share
    δ::Float64 = 0.10           # depreciation rate
    ρ::Float64 = 0.90           # TFP persistence
    σ_ε::Float64 = 0.02         # TFP innovation std dev
    N_z::Int = 7                # productivity states
    N_k::Int = 200              # 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 struct separates two kinds of data. The economic parameters (\(\beta, \alpha, \delta, \rho, \sigma_\varepsilon\)) and grid sizes (\(N_z, N_k\)) have sensible defaults that the user can override. The computed objects (the productivity grid z_grid, transition matrix P_z, and capital grid k_grid) start empty and are filled by a setup! function. This two-stage pattern (create, then initialize) ensures that all computed objects are consistent with the current parameter values:

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)

The ! suffix signals that setup! mutates the struct, following Julia’s convention for functions that modify their arguments. The (; α, β, δ, ...) = m syntax destructures the struct fields into local variables, keeping the code readable. Notice that the capital grid is centered on the steady state \(k_{ss}\), ranging from \(0.2 k_{ss}\) to \(1.8 k_{ss}\). This ensures the grid covers the range where the economy actually spends time during business cycle fluctuations.

If you change a parameter (e.g., m.ρ = 0.95), just call setup!(m) again to rebuild the grids and transition matrix consistently with the new value.

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

RBCModel(0.96, 0.33, 0.1, 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.7065757834312844, 0.7349808400516376, 0.7633858966719906, 0.7917909532923438, 0.8201960099126969, 0.8486010665330501, 0.8770061231534032, 0.9054111797737563, 0.9338162363941096, 0.9622212930144627  …  6.10353654129838, 6.131941597918734, 6.1603466545390875, 6.18875171115944, 6.217156767779794, 6.245561824400147, 6.2739668810205, 6.302371937640853, 6.330776994261205, 6.359182050881559])

Let’s inspect what setup! produced:

println("z_grid: ", round.(m.z_grid, digits=4))
k1 = round(m.k_grid[1], digits=3)
kN = round(m.k_grid[end], digits=3)
println("k_grid: [$k1, ..., $kN]")
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.707, ..., 6.359]
N_k = 200 grid points

The productivity grid matches what we computed manually above. The capital grid runs from well below to well above the steady state, with 200 evenly spaced points providing fine resolution for interpolation.

2.8.4 The Steady State

We also create a helper function that computes the deterministic steady-state values from the model parameters. This is useful for calibration checks, grid construction, and as a starting point for simulation:

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))")
println("y_ss = $(round(ss.y, digits=3))")
println("c_ss = $(round(ss.c, digits=3))")
println("r_ss = $(round(ss.r, digits=3))")

k_ss = 3.533
y_ss = 1.517
c_ss = 1.163
r_ss = 0.142

2.9 Calibration

2.9.1 Why Calibrate?

The model has five parameters: \((\alpha, \beta, \delta, \rho, \sigma_\varepsilon)\). How do we choose their values? In many areas of economics, we would estimate parameters econometrically: maximizing a likelihood function or minimizing a distance metric. The RBC literature instead calibrates: it picks parameter values so that the model’s steady state matches long-run averages in the data.

The key principle behind calibration is separation of evidence. We use well-established long-run facts (growth accounting, interest rates, investment rates) to pin down the parameters, and then ask whether the model (without any further tuning) can reproduce business cycle fluctuations it was not calibrated to match. If the model generates realistic output volatility, consumption smoothing, and investment volatility using only parameters that match long-run averages, that is strong evidence that the model captures something real about how the economy works.

TipCalibration vs. Estimation

Calibration is not estimation. We are not maximizing a likelihood or minimizing a distance from the data. We are using economic theory to set each parameter at a value consistent with a specific, well-established long-run fact, and then evaluating the model on its ability to reproduce short-run dynamics. This makes the exercise a genuine test of the model: can a few structural parameters, pinned down by slow-moving averages, explain fast-moving business cycle fluctuations?

2.9.2 Matching Long-Run Moments

Each parameter is pinned down by a different data moment, as summarized in the table:

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 + \delta \approx 4\%\)/yr	0.96
\(\delta\)	Investment–capital ratio	\(I/K \approx 10\%\)/yr	0.10

• \(\alpha = 0.33\): The National Income and Product Accounts (NIPA) decompose national income into compensation of employees and returns to capital. Labor’s share of income has historically averaged approximately \(2/3\) in the United States, so capital’s share is approximately \(1/3\). In the Cobb–Douglas production function, capital’s share is exactly \(\alpha\), so we set \(\alpha = 0.33\).
• \(\beta = 0.96\): The steady-state Euler equation gives \(r_{ss} = 1/\beta - (1-\delta)\). The average annual real return on capital in the U.S. is about 4%. With \(\delta = 0.10\), solving for \(\beta\) gives \(\beta = 1/(r_{ss} + 1 - \delta) = 1/(0.04 + 0.90) \approx 0.96\). A discount factor of 0.96 means the household values a dollar next year at 96 cents today.
• \(\delta = 0.10\): The Bureau of Economic Analysis estimates that the aggregate capital stock depreciates at roughly 10% per year. This means that about 10% of all machines, buildings, and equipment must be replaced each year just to keep the capital stock constant.

2.9.3 The Productivity Process

The remaining two parameters, \(\rho\) and \(\sigma_\varepsilon\), govern the stochastic productivity process. Unlike \(\alpha\), \(\beta\), and \(\delta\), which come from long-run averages, these are estimated directly from time-series data.

The procedure starts by constructing the Solow residual: the part of output that cannot be explained by measured inputs. Given the production function \(Y_t = z_t K_t^\alpha L_t^{1-\alpha}\), we can back out productivity as:

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

With data on output (\(Y_t\)), the capital stock (\(K_t\)), and hours worked (\(L_t\)), and using our calibrated \(\alpha = 0.33\), we compute the Solow residual series. Regressing \(\log z_t\) on \(\log z_{t-1}\) by OLS gives:

• \(\hat{\rho} \approx 0.90\): productivity shocks are very persistent: a 1% shock today still predicts a 0.53% deviation two periods later (\(0.9^2 \approx 0.81\), adjusted for the regression)
• \(\hat{\sigma}_\varepsilon \approx 0.02\): innovations are small (about 2% of output per period)

r_impl = round(1/m.β - 1 + m.δ, digits=3)
println("Calibrated Parameters:")
println("  α   = $(m.α)   (capital share)")
println("  β   = $(m.β)   (discount factor → r ≈ $r_impl)")
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.142)
  δ   = 0.1   (depreciation rate)
  ρ   = 0.9   (TFP persistence)
  σ_ε = 0.02 (TFP innovation std dev)

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

Implied Steady State:
  K/Y = 2.33   (data: ≈ 3.0)
  I/Y = 0.23   (data: ≈ 0.25-0.30)
  C/Y = 0.77   (data: ≈ 0.65-0.70)

ImportantThe 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?

3 Linear Interpolation

3.1 Why We Need Interpolation

Before we can solve the planner’s Bellman equation, we need a new computational tool. In the McCall and Rust models, the state space was discrete: we stored the value function as a vector and looked up values by index. If the state was wage offer \(j\), we simply read \(V[j]\) from the vector. No approximation was needed.

In the RBC model, capital \(k\) is continuous. We can store \(V\) on a finite grid of \(N_k\) points: \(V(k_1, z), V(k_2, z), \ldots, V(k_N, z)\). But when we optimize the Bellman equation, the optimizer will try many candidate values of \(k'\), and the optimal \(k'\) will almost certainly not land exactly on a grid point. We need to evaluate \(V(k', z)\) for any \(k'\) in the range, not just the grid values.

Interpolation is the solution: use the known grid values to approximate \(V\) at off-grid points. The idea is simple—if we know the value function at \(k = 3.0\) and \(k = 3.5\), we can make a reasonable guess about its value at \(k = 3.27\) by drawing a straight line between the two known points.

We know \(V\) at the blue dots—but optimization may require \(V(2.7)\), which lies between grid points.

3.2 How Linear Interpolation Works

Given two adjacent grid points \(k_i\) and \(k_{i+1}\) with known values \(V_i\) and \(V_{i+1}\), linear interpolation draws a straight line between them and reads off the value at \(k'\):

\[ 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} \in [0,1]\), this becomes a weighted average:

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

The algorithm has three steps:

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}\)

When \(w = 0\), \(k'\) is exactly at \(k_i\) and \(\hat{V} = V_i\). When \(w = 1\), \(k'\) is at \(k_{i+1}\) and \(\hat{V} = V_{i+1}\). When \(w = 0.5\), \(k'\) is at the midpoint and \(\hat{V}\) is the average.

3.3 Interpolations.jl

We could implement linear interpolation by hand—find the interval, compute the weight, take the weighted average. But Julia’s Interpolations.jl package does all of this automatically and returns an object with clean call syntax that behaves like a function:

# 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)

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

The object V_itp behaves like a function: call it with any \(k'\) to get an interpolated value.

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

TipCall 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 that tries many candidate \(k'\) values.

We can visualize the interpolation object over a fine grid using broadcasting:

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. With more grid points, the piecewise-linear approximation better captures curvature. For smooth economic functions, 100–500 grid points typically give excellent accuracy.

4 Solving the RBC Model

4.1 Value Function Iteration

We solve the planner’s Bellman equation using value function iteration (VFI), the same approach we used for the McCall and Rust models. The key difference is that the continuous capital state requires interpolation at every step. The algorithm is:

1. Start with an initial guess \(V^0(k, z) = 0\) for all \((k, z)\)
2. For each grid point \((k_i, z_j)\), find \(k'\) that maximizes \[\log\bigl(z_j k_i^\alpha + (1-\delta)k_i - k'\bigr) + \beta \sum_m P_{jm} V^n(k', z_m)\] using linear interpolation to evaluate \(V^n(k', z_m)\) at non-grid values of \(k'\)
3. Repeat until \(\|V^{n+1} - V^n\|_\infty < \text{tol}\)

ImportantThe Role of Interpolation

In step 2, the optimizer tries many candidate \(k'\) values. Most will not be on our grid. Linear interpolation lets us evaluate the continuation value \(V(k', z')\) for any \(k'\).

4.2 Solving for the Optimal Policy at a Single State

The core building block of our solver is the function optimal_policy. This function takes a single state (current capital \(k\) and current productivity index iz) along with the current value function (stored as interpolation objects V_itp), and finds the optimal savings \(k'\) that maximizes the right-hand side of the Bellman equation at that state.

Think of it this way: for each point \((k_i, z_j)\) on our grid, we need to solve a one-dimensional optimization problem. The planner has total resources \(z_j k_i^\alpha + (1-\delta)k_i\) and must split them between consumption today and capital tomorrow. optimal_policy solves this split for a single \((k, z)\) pair.

The first step is computing total available resources:

function optimal_policy(m, k, iz, V_itp)
    (; β, α, δ, z_grid, P_z, N_z, k_grid) = m
    z = z_grid[iz]
    budget = z * k^α + (1 - δ) * k        # y + undepreciated k
    # ...
end

budget is the maximum the planner could carry forward to next period (setting \(c = 0\)). Any feasible \(k'\) must satisfy \(k' < \text{budget}\).

For a candidate \(k'\), the planner’s payoff is:

    function obj(kp)
        c = budget - kp                              # consumption = what's left
        c <= 0 && return 1e10                        # infeasible: penalty

        EV = sum(P_z[iz, jz] * V_itp[jz](kp) for jz in 1:N_z) # E[V(k', z')]

        return -(log(c) + β * EV)                    # negative (we minimize)
    end

V_itp[jz](kp) calls the interpolation object for productivity state \(j\) to evaluate \(V(k', z_j')\) at any \(k'\). This is where linear interpolation does its work. We sum over all \(N_z\) future productivity states, weighted by transition probabilities \(P_z[iz, jz]\). The c <= 0 && return 1e10 line handles the feasibility constraint: if the candidate \(k'\) exceeds the budget, consumption would be negative, so we return a large penalty to steer the optimizer away. The minus sign converts the problem from maximization to minimization, since Julia’s Optim.jl package minimizes by convention.

The full function uses Brent’s method for one-dimensional optimization. Brent’s method is a reliable bracketed search algorithm that combines the safety of bisection with the speed of inverse quadratic interpolation. It is ideal here because our objective function is smooth, one-dimensional, and we know the bounds: \(k'\) must lie between the minimum grid point and the budget (minus a tiny \(\varepsilon = 10^{-8}\) to ensure positive consumption):

function optimal_policy(m, k, iz, V_itp)
    (; β, α, δ, z_grid, P_z, N_z, k_grid) = m
    z = z_grid[iz]
    budget = z * k^α + (1 - δ) * k                  # total resources

    function obj(kp)
        c = budget - kp                  # consumption = what's left
        c <= 0 && return 1e10        #Penalize negative consumption

        EV = sum(P_z[iz, jz] * V_itp[jz](kp) for jz in 1:N_z)

        return -(log(c) + β * EV) # negative (we minimize)
    end

    k_max = min(k_grid[end], budget - 1e-8)
    res = optimize(obj, k_grid[1], k_max, Brent())
    return (V = -Optim.minimum(res),                 # optimal value
            kp = Optim.minimizer(res),                # optimal next-period capital
            c = budget - Optim.minimizer(res))        # implied consumption
end

optimal_policy (generic function with 1 method)

The function returns a named tuple (V, kp, c) with the optimal value, the optimal capital choice, and implied consumption for that state.

Let’s test it with the initial guess \(V = 0\) at the steady state:

V_init = zeros(m.N_k, m.N_z)
V_init_itp = [linear_interpolation(m.k_grid, V_init[:, jz])
              for jz in 1:m.N_z]

result = optimal_policy(m, ss.k, 4, V_init_itp)
println("V  = $(round(result.V, digits=3))")
println("k' = $(round(result.kp, digits=3))")
println("c  = $(round(result.c, digits=3))")
println("Lower Bound: $(m.k_grid[1]), Upper Bound: $(m.k_grid[end])")

V  = 1.384
k' = 0.707
c  = 3.99
Lower Bound: 0.7065757834312844, Upper Bound: 6.359182050881559

With \(V = 0\) (no future value), the planner behaves as if there were no tomorrow and consumes everything: the optimizer drives \(k'\) all the way to the lower bound of the grid (k_grid[1]) so that consumption equals the entire budget. This makes economic sense: if the future has zero value, there is no reason to save. As VFI iterates and the value function grows, the planner “learns” that the future matters and begins to save more and more, until the value function converges and the savings decision stabilizes.

4.3 The Bellman Operator

The bellman_operator function applies one step of value function iteration by calling optimal_policy at every \((k_i, z_j)\) pair on the grid. It takes the current value function \(V\) (an \(N_k \times N_z\) matrix) and returns the updated value function \(V^{new}\) and the associated policy function \(k'(k, z)\):

function bellman_operator(m::RBCModel, V)
    (; N_k, N_z, k_grid) = m
    V_new = similar(V)
    kp_policy = similar(V)
    c_policy = similar(V)
    V_itp = [linear_interpolation(k_grid, V[:, jz])
             for jz in 1:N_z]
    for iz in 1:N_z
        for ik in 1:N_k
            result = optimal_policy(m, k_grid[ik], iz, V_itp)
            V_new[ik, iz] = result.V
            kp_policy[ik, iz] = result.kp
            c_policy[ik, iz] = result.c
        end
    end
    return V_new, kp_policy, c_policy
end

bellman_operator (generic function with 1 method)

An important efficiency detail: the V_itp interpolation objects are built once per Bellman iteration (before the double loop), not once per \((k_i, z_j)\) state. Since interpolation setup involves sorting and organizing the grid data, doing it once rather than \(N_k \times N_z\) times saves significant computation. The double loop itself is straightforward—all the economic and computational heavy lifting happens inside optimal_policy. We also collect the implied consumption policy c_policy directly here, since optimal_policy already returns it for each state.

4.4 Solving the Model

The solveRBC function ties everything together. It starts with an initial guess of \(V = 0\) (the planner assigns no value to any state), then repeatedly applies the Bellman operator until the value function converges. Convergence is measured by the sup norm: \(\|V^{n+1} - V^n\|_\infty < \text{tol}\), meaning the largest change in the value function across all states is smaller than the tolerance. Each bellman_operator call returns the updated value function alongside the savings and consumption policies, so we do not need a separate pass to compute consumption from the resource constraint. Once converged, the function builds interpolation objects for the policy functions so they can be evaluated at any capital level during simulation:

function solveRBC(m::RBCModel;
                  tol=1e-6, maxiter=1000, verbose=true)
    (; α, δ, z_grid, 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)
    #Iterate the Bellman operator until convergence
    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=3)))")
    end
    #Interpolate the policy functions
    kp_itp = [linear_interpolation(k_grid, kp_policy[:, iz])
              for iz in 1:N_z]
    c_itp  = [linear_interpolation(k_grid, c_policy[:, iz])
              for iz in 1:N_z]
    #Return the solution
    return (V=V, kp=kp_policy, c=c_policy,
            kp_itp=kp_itp, c_itp=c_itp)
end

solveRBC (generic function with 1 method)

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

Converged in 274 iterations (dist = 0.0)

TipPatience

VFI 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.

4.5 Results

With the model solved, we can inspect the value and policy functions to understand how the planner behaves. We plot each function for three productivity levels: low (\(z_1\), the worst state), median (\(z_4\)), and high (\(z_7\), the best state). This shows how the planner’s decisions vary with the current state of the economy.

4.5.1 The Value Function

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

Several features stand out. First, higher productivity \(z\) shifts the value function upward: more output at every capital level means the planner can achieve higher lifetime utility. Second, the value function is concave in \(k\), reflecting diminishing returns to capital: the first unit of capital is extremely valuable, but additional units contribute less and less. Third, the curves are relatively close together, which makes sense because productivity shocks are moderate (\(\pm 10\%\)) and temporary.

4.5.2 The Savings Policy \(k'(k, z)\)

plot(legend=:topleft, size=(700, 400))
for iz in [1, 4, 7]
    plot!(m.k_grid, sol.kp[:, iz], linewidth=2,
          label="z = $(round(m.z_grid[iz], 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="")

The 45-degree line (\(k' = k\)) is the reference line where capital stays constant. When the policy function lies above the 45° line, the planner is saving more than depreciation and capital is growing (\(k' > k\)). When it lies below, the planner is not replacing worn-out capital fast enough and the capital stock is shrinking (\(k' < k\)). The crossing point, where the policy function intersects the 45° line, is the steady state for that productivity level.

Notice how the policy functions shift with productivity. Higher \(z\) means more output, which means the planner has more resources and saves more at every capital level. The high-\(z\) policy function crosses the 45° line at a higher capital level than the low-\(z\) function. This captures an important intuition: during a boom (high \(z\)), the economy accumulates capital; during a recession (low \(z\)), it runs down its capital stock.

4.5.3 The Consumption Policy \(c(k, z)\)

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

Consumption is increasing in both \(k\) and \(z\). Richer economies (higher \(k\)) and more productive economies (higher \(z\)) consume more. Importantly, when productivity is high, the planner can afford to both consume more and save more simultaneously: the extra output is split between current consumption and future capital. The consumption function is relatively smooth and nearly linear, reflecting the household’s desire to smooth consumption over time rather than responding dramatically to each productivity realization.

4.5.4 The Investment Policy

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

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

Investment responds much more strongly to productivity than consumption does. The gap between the high-\(z\) and low-\(z\) investment curves is large relative to the corresponding gap for consumption. Why? Because consumption smoothing means the household absorbs productivity shocks mainly through the investment margin. When \(z\) is high, the household consumes a bit more but channels most of the extra output into investment. When \(z\) is low, the household cuts investment sharply to protect its consumption. This asymmetry between consumption and investment responses foreshadows a key business cycle fact: investment is the most volatile component of GDP.

5 Simulation

5.1 Simulating the Economy

Solving the model gives us policy functions: for any \((k, z)\), we know the optimal \(k'\), \(c\), and \(i\). But to understand what the model implies for business cycles, we need to see these policies in action. Simulation generates artificial time series—quarter-by-quarter paths of output, consumption, investment, and capital—that we can compare to actual macroeconomic data.

The simulation procedure is straightforward. We start the economy at the deterministic steady state \((k_0 = k_{ss}, z_0 = z_{median})\). Each period, we:

1. Compute output, consumption, and investment using the policy function at the current state \((k_t, z_t)\)
2. Draw next period’s productivity \(z_{t+1}\) from the Markov chain (using the transition probabilities from the current state)
3. Set \(k_{t+1}\) to the optimal savings from the policy function
4. Repeat

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); z_idx = 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
    z_idx[1] = (N_z + 1) ÷ 2        # start at median z
    for t in 1:T
        iz = z_idx[t]
        z = z_grid[iz]
        kp = sol.kp_itp[iz](k_sim[t])
        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
            z_idx[t+1] = rand(Categorical(P_z[iz, :]))
        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)

Several implementation details are worth noting. First, sol.kp_itp[iz](k_sim[t]) calls the interpolated policy function at the current capital level, which may not be on our grid. It could be anywhere in the grid range. Interpolation is thus used in two places: inside VFI to evaluate the continuation value, and here in simulation to evaluate the policy function at off-grid capital levels.

Second, rand(Categorical(P_z[iz, :])) draws the next productivity state from the appropriate row of the transition matrix. The Categorical distribution (from Distributions.jl) takes a probability vector and returns a random integer according to those probabilities—this is how we simulate the Markov chain.

Third, we store all results in a DataFrame for easy analysis. The DataFrame constructor takes column name-value pairs and produces a tabular data structure that supports filtering, transformation, and statistical operations.

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	3.53288	1.0	1.51664	1.16335	0.353288
2	2	3.53288	1.0	1.51664	1.16335	0.353288
3	3	3.53288	1.0	1.51664	1.16335	0.353288
4	4	3.53288	1.0	1.51664	1.16335	0.353288
5	5	3.53288	1.0	1.51664	1.16335	0.353288
6	6	3.53288	1.0	1.51664	1.16335	0.353288

Each row is one period, and columns contain all the macro aggregates. The first(sim, 6) command displays the first six rows, letting us verify that the simulation starts at the steady state and begins responding to productivity shocks.

5.1.1 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")

The simulation reveals the signature patterns of business cycles. Consumption is visibly smoother than output—it fluctuates but with smaller amplitude. This is consumption smoothing in action: the household spreads income over time using savings as a buffer rather than consuming all of a good shock immediately. Investment is much more volatile than output—it swings widely above and below its mean. When output rises, most of the extra income goes to investment; when output falls, investment bears the brunt of the cut. Plotting the three series on separate panels (rather than overlaid) makes these amplitude differences easy to read off directly. This is exactly what we see in U.S. macro data: consumption moves gently while investment swings dramatically.

5.1.2 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 rises gradually during expansions and falls slowly during contractions. This is because capital is a stock variable: it takes many periods of high investment to significantly increase the capital stock, and many periods of low investment for it to depreciate away. This sluggishness is the model’s propagation mechanism: even if a productivity shock lasts only one period, the capital stock takes many periods to adjust back to steady state, generating the persistent dynamics that characterize business cycles. The title of the original paper, “Time to Build and Aggregate Fluctuations,” refers precisely to this insight: the time it takes to build and deplete capital is what turns transitory productivity shocks into sustained macroeconomic fluctuations.

5.2 Business Cycle Statistics

The visual impression from the plots is informative, but the real test of the model is quantitative: can it reproduce the specific statistical patterns observed in U.S. data? The standard approach is to compute business cycle statistics: the standard deviations and correlations of log deviations from trend for each macro aggregate. We simulate a long series (10,000 periods to minimize sampling noise) and compute these moments.

The transform! function from DataFrames.jl adds new columns computed from existing ones. Here we compute log deviations from mean: \(\hat{x}_t = \log x_t - \overline{\log x}\). These are percentage deviations from trend (since \(\log(1 + x) \approx x\) for small \(x\)), which is the standard way to measure business cycle fluctuations:

sim_long = simulateRBC(m, sol, 10_000)

transform!(sim_long,
    :y => (x -> log.(x) .- mean(log.(x))) => :y_dev,
    :c => (x -> log.(x) .- mean(log.(x))) => :c_dev,
    :i => (x -> log.(x) .- mean(log.(x))) => :i_dev,
    :k => (x -> log.(x) .- mean(log.(x))) => :k_dev)

10000×10 DataFrame

9975 rows omitted

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	3.53288	1.0	1.51664	1.16335	0.353288	-0.00353789	-0.00357166	-0.00165563	-0.00479484
2	2	3.53288	1.0	1.51664	1.16335	0.353288	-0.00353789	-0.00357166	-0.00165563	-0.00479483
3	3	3.53288	1.0	1.51664	1.16335	0.353288	-0.00353788	-0.00357166	-0.00165562	-0.00479483
4	4	3.53288	1.0	1.51664	1.16335	0.353288	-0.00353788	-0.00357166	-0.00165562	-0.00479483
5	5	3.53288	1.0	1.51664	1.16335	0.353288	-0.00353788	-0.00357166	-0.00165561	-0.00479482
6	6	3.53288	1.0	1.51664	1.16335	0.353288	-0.00353788	-0.00357166	-0.00165561	-0.00479482
7	7	3.53288	1.0	1.51664	1.16335	0.353288	-0.00353788	-0.00357166	-0.0016556	-0.00479481
8	8	3.53288	1.0	1.51664	1.16335	0.353288	-0.00353788	-0.00357166	-0.0016556	-0.00479481
9	9	3.53288	1.0	1.51664	1.16335	0.353288	-0.00353788	-0.00357166	-0.0016556	-0.0047948
10	10	3.53288	1.0	1.51664	1.16335	0.353288	-0.00353788	-0.00357166	-0.00165559	-0.0047948
11	11	3.53288	1.0	1.51664	1.16335	0.353288	-0.00353787	-0.00357166	-0.00165559	-0.0047948
12	12	3.53288	1.0	1.51664	1.16335	0.353288	-0.00353787	-0.00357166	-0.00165558	-0.00479479
13	13	3.53288	1.0	1.51664	1.16335	0.353288	-0.00353787	-0.00357166	-0.00165558	-0.00479479
⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮
9989	9989	3.56075	1.03817	1.57862	1.19233	0.386299	0.0365189	0.0210287	0.0876725	0.00306369
9990	9990	3.59098	1.0	1.52483	1.17318	0.351649	0.00184465	0.00483755	-0.00630495	0.0115159
9991	9991	3.58353	1.0	1.52378	1.17163	0.352151	0.00115946	0.00351832	-0.00487775	0.00943955
9992	9992	3.57733	0.96323	1.46691	1.14967	0.317239	-0.0368755	-0.015399	-0.109285	0.00770755
9993	9993	3.53683	0.96323	1.46141	1.14071	0.320704	-0.0406323	-0.0232284	-0.0984199	-0.00367664
9994	9994	3.50385	0.96323	1.4569	1.13552	0.321384	-0.0437238	-0.0277893	-0.0963025	-0.0130448
9995	9995	3.47485	0.96323	1.45291	1.13251	0.320403	-0.0464666	-0.030443	-0.0993602	-0.0213563
9996	9996	3.44777	0.96323	1.44916	1.12806	0.321107	-0.0490486	-0.0343806	-0.0971655	-0.0291806
9997	9997	3.4241	0.96323	1.44587	1.12328	0.322597	-0.051322	-0.0386278	-0.0925355	-0.0360697
9998	9998	3.40429	0.96323	1.44311	1.11926	0.323843	-0.0532371	-0.0422059	-0.0886813	-0.0418728
9999	9999	3.3877	0.96323	1.44078	1.11588	0.324902	-0.0548488	-0.0452337	-0.085414	-0.0467568
10000	10000	3.37383	0.96323	1.43883	1.11304	0.325798	-0.0562024	-0.0477864	-0.0826603	-0.0508587

σ_y = std(sim_long.y_dev)
σ_c = std(sim_long.c_dev)
σ_i = std(sim_long.i_dev)
σ_k = std(sim_long.k_dev)
stats = DataFrame(
    Variable   = ["Output", "Consumption",
                  "Investment", "Capital"],
    σ_pct      = 100 .* [σ_y, σ_c, σ_i, σ_k],
    σ_relative = [1.0, σ_c/σ_y, σ_i/σ_y, σ_k/σ_y],
    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.20502	1.0	1.0
2	Consumption	5.42276	0.873931	0.97211
3	Investment	10.2088	1.64525	0.910738
4	Capital	6.47982	1.04429	0.86667

ImportantKey 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)

TipThe 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.

5.3 Comparative Statics

Comparative statics, changing one parameter while holding the others fixed, helps us understand what drives the model’s predictions. We examine three parameters: the persistence of productivity shocks \(\rho\), capital’s share \(\alpha\), and the discount factor \(\beta\). For each experiment, we re-solve and re-simulate the model with a smaller grid (\(N_k = 100\)) for speed.

5.3.1 Effect of Persistence (\(\rho\))

The persistence parameter \(\rho\) controls how long productivity shocks last. How does this affect the character of business cycles?

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 produces longer, more pronounced business cycles. With \(\rho = 0.50\), output fluctuates rapidly but reverts to trend quickly: the economy forgets shocks fast, so booms and busts are short-lived. With \(\rho = 0.95\), a single positive shock sustains an expansion for many periods because the economy expects high productivity to continue. This comparison underscores why the calibrated value of \(\rho\) matters so much: the model needs \(\rho \approx 0.90\) to generate the persistent cycles we observe in the data.

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

Capital’s share \(\alpha\) determines how productive capital is relative to labor. A higher \(\alpha\) means capital contributes more to output, which affects both the steady state and the dynamics.

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)
    iz_mid = (m_trial.N_z + 1) ÷ 2
    plot!(m_trial.k_grid, sol_trial.kp[:, iz_mid],
          linewidth=2, label="α = $α_val")
end
xlabel!("Current Capital k")
ylabel!("Next Period Capital k'")

Higher \(\alpha\) makes capital more productive at the margin, which increases the return to saving. The planner responds by saving more and accumulating a higher steady-state capital stock. Notice that the policy functions shift upward with \(\alpha\) (the planner chooses a higher \(k'\) for any given \(k\)) and the crossing point with the 45° line moves right, confirming a higher steady state.

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

The discount factor \(\beta\) measures how much the household values the future relative to the present. A patient household (\(\beta\) close to 1) is willing to defer consumption in exchange for higher future income from capital accumulation.

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)
    iz_mid = (m_trial.N_z + 1) ÷ 2
    plot!(m_trial.k_grid, sol_trial.kp[:, iz_mid],
          linewidth=2, label="β = $β_val")
end
xlabel!("Current Capital k")
ylabel!("Next Period Capital k'")

More patient agents (\(\beta\) closer to 1) save more at every capital level and accumulate more capital in the long run. The policy function for \(\beta = 0.99\) lies well above the one for \(\beta = 0.90\), and the steady state shifts dramatically rightward. Impatient agents (\(\beta = 0.90\)) consume heavily today, save little, and arrive at a much lower steady state. The difference is quantitatively large: the steady-state capital stock roughly doubles between \(\beta = 0.90\) and \(\beta = 0.99\).

6 Linearized Solution

6.1 The Euler Equation

Value function iteration gives an accurate global solution (it works at any \((k, z)\) in the state space), but it is computationally expensive. Each Bellman iteration solves \(N_k \times N_z = 200 \times 7 = 1{,}400\) optimization problems, and convergence takes many iterations. For a model with one continuous state variable this is manageable, but as models grow larger (multiple capital stocks, heterogeneous agents, many shocks), VFI becomes prohibitively slow.

A linearized solution provides a fast, analytical alternative. Instead of solving the Bellman equation globally, we take a first-order Taylor expansion of the model’s optimality conditions around the deterministic steady state. This gives a linear policy rule that maps today’s state into tomorrow’s capital—no iteration, no optimization, no interpolation. The cost is accuracy: the linearized solution is exact at the steady state and increasingly approximate as we move away from it. But since business cycle fluctuations are small (a few percent of GDP), the linearized solution is typically excellent.

The approach is to linearize the first-order conditions around the steady state. 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 take first-order Taylor expansions around the steady state, dropping products of deviations like \(\hat{k}_t \cdot \hat{z}_t\) as second-order small.

TipWhy Linearize?

VFI is accurate but slow (\(N_k \times N_z\) optimization problems per iteration). Linearization gives a policy rule that is nearly exact near the steady state—and almost instant to compute. Most of the macro literature uses linearization or higher-order perturbation.

The linearized system is built from two equations: the Euler equation (the household’s optimality condition for savings) and the resource constraint. The Euler equation says that the marginal cost of saving one more unit today equals the expected discounted marginal benefit of the extra capital tomorrow:

\[ \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] \]

The left side, \(1/c_t\), is the marginal utility of consumption today—the cost of forgoing one unit of consumption. The right side is the discounted expected marginal utility of the return from investing that unit: each unit saved yields \(\alpha z_{t+1} k_{t+1}^{\alpha-1} + 1 - \delta\) units next period (the marginal product of capital plus undepreciated capital), valued at \(1/c_{t+1}\). At the optimum, the household is indifferent between consuming one more unit today and saving it to consume the return tomorrow. This is the same intertemporal tradeoff that generates the steady-state condition \(1 = \beta(\alpha k_{ss}^{\alpha-1} + 1 - \delta)\), but now applied to every period and every state of the world.

6.2 Linearizing the Resource Constraint

We now linearize the two equations around the steady state, starting with the simpler one: the aggregate resource constraint.

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

We write each variable as its steady-state value plus a deviation: \(k_t = k_{ss} + \hat{k}_t\), \(c_t = c_{ss} + \hat{c}_t\), \(z_t = 1 + \hat{z}_t\). The nonlinear term \(z_t k_t^\alpha\) is linearized using a first-order Taylor expansion: \((1 + \hat{z}_t)(k_{ss} + \hat{k}_t)^\alpha \approx k_{ss}^\alpha + \alpha k_{ss}^{\alpha-1} \hat{k}_t + k_{ss}^\alpha \hat{z}_t\), dropping products of deviations as second-order small. Substituting into the resource constraint and subtracting the steady-state equation (\(c_{ss} + k_{ss} = y_{ss} + (1-\delta)k_{ss}\)) to cancel the constants:

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

Rearranging:

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

This expresses \(\hat{c}_t\) as a linear function of \(\hat{k}_t\), \(\hat{k}_{t+1}\), and \(\hat{z}_t\), using the steady-state condition \(r_{ss} + 1 - \delta = 1/\beta\).

6.3 Linearizing the Euler Equation

The gross return to capital is \(R_{t+1} = \alpha z_{t+1} k_{t+1}^{\alpha-1} + 1 - \delta\), with steady-state value \(R_{ss} = 1/\beta\). Linearizing the Euler equation \(1/c_t = \beta \mathbb{E}_t[R_{t+1}/c_{t+1}]\) around the steady state and using \(\beta R_{ss} = 1\) to cancel the constants:

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

The linearized return is:

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

More capital lowers the return (diminishing returns, \(\alpha - 1 < 0\)), while higher productivity raises it.

6.4 Solving the Linearized System

We now have two linearized equations (resource constraint and Euler equation) in three unknowns (\(\hat{c}_t\), \(\hat{k}_{t+1}\), and exogenous \(\hat{z}_t\)). The method of undetermined coefficients gives us the solution.

6.4.1 Guessing a Linear Policy Rule

Since both equations and the exogenous process \(\hat{z}_{t+1} = \rho \hat{z}_t + \sigma_\varepsilon \varepsilon_{t+1}\) are linear, we guess that the optimal policy rule is also linear:

\[ \hat{k}_{t+1} = \phi_k \, \hat{k}_t + \phi_z \, \hat{z}_t \]

Substituting into the linearized resource constraint:

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

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

6.4.2 Finding \(\phi_k\)

The crucial step is to substitute our guess into the linearized Euler equation and require that the equation holds for all values of \((\hat{k}_t, \hat{z}_t)\). Since both sides are linear in \(\hat{k}_t\) and \(\hat{z}_t\), this means the coefficients on \(\hat{k}_t\) must match on both sides, and separately the coefficients on \(\hat{z}_t\) must match. 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\), substituting everything into the Euler equation, and matching coefficients on \(\hat{k}_t\):

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

Since \(\psi_k = 1/\beta - \phi_k\), this is one equation in one unknown:

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

This is a quadratic equation in \(\phi_k\), so it has two roots. We pick the stable root (\(|\phi_k| < 1\)), which ensures that capital mean-reverts toward the steady state after a shock. The unstable root (\(|\phi_k| > 1\)) would imply that deviations from steady state grow without bound, which is economically nonsensical for a stationary model.

6.4.3 Finding \(\phi_z\)

Matching coefficients on \(\hat{z}_t\) gives a linear equation for \(\phi_z\):

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

This is arithmetic once \(\phi_k\) (and hence \(\psi_k\)) is known.

ImportantSummary 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 resource constraint
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\)

6.5 Julia Implementation

The implementation is short because linearization reduces the problem to solving one nonlinear equation (the quadratic for \(\phi_k\)) and one linear equation (for \(\phi_z\)). We use NonlinearSolve.jl with Newton’s method for the quadratic, initializing at \(\phi_k = 0.9\) (close to the stable root we expect).

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
    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
    num = y_ss*(1 - ρ) + β*c_ss*r_ss*ρ
    den = ψ_k - β*c_ss*γ + (1 - ρ)
    ϕ_z = num / den

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

linearize_rbc (generic function with 1 method)

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

ϕ_k = 0.8589
ϕ_z = 0.9403

The linearized capital policy rule is:

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

The coefficient \(\phi_k\) tells us how persistent capital deviations are: if capital is 1% above steady state today, it will be roughly \(\phi_k \times 1\%\) above steady state tomorrow (holding productivity constant). Since \(|\phi_k| < 1\), capital mean-reverts: deviations shrink over time as the economy returns to steady state. The coefficient \(\phi_z\) tells us how much capital responds to a productivity shock: since \(\phi_z > 0\), higher productivity leads to more saving. A 1% productivity shock increases next period’s capital by roughly \(\phi_z \times 1\%\).

6.6 Comparing to VFI

How good is the linearized approximation? We compare the linearized policy functions and simulated time series directly against the VFI solution. If the two are close, linearization gives us the same economic insights at a fraction of the computational cost.

6.6.1 Policy Functions

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

plot(m.k_grid, sol.kp[:, iz_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: the linearized policy function is almost indistinguishable from the VFI solution. Far from the steady state, the linearized solution diverges because it is a straight line that cannot capture the curvature of the true policy function. The VFI solution curves gently (reflecting diminishing returns), while the linearized solution continues in a straight line. For typical business cycle fluctuations, which stay within \(\pm 10\%\) of the steady state, this divergence is negligible.

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

6.6.2 Simulating the Linearized Model

The linearized simulation function is nearly identical to the VFI simulation, except the policy function is the simple linear rule \(k_{t+1} = k_{ss} + \phi_k(k_t - k_{ss}) + \phi_z(z_t - 1)\) instead of an interpolated object. We use the same random seed so both simulations face identical productivity shocks, making the comparison apples-to-apples:

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);  z_idx = zeros(Int, T)
    y_sim = zeros(T);  c_sim = zeros(T);  i_sim = zeros(T)
    k_sim[1] = k_ss
    z_idx[1] = (N_z + 1) ÷ 2
    for t in 1:T
        z = z_grid[z_idx[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
            z_idx[t+1] = rand(Categorical(P_z[z_idx[t], :]))
        end
    end
    return DataFrame(t=1:T, k=k_sim,
        z=z_grid[z_idx], y=y_sim, c=c_sim, i=i_sim)
end

simulateRBC_linear (generic function with 1 method)

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")

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))

The two simulations track each other very closely—the solid and dashed lines are nearly superimposed throughout the 200-period simulation. This confirms that for the moderate fluctuations generated by the calibrated productivity process, linearization is an excellent approximation to the full VFI solution. The small differences that do appear tend to emerge during the most extreme episodes, where the economy moves furthest from the steady state and the linear approximation becomes less accurate.

ImportantLinearization: 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.

7 Summary

These notes developed the Kydland & Prescott (1982) Real Business Cycle model from first principles and solved it using two complementary methods: value function iteration with linear interpolation, and linearization of the first-order conditions.

Economic Concepts:

• The RBC model shows that a competitive economy with technology shocks can generate realistic business cycle fluctuations without any market failures.
• The First Welfare Theorem lets us replace the competitive equilibrium with the simpler social planner’s problem.
• The planner’s problem is a stochastic Bellman equation with state variables \((k, z)\) and control \(k'\).
• Calibration sets parameters to match long-run data moments (\(\alpha\), \(\beta\), \(\delta\)) and the estimated Solow residual process (\(\rho\), \(\sigma_\varepsilon\)).
• The model reproduces key stylized facts: consumption smoothing (\(\sigma_c / \sigma_y < 1\)), volatile investment (\(\sigma_i / \sigma_y > 1\)), and procyclical comovement.
• The linearized solution approximates the Euler equation around the steady state and delivers a closed-form policy rule \(\hat{k}_{t+1} = \phi_k \hat{k}_t + \phi_z \hat{z}_t\) that is highly accurate near the steady state.

Julia Tools:

Key Julia tools used in these notes

Tool	Purpose
Interpolations.jl: linear_interpolation	Evaluate \(V(k', z)\) at off-grid \(k'\)
Optim.jl: Brent()	1D optimization of the Bellman equation
NonlinearSolve.jl: NonlinearProblem, NewtonRaphson()	Solve for the linearized coefficient \(\phi_k\)
DataFrames.jl: DataFrame, transform!	Tabular simulation output and in-place column operations
@kwdef mutable struct + setup!	Bundle parameters; rebuild grids consistently
(; ...) = m destructuring	Extract struct fields cleanly
Categorical	Draw next state from Markov transition probabilities
Broadcasting .	Element-wise operations on policy arrays