The BKM Algorithm

Finding solutions to economic models with heterogeneity and aggregate uncertainty is a notoriously difficult task. It is an active area of research in Mathematics (see mean field games with aggregate uncertainty). Because such models naturally arise when considering economic situations, Economists have developed a battery of techniques to (numerically) solve them. In this post, I would like to describe a recent algorithm by Boppart, Krusell and Mitman (BKM) that is both fast and accurate. I will first describe the problem that Economists face when working with heterogeneous model with aggregate uncertainty, heuristically discuss the BKM algorithm (based on the presentation of Reiter (2018)) and show an application in Julia.

BKM in a nutshell

It is common to use Bellman’s principle of optimality to characterize solutions of a multi-stage decision processes. The principle of optimality leads to a solution in a recursive form $d_{t} = d(S_{t})$, where $S_t$ is a vector of state variables and $d(.)$ a policy function describing the optimal action of a decision-maker when faced with any given state.

An alternative representation of the problem is to consider a solution in the sequence form. At each step in time, a decision-maker observes a new realization of a random process $z_t$ and taking into account the full history of past shocks, the agent makes a new choice to maximize her expected discounted sum of returns: $d_t = d(z_t, z_{t-1}, z_{t-2}, …)$.

While most of the time the recursive form is a much more parsimonious approach, it fails when $S_t$ is infinite-dimensional. In models with heterogeneous agents (HA) and aggregate uncertainty, this is generally the case because the distribution of agents over certain variables will end up being in $S_t$. While this a problem with is the recursive approach (how can we discretize $S_t$ to put it on a computer?), the sequence form is immune to this problem. The BKM algorithm uses this insight, adding the assumption of linearity of $d(.)$ with respect to the aggregate state $z_t$:

$$ d_t = z_t d(1, 0, 0, …) + z_{t-1}d(0, 1, 0, …) + z_{t-2}d(0, 0, 1, …) + … $$ $$ d_t = \sum_{k=0}^{+\infty} z_{t-k} d_{k} $$

where $$ d_{1} = d(1,0,0,…)$$ $$ d_{2} = d(0,1,0,…)$$ $$ d_{3} = d(0,0,1,…)$$

The series of $d_{k}$ describes the trajectory of the economy after an “MIT” shock (when the economy is hit by an unexpected single-period aggregate shock). If the linearity assumption holds, the evolution of equilibrium variables are simply a moving average of past shocks. We made progress because solving for the trajectory after an MIT shock is perfectly feasible in a RA model. As long as one can calculate the steady-state with no aggregate uncertainty and solve the deterministic perfect foresight path, BKM can be used.

Numerical implementation

Here is the implementation of the method using a toy example. I intentionally circumvent the problem of finding the perfect foresight transition path, which is (potentially) the complicated part of BKM. This is the example given in Reiter (2018): the exact model is the non-linear model $x_{t} = a x_{t-1} + b x_{t-1}^2 + z_{t}$

using Distributions
using Plots
# Parameters
max_iter=100 #number of iterations for the simulation
a = 0.5
b = 0.1
sigma_shock=1.0
mu_shock=0.
d = Normal(mu_shock, sigma_shock)

# transition function
function iter_x(x_min1::Float64, a::Float64, b::Float64)
    a*x_min1 + b*x_min1^2
end
iter_x (generic function with 1 method)

Impulse response function

Let us assume that the economy is at the non-stochastic steady-state and a shock occurs:

# We assume that after 100 periods, the economy is back to the steady-state
max_iter_mit = 100
x_mit=zeros(max_iter_mit)
# Initial shock
z_t[1] = sigma_shock #a 1 std. deviation
x_mit[1] = 0 #steady-state

for i=2:max_iter_mit
    x_mit[i] = iter_x(x_mit[i-1], a, b) + z_t[i-1]
end

# Scaled-version of the impulse response:
x_mit_scaled = x_mit./z_t[1];

I define two functions to calculate the moving average:

function calculate_Xt(x_scaled::Array{Float64,1}, shocks::Array{Float64,1}, t::Int64, kk::Int64)

    output = 0.

    for k=1:kk
        output+=x_scaled[k]*shocks[t-k+1]
    end

    return output
end

function BPM_path!(XT::Array{Float64,1}, max_iter::Int64, x_scaled::Array{Float64,1}, shocks::Array{Float64,1})

    for t=2:max_iter
        XT[t] = calculate_Xt(x_scaled, shocks, t, t)
    end

end
BPM_path! (generic function with 1 method)
# Initialization
XT = zeros(max_iter);
# Series of shocks
shocks_t = rand(d, max_iter).*0.5;
# Solving using BKM:
@time BPM_path!(XT, max_iter, x_mit_scaled, shocks_t)
# True value:
x_true = zeros(max_iter)
for i=2:max_iter
    x_true[i] = iter_x(x_true[i-1], a, b) + shocks_t[i-1]
end
  0.025906 seconds (29.30 k allocations: 1.489 MiB, 26.80% gc time)
plot(XT[2:end], label="x_t BKM")
plot!(x_true[2:end], label="x_t true")
title!("b=$(b)")

png

The previous plot shows how the BKM algorithm approximates the true model and it does quite a good job. Of course, the more the model is linear with respect to $z_t$ (captured by the value of $b_2$), the better the approximation. To illustrate this idea, I use BKM on a perfectly linear model ($b=0$) and on a model with stronger non-linearities ($b=0.2$). As expected, the approximation is perfect when the model is linear and the approximation deteriorates when strong non-linearities are present.

plot(p1,p2)

png

Conclusion

The BKM algorithm is a new addition to the toolbox of methods to solve HA models. The sequence representation of the problem seems to be a fruitful area of research, as it has also been used in Le Grand and Ragot (2017) and (2019) to develop fast and reliable method to solve HA models. The recent contribution of Auclert at al (2019) also uses a similar approach.

References

  • Achdou, Yves, et al. “Partial differential equation models in macroeconomics.” Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 372.2028 (2014): 20130397.
  • Auclert, Adrien, et al. Using the Sequence-Space Jacobian to Solve and Estimate Heterogeneous-Agent Models. No. w26123. National Bureau of Economic Research, 2019.
  • Boppart, Timo, Per Krusell, and Kurt Mitman. “Exploiting MIT shocks in heterogeneous-agent economies: the impulse response as a numerical derivative.” Journal of Economic Dynamics and Control 89 (2018): 68-92.
  • Le Grand, François, and Xavier Ragot. “Optimal fiscal policy with heterogeneous agents and aggregate shocks.” Document de travail (2017).
  • Le Grand, François, and Ragot, Xavier. “Managing Inequality over the Business Cycles: Optimal Policies with Heterogeneous Agents and Aggregate Shocks”. No. 1090. Society for Economic Dynamics, 2019.
  • Reiter, Michael. “Comments on” Exploiting MIT Shocks in Heterogeneous-Agent Economies: The Impulse Response as a Numerical Derivative” by T. Boppart, P. Krusell and K. Mitman.” Journal of Economic Dynamics and Control 89 (2018): 93-99.

Appendix

versioninfo()
Julia Version 1.0.3
Commit 099e826241 (2018-12-18 01:34 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Core(TM) i7-8850H CPU @ 2.60GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.0 (ORCJIT, skylake)
### To create plot 2

#### Linear Model

# Parameters
max_iter=100 #number of iterations for the simulation
a = 0.5
b = 0.0
sigma_shock=1.0
mu_shock=0.
d = Normal(mu_shock, sigma_shock)

# IRF
x_mit=zeros(max_iter_mit)
# Initial shock
z_t[1] = sigma_shock #a 1 std. deviation
x_mit[1] = 0 #steady-state

for i=2:max_iter_mit
    x_mit[i] = iter_x(x_mit[i-1], a, b) + z_t[i-1]
end

# Scaled-version of the impulse response:
x_mit_scaled = x_mit./z_t[1];


# Initialization
XT = zeros(max_iter);
# Series of shocks
shocks_t = rand(d, max_iter).*0.5;
# Solving using BKM:
@time BPM_path!(XT, max_iter, x_mit_scaled, shocks_t)
# True value:
x_true = zeros(max_iter)
for i=2:max_iter
    x_true[i] = iter_x(x_true[i-1], a, b) + shocks_t[i-1]
end

p1 = plot(XT[2:end], label="x_t BKM")
plot!(x_true[2:end], label="x_t true")
title!("b=$(b)")

0.000009 seconds (4 allocations: 160 bytes)

### Plot 2

#### Linear Model

# Parameters
max_iter=100 #number of iterations for the simulation
a = 0.5
b = 0.2
sigma_shock=1.0
mu_shock=0.
d = Normal(mu_shock, sigma_shock)

# IRF
x_mit=zeros(max_iter_mit)
# Initial shock
z_t[1] = sigma_shock #a 1 std. deviation
x_mit[1] = 0 #steady-state

for i=2:max_iter_mit
    x_mit[i] = iter_x(x_mit[i-1], a, b) + z_t[i-1]
end

# Scaled-version of the impulse response:
x_mit_scaled = x_mit./z_t[1];


# Initialization
XT = zeros(max_iter);
# Series of shocks
shocks_t = rand(d, max_iter).*0.5;
# Solving using BKM:
@time BPM_path!(XT, max_iter, x_mit_scaled, shocks_t)
# True value:
x_true = zeros(max_iter)
for i=2:max_iter
    x_true[i] = iter_x(x_true[i-1], a, b) + shocks_t[i-1]
end

p2 = plot(XT[2:end], label="x_t BKM")
plot!(x_true[2:end], label="x_t true")
title!("b=$(b)")
  0.000007 seconds (4 allocations: 160 bytes)

png