Solving DSGE Models Numerically: Perturbation and Global Methods

A dynamic stochastic general equilibrium model is not "solved" by evaluating its equations at a few future values. The numerical task is to recover functions that map the state of the economy into decisions. If the state is

a solution typically consists of policy functions such as

together with laws of motion for the exogenous states. That distinction matters because a piece of code can satisfy one Euler equation at one pair of points without solving the dynamic model.

The earlier version of this article made exactly that mistake. It called a root finder on a static expression "first-order perturbation" and called a forward numerical derivative "finite-difference solution." Neither operation computes the policy functions of a DSGE model.

This revision starts from the actual numerical problem.

A benchmark stochastic growth model

Consider a planner with preferences

subject to

with productivity following a Markov process

The state is

The control can be written as next-period capital $k_{t+1}$, because consumption is determined by the resource constraint:

A solution is a policy

that satisfies optimality for every relevant state, not merely along one simulated path.

The Bellman equation

The recursive problem is

where

The expectation is

This equation gives us one route to a global numerical solution: value-function iteration.

The Euler equation

For an interior optimum, the policy also satisfies

A correct numerical solution should make the Euler-equation residual small over the region where the policy will be used. That residual is a diagnostic. It is not, by itself, a solution algorithm.

Deterministic steady state

At

and a non-stochastic steady state,

so the Euler equation becomes

Therefore,

Consumption is

The steady state is important for both local and global algorithms. Perturbation expands the solution around it. A global method often uses it to define a sensible computational domain.

What perturbation actually means

Let

collect the equilibrium conditions, where $\sigma$ scales shock size. A perturbation method treats the equilibrium policy function as an unknown smooth function of the state and shock scale and computes derivatives of that function around the deterministic steady state,

A first-order approximation has the form

where hats denote deviations from steady state, often in logs. The matrices $A$ and $B$ are not obtained by calling a generic scalar root finder on the Euler equation. They come from differentiating the complete equilibrium system and solving the resulting linear rational-expectations problem. At first order, certainty equivalence commonly appears: shock variances do not change the mean policy rule.

At second order, curvature introduces terms involving variances and interactions, allowing uncertainty to affect expected decisions and welfare. That is one reason second-order perturbation is used for risk premia and welfare calculations.

Local accuracy is the main trade-off

Perturbation is attractive because it is fast. For large macroeconomic models with many state and control variables, a first- or second-order local solution can be dramatically cheaper than constructing a high-dimensional global grid. But the approximation is local. If the economy moves far from the expansion point, or if occasionally binding constraints matter, the truncated Taylor expansion can become inaccurate or even imply impossible decisions.

This is not a defect in Taylor series. It is the consequence of asking a local approximation to describe a global nonlinear problem.

Finite differences are not a competing DSGE solution method by themselves

A finite difference such as

approximates a derivative. That tool can appear inside many numerical algorithms. For example, finite differences can approximate derivatives in a Hamilton-Jacobi-Bellman equation, compute Jacobians for Newton methods, or discretize a continuous-state problem. But evaluating finite differences of the production function does not solve the DSGE model.

The numerical method is defined by the equation being discretized and the policy or value function being recovered. So the meaningful comparison is not

It is closer to

with finite differences being one possible numerical ingredient.

A global alternative: value-function iteration

For the growth model above, discretize capital on a grid

and productivity on states

For each current state $(k_i,z_s)$ and each candidate $k_j'$, compute feasible consumption

The Bellman update is

Repeat until

The maximizing index at each state is the discrete policy function.

Reproducible Python implementation

The following code solves a two-state stochastic growth model by value-function iteration.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
from __future__ import annotations

import numpy as np
from numpy.typing import NDArray

FloatArray = NDArray[np.float64]
IntArray = NDArray[np.int64]

beta: float = 0.96
alpha: float = 0.36
delta: float = 0.08

productivity: FloatArray = np.exp(
    np.array([-0.02, 0.02], dtype=float)
)

transition: FloatArray = np.array(
    [
        [0.95, 0.05],
        [0.05, 0.95],
    ],
    dtype=float,
)

if transition.shape != (2, 2):
    raise ValueError("Transition matrix must be 2 x 2.")

if not np.allclose(
    transition.sum(axis=1),
    1.0,
):
    raise ValueError(
        "Each transition row must sum to one."
    )

k_ss: float = (
    alpha
    / (1.0 / beta - 1.0 + delta)
) ** (1.0 / (1.0 - alpha))

capital_grid: FloatArray = np.linspace(
    0.5 * k_ss,
    1.5 * k_ss,
    250,
)

value: FloatArray = np.zeros(
    (
        productivity.size,
        capital_grid.size,
    ),
    dtype=float,
)

policy_index: IntArray = np.zeros(
    value.shape,
    dtype=np.int64,
)

tolerance: float = 1e-8
max_iterations: int = 1_000

for iteration in range(max_iterations):
    value_new = np.empty_like(value)
    policy_new = np.empty_like(policy_index)

for z_index, z_value in enumerate(
        productivity
    ):
        expected_value: FloatArray = (
            transition[z_index] @ value
        )

for k_index, capital in enumerate(
            capital_grid
        ):
            resources: float = (
                z_value * capital**alpha
                + (1.0 - delta) * capital
            )

consumption: FloatArray = (
                resources - capital_grid
            )

objective: FloatArray = np.full(
                capital_grid.shape,
                -np.inf,
                dtype=float,
            )

feasible: NDArray[np.bool_] = (
                consumption > 0.0
            )

objective[feasible] = (
                np.log(consumption[feasible])
                + beta
                * expected_value[feasible]
            )

best_index: int = int(
                np.argmax(objective)
            )

value_new[
                z_index,
                k_index,
            ] = objective[best_index]

policy_new[
                z_index,
                k_index,
            ] = best_index

sup_norm: float = float(
        np.max(np.abs(value_new - value))
    )

value = value_new
    policy_index = policy_new

if sup_norm < tolerance:
        break
else:
    raise RuntimeError(
        "Value iteration did not converge."
    )

policy_capital: FloatArray = (
    capital_grid[policy_index]
)

print(
    "iterations:",
    iteration + 1,
)

print(
    "sup-norm:",
    f"{sup_norm:.3e}",
)

print(
    "steady-state capital:",
    f"{k_ss:.6f}",
)

With the parameters above, the code converges on the specified grid. The policy moves next-period capital upward in the high-productivity state and downward in the low-productivity state around the deterministic steady state, which is the direction economic intuition predicts.

Grid error is different from model error

A value-function solution on a finite grid contains discretization error. If the true optimum lies between $k_j$ and $k_{j+1}$, a discrete policy must choose one grid point. Increasing the grid density reduces this source of error but increases computational cost. Interpolation can improve the approximation without making the grid prohibitively dense.

The distinction is worth making explicit:

Perturbation mainly introduces truncation error from a local Taylor expansion. Grid methods introduce discretization and interpolation error. Both need diagnostics.

Euler-equation errors are a useful common diagnostic

After obtaining a policy $g(k,z)$, compute consumption from the resource constraint and evaluate the Euler residual,

A small residual over the relevant state space indicates that the approximate policy nearly satisfies the first-order condition. This allows different numerical methods to be compared on a common economic equation rather than on implementation-specific convergence criteria alone.

Projection and collocation methods

Value-function iteration is not the only global method. Projection methods approximate an unknown policy or value function by basis functions,

then choose the coefficients $a_m$ so that equilibrium residuals are small at selected collocation points or in a weighted integral sense. Chebyshev polynomials are a common basis because they have good approximation properties over bounded intervals. Projection can be much faster than dense grids in smooth low-dimensional problems. The curse of dimensionality remains important.

Occasionally binding constraints change the method choice

Consider a borrowing constraint,

or a policy rate constrained by

Near a point where the constraint never binds, a local perturbation can completely miss the kink created when it becomes active. Piecewise-linear methods, occasionally binding constraint algorithms, endogenous-grid methods, projection, or other global approaches can be more appropriate. The numerical method should follow the economic structure.

Choosing a method

A useful summary is:

There is no universally best solver. There is a model, a region of the state space that matters, and an accuracy requirement.

Conclusion

A DSGE solution is a set of decision rules satisfying equilibrium conditions over the relevant states. Perturbation obtains local derivatives of those rules around a steady state. Global methods approximate the functions over a larger region. Finite differences can help approximate derivatives inside such methods, but a finite-difference formula is not itself a DSGE solution.

That distinction is the difference between numerical analysis and code that merely produces numbers.

References

• Schmitt-Grohé, S., & Uribe, M. (2004). Solving dynamic general equilibrium models using a second-order approximation to the policy function. Journal of Economic Dynamics and Control, 28(4), 755–775.
• Judd, K. L. (1998). Numerical Methods in Economics. MIT Press.
• Miranda, M. J., & Fackler, P. L. (2002). Applied Computational Economics and Finance. MIT Press.
• Heer, B., & Maussner, A. (2009). Dynamic General Equilibrium Modeling (2nd ed.). Springer.