The standard incomplete markets model: steady state

NBER Heterogeneous-Agent Macro Workshop

Matthew Rognlie

Spring 2022

Standard incomplete markets model

  • Come in with assets including return $(1+r)a$ and stochastic income $y(s)$
  • Allocate to assets for tomorrow $a'$ and consumption today $c$
  • ... all subject to a borrowing constraint $a'\geq \underline{a}$ and no ability to insure against the shock $s$
$$ \begin{gather} V(s,a)=\max_{c,a'} u(c)+\beta\mathbb{E}[V(s',a')|s] \tag{1}\\ \text{s.t. }a'+c=(1+r)a+y(s) \tag{2}\\ a'\geq \underline{a} \tag{3} \\ \end{gather} $$
  • (aka Bewley-Huggett-Aiyagari-Imrohoroglu-Zeldes-Deaton-Carroll model)

Nice features of this model

  • Simple way to think about consumption-saving problem given risk (income smoothing, precautionary savings, etc.)
  • Consistent with high marginal propensities to consume (MPCs) out of unanticipated one-time shocks, in line with data (crucial for us!)
  • Endogenously generates wealth distribution (and inequality)
  • Unlike representative-agent model, steady-state asset demand not infinitely elastic in $r$
  • Easy to extend
    • other shocks, preference heterogeneity, endogenous labor, life-cycle structure, other assets...
  • Backbone of everything we'll do!

Solving the model in steady state

Four steps, middle two the hardest:

  1. Discretizing state space: get Markov chain for $s$ and grid for $a$
  1. Solving for steady-state policy function: get policies $a'(s,a)$ and $c(s,a)$ on grid
  1. Solving for steady-state distribution: get discretized distribution $D(s,a)$ of agents on grid
  1. Aggregating: combine steps 2 and 3 to get aggregate assets and consumption

Preliminaries

Before we start, we need to import the Python libraries that we'll be using:

In [1]:
import numpy as np
import numba
from scipy import optimize
import matplotlib.pyplot as plt

# some useful plot defaults
plt.rcParams.update({'font.size' : 20, 'lines.linewidth' : 3.5, 'figure.figsize' : (13,7)})

Part 1 of computation: discretizing state space

Part 1A: discretizing the asset space

  • In theory, households pick assets from $[\underline{a},\infty)$, need to discretize into a grid
  • Will use linear interpolation on policy function, so we want more gridpoints where policy function is more nonlinear in assets
    • This turns out to be near $\underline{a}$
  • But also want grid going up to very high level of assets, to capture wealthiest
  • Solution: double-exponential transformation of a uniformly spaced grid
  • Start with uniformly spaced points $u_i$ on an interval $[0,\bar{u}]$, construct gridpoints $a_i$
$$ a_i = \underline{a} + e^{e^{u_i}-1}-1 \tag{4} $$

Discretizing asset space continued

Uniform points $u_i$ on $[0,\bar{u}]$, construct $a_i$ as

$$ a_i = \underline{a} + e^{e^{u_i}-1}-1 \tag{4} $$

This gives min of $\underline{a}$, max of $e^{e^{\bar{u}}-1}-1$. Equating this max with $\overline{a}$ gives

$$ \bar{u} = \log(1+\log(1+ \overline{a} - \underline{a})) $$

Implement in function below:

In [2]:
def discretize_assets(amin, amax, n_a):
    # find maximum ubar of uniform grid corresponding to desired maximum amax of asset grid
    ubar = np.log(1 + np.log(1 + amax - amin))
    
    # make uniform grid
    u_grid = np.linspace(0, ubar, n_a)
    
    # double-exponentiate uniform grid and add amin to get grid from amin to amax
    return amin + np.exp(np.exp(u_grid) - 1) - 1

Visualizing results

Example of 50 points between 0 and 10,000. (Normally we want more than this, but easier to visualize.)

In [3]:
a_grid = discretize_assets(0, 10_000, 50)
with np.printoptions(precision=2):
    print(a_grid)

[0.00e+00 4.98e-02 1.05e-01 1.65e-01 2.32e-01 3.07e-01 3.90e-01 4.82e-01
 5.86e-01 7.03e-01 8.34e-01 9.83e-01 1.15e+00 1.34e+00 1.57e+00 1.82e+00
 2.11e+00 2.45e+00 2.85e+00 3.31e+00 3.86e+00 4.51e+00 5.29e+00 6.21e+00
 7.34e+00 8.70e+00 1.04e+01 1.24e+01 1.50e+01 1.82e+01 2.23e+01 2.75e+01
 3.42e+01 4.29e+01 5.43e+01 6.96e+01 9.01e+01 1.18e+02 1.57e+02 2.11e+02
 2.87e+02 3.97e+02 5.58e+02 7.97e+02 1.16e+03 1.71e+03 2.58e+03 3.97e+03
 6.23e+03 1.00e+04]

Note vastly greater concentration near the borrowing constraint of 0: out of 50, first 12 points are less than 1, next 10 points are less than 5

Many variations on this (e.g. could add scalar parameters to double exponential, or try triple exponential as suggested by Carroll)

Part 1B: discretizing incomes

  • We had arbitrary Markov process $s$ and mapping to incomes $y(s)$
  • Lots of possible choices: could have two states for employed vs. unemployed, etc.
  • Simplest case: assume $\log y_t$ follows an AR(1) process with some persistence $\rho$ and normal innovations $\epsilon_t$
$$ \log y_t = \rho \log y_{t-1} + \epsilon_t \tag{5} $$
  • Try to approximate (5) with $n$ discrete states
  • Basis for Rouwenhorst method: models sum $s$ of $n-1$ independent hidden binary states (can be 0 or 1), each of which have probability $p$ of staying at current value, $1-p$ of switching
  • $n$ possible values of this sum, from 0 to $n-1$; as $n$ gets higher, sum $s$ behaves more and more like AR(1) with normal innovations

Getting Markov matrix (for completeness, won't cover in detail)

Can show that the Markov transition matrix $\Pi_n$ over values of $s$ obeys the recursion

$$ \Pi_{n}=p\begin{bmatrix}\Pi_{n-1} & \mathbf{0}\\ \mathbf{0}' & 0 \end{bmatrix}+(1-p)\begin{bmatrix}\mathbf{0} & \Pi_{n-1}\\ 0 & \mathbf{0}' \end{bmatrix}+(1-p)\begin{bmatrix}\mathbf{0}' & 0\\ \Pi_{n-1} & \mathbf{0} \end{bmatrix}+p\begin{bmatrix}0 & \mathbf{0}'\\ \mathbf{0} & \Pi_{n-1} \end{bmatrix} \tag{6} $$

with base case

$$ \Pi_{2}=\begin{bmatrix}p & 1-p\\ 1-p & p \end{bmatrix} \tag{7} $$

Implement:

In [4]:
def rouwenhorst_Pi(N, p):
    # base case Pi_2
    Pi = np.array([[p, 1 - p],
                   [1 - p, p]])
    
    # recursion to build up from Pi_2 to Pi_N
    for n in range(3, N + 1):
        Pi_old = Pi
        Pi = np.zeros((n, n))
        
        Pi[:-1, :-1] += p * Pi_old
        Pi[:-1, 1:] += (1 - p) * Pi_old
        Pi[1:, :-1] += (1 - p) * Pi_old
        Pi[1:, 1:] += p * Pi_old
        Pi[1:-1, :] /= 2
        
    return Pi

Using this to match an AR(1)

  • Drawing a random sequence $s_t$ using $\Pi_n$, one can show "persistence" $\rho$ from regressing $s_t$ on $s_{t-1}$ is $2p-1$, so to match desired persistence $\rho$ set
$$ p = \frac{1+\rho}{2} $$
  • $s_t$ is sum of $n-1$ independent states with 50-50 probability of being 0 or 1, so its stationary distribution is the binomial distribution with $n-1$ trials and success probability 1/2, which has standard deviation $\sigma = \sqrt{n-1}/2$
  • So to match desired cross-sectional standard deviation of log income, scale $s_t$ by constant
$$ \alpha \equiv \frac{2 \sigma}{\sqrt{n-1}} $$
  • Could add $C$ constant to $s_t$, but this will scale income by $e^C$, we'll pick scaling factor directly to normalize income

Technical preliminary: find stationary distribution of Markov chain

Given transition matrix $\Pi$, simple and robust way to find stationary distribution is to repeatedly apply transition matrix, updating distribution $\pi$ to $\Pi'\pi$, until convergence up to some tolerance:

In [5]:
def stationary_markov(Pi, tol=1E-14):
    # start with uniform distribution over all states
    n = Pi.shape[0]
    pi = np.full(n, 1/n)
    
    # update distribution using Pi until successive iterations differ by less than tol
    for _ in range(10_000):
        pi_new = Pi.T @ pi
        if np.max(np.abs(pi_new - pi)) < tol:
            return pi_new
        pi = pi_new

Now use our insights to obtain discretized income process

Given a target for persistence $\rho$, cross-sectional standard deviation $\sigma$, and the number of discrete states $n_s$, this function gives us the $y(s)$ for each state $s=0,\ldots,n_s-1$ (normalized to have mean of 1), the stationary distribution $\pi(s)$ over states, and Markov transition matrix $\Pi(s,s')$ between states.

In [6]:
def discretize_income(rho, sigma, n_s):
    # choose inner-switching probability p to match persistence rho
    p = (1+rho)/2
    
    # start with states from 0 to n_s-1, scale by alpha to match standard deviation sigma
    s = np.arange(n_s)
    alpha = 2*sigma/np.sqrt(n_s-1)
    s = alpha*s
    
    # obtain Markov transition matrix Pi and its stationary distribution
    Pi = rouwenhorst_Pi(n_s, p)
    pi = stationary_markov(Pi)
    
    # s is log income, get income y and scale so that mean is 1
    y = np.exp(s)
    y /= np.vdot(pi, y)
    
    return y, pi, Pi

Example

Let's try $\rho=0.975$, $\sigma=0.7$, and $n_s=7$. The former two values roughly match some estimates of the quarterly persistence, and the standard deviation, of the persistent component of log income in the US.

In [7]:
y, pi, Pi = discretize_income(0.975, 0.7, 7)
y, pi
Out[7]:
(array([0.1413694 , 0.25036602, 0.44339966, 0.78526334, 1.3907059 ,
        2.46294815, 4.36189534]),
 array([0.015625, 0.09375 , 0.234375, 0.3125  , 0.234375, 0.09375 ,
        0.015625]))

Verify mean of $y$ is 1 as desired:

In [8]:
np.vdot(pi, y)
Out[8]:
1.0

Verify standard deviation is $\sigma=0.7$:

In [9]:
mean_log_y = np.vdot(pi, np.log(y))
sd_log_y = np.sqrt(np.vdot(pi, (np.log(y) - mean_log_y)**2))
sd_log_y
Out[9]:
0.7000000000003479

Part 2 of computation: backward iteration to obtain the policy function

Back to standard incomplete markets model, with $t$ subscripts added:

$$ \begin{gather} V_t(s,a)=\max_{c,a'} u(c)+\beta\mathbb{E}[V_{t+1}(s',a')|s] \tag{1}\\ \text{s.t. }a'+c=(1+r_t)a+y_t(s) \tag{2}\\ a'\geq \underline{a} \tag{3} \\ \end{gather} $$

Optimal policy characterized by two conditions, which only require derivative of value function $V_a$ with respect to assets ("marginal value function"), not value function itself.

Envelope condition: $$ V_{a,t}(s,a)=(1+r_t)u'(c_t(s,a)) \tag{8} $$

First-order condition (inequality binds unless borrowing constraint binds):

$$ u'(c_t(s,a)) \geq \beta \mathbb{E}[V_{a,t+1}(s',a_t'(s,a))|s] \tag{9} $$

More efficient and accurate to use (8) and (9) to iterate directly on $V_a$ than to do standard value function iteration (Carroll 2006, "endogenous gridpoints method")

Backward iteration in time: take tomorrow's $V_{a,t+1}$ on right in (9), solve for today's policies $a_t'(s,a)$ and $c_t(s,a)$ and marginal value function $V_{a,t}$

Step 1: discounting and expectations

Define "end-of-period" value function $W_t$ as

$$ W_t(s,a') \equiv \beta \mathbb{E}[V_{t+1}(s',a')|s] \tag{10} $$

so that (9) simplifies to

$$ u'(c_t(s,a)) \geq W_{a,t}(s,a_t'(s,a)) \tag{11} $$

Take expectations in (10) by multiplying by Markov matrix $\Pi$, so if $V_{a,t+1}(s,a)$ is in $n_s\times n_a$ array Va, get $n_s\times n_a$ array Wa of $W_{a,t}$ by writing

Wa = (beta * Pi) @ Va

Step 2: solving for asset policy using first-order condition (hardest!)

  • Simplified first-order condition using $W_t$:
$$u'(c_t(s,a)) \geq W_{a,t}(s,a_t'(s,a)) \tag{11}$$
  • Assume for now this holds with equality, and define $a_t(s,a')$ to be inverse of $a_t'(s,a)$, i.e. $a_t(s,a_t'(s,a))\equiv a$, so (11) can be rewritten as
$$ c_t(s,a_t(s,a')) = \underbrace{(u')^{-1}(W_{a,t}(s,a'))}_{\equiv c_t^{endog}(s,a')} \tag{12} $$
  • This is time-$t$ consumption at the so-called endogenous gridpoint $(s,a_t(s,a'))$, which is at whatever asset level today that makes the household choose assets $a'$ tomorrow
  • This must obey budget constraint, where we define cash on hand to be $coh_t(s,a)\equiv (1+r_t)a + y_t(s)$:
$$ c_t^{endog}(s,a')+a' = coh_t(s,a(s,a')) \tag{13} $$
  • Calculating $c_t^{endog}$ and $coh_t$ in code, making the additional assumption that the utility function is isoelastic, i.e. $u(c) \equiv \frac{c^{1-\sigma^{-1}}}{1-\sigma^{-1}}$, with elasticity of intertemporal substitution $\sigma$ (eis in the code), so that $(u')^{-1}(x) = x^{-\sigma}$:
    c_endog = Wa**(-eis)
coh = y[:, np.newaxis] + (1+r)*a_grid

Step 2, continued

  • Take budget constraint (2)
$$ c_t^{endog}(s,a')+a' = coh_t(s,a(s,a')) = (1+r_t)a(s,a') + y_t(s) \tag{13} $$
  • Could isolate $a_t(s,a')$ by writing
$$ a_t(s,a') = \frac{c_t^{endog}(s,a') + a' - y_t(s)}{1+r_t} $$
  • Ultimately, though, we want the inverse (holding $s$ fixed) of this function, $a_t'(s,a)$, defined on our regular grid
  • Given $a_t(s,a')$, the easiest way to get that inverse is to evaluate $(a_t(s,a'),a')$ on the grid of all $a'$ for a given $s$, treating these as sample points from the function $a_t'(s,a)$, and then to linearly interpolate to evaluate $a_t'(s,a)$ on our grid of $a$

  • Turns out an equivalent (thanks to linearity) and slightly more convenient alternative is to think of (13) as giving a collection of sample points $(coh(s,a_t(s,a')), a')$ from the function mapping cash on hand to asset choice $a'$, then linearly interpolate to evaluate $a'$ for cash-on-hand $coh_t(s,a)$ on the actual grid of $a$, giving us $a_t'(s,a)$ on our grid of $a$

  • Need to do this separately for each $s$, so code iterates over them:

a = np.empty_like(coh)
for s in range(len(y)):
    a[s, :] = np.interp(coh[s, :], c_endog[s, :] + a_grid, a_grid)

Step 3: enforcing the borrowing constraint and backing out consumption

  • We just obtained policy function $a_t'(s,a)$ assuming Euler equation (11) held with equality, ignoring borrowing constraint $a'\geq \underline{a}$
  • Simple to fix: if $a_t'(s,a)$ we calculated is below $\underline{a}$, then (11) doesn't hold with equality after all, and instead constraint binds so $a_t'(s,a)=\underline{a}$
  • Hence, we should replace $a_t'(s,a)$ with $\max(a_t'(s,a),\underline{a})$
  • Then back out consumption from budget contraint (2), as whatever is left over from cash on hand after we choose assets: $c_t(s,a) = coh_t(s,a) - a_t'(s,a)$
  • Implementing in code (assuming first point on grid is $\underline{a}$):
a = np.maximum(a, a_grid[0])
c = coh - a

Step 4: use envelope condition to recover derivative of value function

Implement envelope condition (8), again using parametrization of $u$ where $u'(c)=c^{-1/\sigma}$, and $\sigma$ is eis in code:

Va = (1+r) * c**(-1/eis)

Combining everything into a single backward iteration

Combine the 9 lines of code we wrote for steps 1 to 4 into the body of a single backward iteration function

In [10]:
def backward_iteration(Va, Pi, a_grid, y, r, beta, eis):
    # step 1: discounting and expectations
    Wa = (beta * Pi) @ Va
    
    # step 2: solving for asset policy using the first-order condition
    c_endog = Wa**(-eis)
    coh = y[:, np.newaxis] + (1+r)*a_grid
    
    a = np.empty_like(coh)
    for s in range(len(y)):
        a[s, :] = np.interp(coh[s, :], c_endog[s, :] + a_grid, a_grid)
        
    # step 3: enforcing the borrowing constraint and backing out consumption
    a = np.maximum(a, a_grid[0])
    c = coh - a
    
    # step 4: using the envelope condition to recover the derivative of the value function
    Va = (1+r) * c**(-1/eis)
    
    return Va, a, c

This function takes in an $n_s\times n_a$ array Va containing tomorrow's marginal value function $V_{a,t+1}$ at every point on the grid. It also takes in the $n_s\times n_s$ Markov matrix Pi ($\Pi$), the length-$n_a$ asset grid a, the length-$n_s$ array y of incomes $y_t(s)$, and the scalars r ($r_t$), beta ($\beta$), and eis ($\sigma$).

It outputs $n_s\times n_a$ arrays containing today's marginal value function $V_{a,t}$ and the policies $a'_t$ and $c_t$ at every point on the grid.

Obtaining steady-state policies

To obtain steady-state $V_a$ and policies $a'$ and $c$, we need:

  • Initial guess for steady-state $V_a$
  • Iterate backward using backward_iteration until convergence criterion hit

For now we try:

  • $c$ equals 5% of cash-on-hand, use envelope condition to turn into guess for $V_a$
  • terminate backward iteration once the max distance between the asset policy function on two successive iterations falls below some tolerance, say $10^{-9}$
    • (also max number of iterations, say 10,000, as a fail-safe so that the computer doesn't run forever if it isn't converging)
In [11]:
def policy_ss(Pi, a_grid, y, r, beta, eis, tol=1E-9):
    # initial guess for Va: assume consumption 5% of cash-on-hand, then get Va from envelope condition
    coh = y[:, np.newaxis] + (1+r)*a_grid
    c = 0.05*coh
    Va = (1+r) * c**(-1/eis)
    
    # iterate until maximum distance between two iterations falls below tol, fail-safe max of 10,000 iterations
    for it in range(10_000):
        Va, a, c = backward_iteration(Va, Pi, a_grid, y, r, beta, eis)
        
        # after iteration 0, can compare new policy function to old one
        if it > 0 and np.max(np.abs(a - a_old)) < tol:
            return Va, a, c
        
        a_old = a

Calibrated example

Make grids for assets and income, with 500 points and 7 points respectively, and a borrowing constraint of 0. Use same calibration for income as we did above:

In [12]:
a_grid = discretize_assets(0, 10_000, 500)
y, pi, Pi = discretize_income(0.975, 0.7, 7)

Choose some values for the other parameters, with the idea that this is a quarterly calibration:

In [13]:
r = 0.01/4
beta = 1-0.08/4
eis = 1

Solve for steady-state policies:

In [14]:
Va, a, c = policy_ss(Pi, a_grid, y, r, beta, eis)

Visualizing consumption policy function

Looking at entire function, it seems linear, because at high assets it becomes asymptotically linear, and our grid includes really high asset levels to make sure we don't miss any rich households:

In [15]:
for s, ys in enumerate(y):
    plt.plot(a_grid, c[s, :], label=f'y={ys:.2f}')
plt.legend()
plt.xlabel('assets')
plt.ylabel('consumption');

Visualizing consumption policy function, continued

Looking only at lower assets, we see a more interesting picture, with some of the well-known concavity of the consumption function for lower income states:

In [16]:
for s, ys in enumerate(y):
    plt.plot(a_grid[:120], c[s, :120], label=f'y={ys:.2f}')
plt.legend()
plt.xlabel('assets')
plt.ylabel('consumption');

Visualizing consumption policy function, continued

It's even more visible if we zoom in on only the lowest income states:

In [17]:
for s, ys in enumerate(y[:4]):
    plt.plot(a_grid[:120], c[s, :120], label=f'y={ys:.2f}')
plt.legend()
plt.xlabel('assets')
plt.ylabel('consumption');

Visualizing consumption policy function, continued

It's also interesting to plot net saving—the increase or decrease in assets. For each income state, there is a "target" level of assets where net saving is zero. For lower incomes, this is zero (the borrowing constraint). For higher incomes, people are net savers until they reach a very high level of assets. You can see how this might produce a right-skewed distribution of assets!

In [18]:
for s, ys in enumerate(y):
    plt.plot(a_grid[:300], a[s, :300] - a_grid[:300], label=f'y={ys:.2f}')
plt.axhline(0, linestyle='--', color='k')
plt.legend()
plt.xlabel('assets')
plt.ylabel('savings');

Marginal propensities to consume

Another important concept is the marginal propensity to consume (MPC) out of a one-time, unexpected increase in cash on hand. This equals the derivative of consumption function with respect to assets, divided by asset income $1+r$.

Calculating it by taking differences of the consumption function to approximate the derivative:

In [19]:
mpcs = np.empty_like(c)
    
# symmetric differences away from boundaries
mpcs[:, 1:-1] = (c[:, 2:] - c[:, 0:-2]) / (a_grid[2:] - a_grid[:-2]) / (1+r)

# asymmetric first differences at boundaries
mpcs[:, 0]  = (c[:, 1] - c[:, 0]) / (a_grid[1] - a_grid[0]) / (1+r)
mpcs[:, -1] = (c[:, -1] - c[:, -2]) / (a_grid[-1] - a_grid[-2]) / (1+r)

# special case of constrained
mpcs[a == a_grid[0]] = 1

Visualizing marginal propensities to consume

Since the consumption function is concave, MPCs are declining in assets.

We see here that it seems to jump downward, e.g. from 1 to around 0.5 as households go from constrained to "not constrained, but probably constrained tomorrow, so we'll smooth any resources at the margin between 2 periods". These stark jumps exist in discrete time as long as we have discretized income states.

In [20]:
for s in range(7):
    plt.plot(a_grid[:50], mpcs[s, :50], linewidth=3, label=f'y = {y[s]:.2f}')
plt.xlabel('Assets')
plt.ylabel('MPC')
plt.title('Quarterly marginal propensities to consume by income state y(s)')
plt.legend();

Part 3 of computation: forward iteration to obtain the distribution

  • We'd like to obtain a discretized distribution giving the mass of agents at each point on our grid
  • Problem: in general, if we start at a gridpoint, the asset policy function takes us off the grid
  • For instance, look at asset policy for second-highest income state ($y = 2.46$) and assets of zero:
In [21]:
a[5, 0]
Out[21]:
0.4364460365195778

This turns out to be in between points 66 and 67 on the grid:

In [22]:
a_grid[60:70]
Out[22]:
array([0.38029038, 0.38883467, 0.3974721 , 0.40620394, 0.41503147,
       0.423956  , 0.43297886, 0.44210141, 0.45132502, 0.46065108])

Where should we assume this mass of households goes?

Solution: lotteries

  • Assume household follows a "lottery", going to each of the nearest two gridpoints with some probability, such that on average, the household chooses the right asset value
  • We can think of this as a local first-order approximation
  • In particular, let the asset grid be $(a_0,\ldots,a_{n_a-1})$. If some policy $a'$ lies in between gridpoints $a_i$ and $a_{i+1}$, then we want a lottery where the household goes to $a_i$ with probability $\pi$ and $a_{i+1}$ with probability $1-\pi$, such that
$$ \pi a_i + (1-\pi) a_{i+1} = a' $$
  • This can easily be solved to give
$$ \pi = (a_{i+1} - a')/(a_{i+1} - a_i) \tag{14} $$
  • Analogous to linear interpolation: replacing $a'$ with mix of $\pi$ at $a_i$ and $1-\pi$ at $a_{i+1}$ means replacing $f(a')$ with $\pi f(a_i)+(1-\pi)f(a_{i+1})$, for any function $f$

Obtaining the lottery

For any individual policy $a'$, we can (1) identify bracketing gridpoints $a_i$ and $a_{i+1}$ and (2) calculate $\pi$ from (14). Following code takes in $a'$ as a and returns $i$ and $\pi$ as a_i and a_pi.

In [23]:
def get_lottery(a, a_grid):
    # step 1: find the i such that a' lies between gridpoints a_i and a_(i+1)
    a_i = np.searchsorted(a_grid, a) - 1
    
    # step 2: obtain lottery probabilities pi
    a_pi = (a_grid[a_i+1] - a)/(a_grid[a_i+1] - a_grid[a_i])
    
    return a_i, a_pi

Try it for single point a[5,0] we saw earlier, which was between gridpoints 66 and 67

In [24]:
get_lottery(a[5,0], a_grid)
Out[24]:
(66, 0.6199338577164767)

This gives us $i=66$ and $\pi\approx 0.62$, with latter indicating that a[5,0] is closer to gridpoint 66 than 67, so it gets higher weight in lottery.

This function also works (NumPy "broadcasts" it) if we apply it to entire policy function a, with outputs having same $n_s\times n_a$ dimensionality as a:

In [25]:
a_i, a_pi = get_lottery(a, a_grid)
a_i.shape, a_i.dtype, a_pi.shape, a_pi.dtype
Out[25]:
((7, 500), dtype('int64'), (7, 500), dtype('float64'))

Using lottery to iterate on distribution

  • Suppose we have discretized distribution $D_t(s,a)$ with mass at each gridpoint
  • Send $\pi(s,a)$ of mass to gridpoint at $i(s,a)$, $1-\pi(s,a)$ to gridpoint at $i(s,a)+1$
  • Easy to implement directly, but too slow using for loops in Python, so use @numba.njit decorator to accelerate with Numba's just-in-time compilation
In [26]:
@numba.njit
def forward_policy(D, a_i, a_pi):
    Dend = np.zeros_like(D)
    for s in range(a_i.shape[0]):
        for a in range(a_i.shape[1]):
            # send pi(s,a) of the mass to gridpoint i(s,a)
            Dend[s, a_i[s,a]] += a_pi[s,a]*D[s,a]
            
            # send 1-pi(s,a) of the mass to gridpoint i(s,a)+1
            Dend[s, a_i[s,a]+1] += (1-a_pi[s,a])*D[s,a]
            
    return Dend
  • This updates beginning-of-period distribution $D_(s,a)$ to end-of-period distribution $D^{end}_t(s,a')$, after asset choice $a'$ made but before new draw of $s'$
  • Note: could also try writing lottery as a sparse Markov matrix using sparse matrix libraries, but this turns out to be both slower and more cumbersome

Iterating on the state $s$ as well

  • Previous function takes us from $D_(s,a)$ to "end-of-period" $D_t^{end}(s,a')$
  • We want $D_{t+1}(s',a')$, i.e. to update $s$ to new draw $s'$
  • Like in stationary_markov earlier, multiply by transpose $\Pi'$ of Markov matrix $\Pi(s,s')$
    • note contrast with first step of backward_iteration, where we multiply by $\Pi$ to take expectations
  • Combine with forward_policy function to get complete forward iteration from $D_t(s,a)$ to $D_{t+1}(s',a')$
In [27]:
def forward_iteration(D, Pi, a_i, a_pi):
    Dend = forward_policy(D, a_i, a_pi)    
    return Pi.T @ Dend

Iterating to find steady-state distribution of $(s,a)$

  • As initial guess for steady-state $D(s,a)$, use $D(s,a)=\frac{1}{n_a}\pi(s)$
    • Right distribution for $s$, uniform on $a$
  • Iterate until tolerance is reached
In [28]:
def distribution_ss(Pi, a, a_grid, tol=1E-10):
    a_i, a_pi = get_lottery(a, a_grid)
    
    # as initial D, use stationary distribution for s, plus uniform over a
    pi = stationary_markov(Pi)
    D = pi[:, np.newaxis] * np.ones_like(a_grid) / len(a_grid)
    
    # now iterate until convergence to acceptable threshold
    for _ in range(10_000):
        D_new = forward_iteration(D, Pi, a_i, a_pi)
        if np.max(np.abs(D_new - D)) < tol:
            return D_new
        D = D_new

Applying to our calibrated example

We get an $n_s\times n_a$ array giving mass at each $D(s,a)$

In [29]:
D = distribution_ss(Pi, a, a_grid)
D.shape
Out[29]:
(7, 500)

How to visualize? Can try plotting CDF of assets, obtained by summing D across $s$ and taking cumulative sum across $a$.

Turns out not to be so informative, since virtually everyone has assets far below the top of the distribution (we made grid maximum higher than necessary):

In [30]:
plt.plot(a_grid, D.sum(axis=0).cumsum());

Visualizing steady-state distribution

  • Better to plot up to a much lower level of assets, e.g. up to $a=20$
    • (20x average quarterly income; 5x average annual income)
  • Only small fraction has assets above $a=20$ (but they will turn out to hold a lot of assets!); half hold nothing at all
In [31]:
i = np.argmax(a_grid > 20) # first gridpoint above 20
plt.plot(a_grid[:i], D.sum(axis=0)[:i].cumsum())
plt.hlines(1,0, a_grid[i-1], colors='k', linestyles='dotted');

Visualizing steady-state distribution, continued

  • Zoom in further, up to $a=2$
  • Start to see some interesting features
    • Jumps in CDF near zero that look like mass points of distribution
    • Kink just above 1 that looks like a sudden drop in density
In [32]:
i = np.argmax(a_grid > 2) # first gridpoint above 2
plt.plot(a_grid[:150], D.sum(axis=0)[:150].cumsum())
plt.hlines(1,0, a_grid[149], colors='k', linestyles='dotted');

Visualizing steady-state distribution

  • Plot CDF separately by income state, highly right-skewed in upper states, interesting behavior mostly in $y=1.39$ middle state
  • Why? Common occurences with discrete time and discrete states:
    • Mass point at $a=0$ for $y < 1.39$ leads to mass points above 0, as households whose income has increased to $y=1.39$ save for one, two, three, ... periods
    • Target assets for $y=1.39$ is slightly above 1; $y=1.39$ households converge to this level from both sides but don't cross, more converge from below so density discontinuous
In [33]:
for s, ys in enumerate(y):
    plt.plot(a_grid[:150], D[s][:150].cumsum()/pi[s], label=f'y={ys:.2f}')
plt.legend();

How many assets held in each state?

Let's get a sense of extent to which higher incomes hold far higher assets, even relative to income:

In [34]:
for s, ys in enumerate(y):
    print(f'Ave assets at y={ys:.2f}: {np.vdot(a_grid, D[s, :]) / pi[s]:.2f}')

Ave assets at y=0.14: 0.02
Ave assets at y=0.25: 0.05
Ave assets at y=0.44: 0.16
Ave assets at y=0.79: 0.56
Ave assets at y=1.39: 2.19
Ave assets at y=2.46: 7.01
Ave assets at y=4.36: 17.67

What fraction of assets held in each state?

Higher income states hold a lot more, but rarer. Still, top few hold most assets:

In [35]:
for s, ys in enumerate(y):
    print(f'Total assets at y={ys:.2f} ({100*pi[s]:.1f}% of hh)\t: {np.vdot(a_grid, D[s, :]):.2f}')

Total assets at y=0.14 (1.6% of hh)	: 0.00
Total assets at y=0.25 (9.4% of hh)	: 0.00
Total assets at y=0.44 (23.4% of hh)	: 0.04
Total assets at y=0.79 (31.2% of hh)	: 0.17
Total assets at y=1.39 (23.4% of hh)	: 0.51
Total assets at y=2.46 (9.4% of hh)	: 0.66
Total assets at y=4.36 (1.6% of hh)	: 0.28

Sum of these is total assets held in economy:

In [36]:
np.sum(a_grid*D) # use array broadcasting since a_grid missing 's' dimension
Out[36]:
1.6645070676480889

which should equal (up to numerical error) total assets chosen for tomorrow:

In [37]:
np.vdot(a, D)
Out[37]:
1.6645070661939019

Part 4 of computation: aggregating to full steady state

  • Simplest part: combine parts 2 and 3 and aggregate to get total consumption and assets in one big function
  • Return a dict with all inputs to the problem (Pi, a_grid, y, r, beta, eis), the individual level policies and marginal value function (a, c, Va), the distribution (D), and aggregate policies (A and C)
In [38]:
def steady_state(Pi, a_grid, y, r, beta, eis):
    Va, a, c = policy_ss(Pi, a_grid, y, r, beta, eis)
    D = distribution_ss(Pi, a, a_grid)
    
    return dict(D=D, Va=Va, 
                a=a, c=c,
                A=np.vdot(a, D), C=np.vdot(c, D), # aggregation
                Pi=Pi, a_grid=a_grid, y=y, r=r, beta=beta, eis=eis)
In [39]:
ss = steady_state(Pi, a_grid, y, r, beta, eis)
ss.keys()
Out[39]:
dict_keys(['D', 'Va', 'a', 'c', 'A', 'C', 'Pi', 'a_grid', 'y', 'r', 'beta', 'eis'])

Application: comparative statics in $r$

Steady-state assets slope sharply upward as we increase the steady-state real interest rate $r$, and indeed asymptote to infinity as we approach $r=\beta^{-1}-1$, which in this case is just above 0.02:

In [40]:
rs = r + np.linspace(-0.02, 0.015, 15)
As = [steady_state(Pi, a_grid, y, r, beta, eis)['A'] for r in rs]
plt.plot(rs, As)
plt.xlabel('Real interest rate')
plt.ylabel('Aggregate assets over quarterly income');

Application: comparative statics in income risk

Steady-state assets also slope strongly upward in income risk, which we'll manipulate below by calibrating to different standard deviations $\sigma$ of cross-sectional log income. (This turns out to change only the vector y of incomes, making them more dispersed, not the Markov matrix Pi.)

In [41]:
sigmas = np.linspace(0.3, 1.2, 8) # our benchmark was sigma=0.7
As = []
for sigma in sigmas:
    y_new, pi_new, Pi_new = discretize_income(0.975, sigma, 7)
    As.append(steady_state(Pi_new, a_grid, y_new, r, beta, eis)['A'])
plt.plot(sigmas, As)
plt.xlabel('Cross-sectional standard deviation of income')
plt.ylabel('Assets over quarterly income');

Application: comparative statistics in intertemporal substitution

  • Steady-state assets are also decreasing in the elasticity of intertemporal substitution, or equivalently increasing in the coefficient of risk aversion (its inverse)
  • This is because with lower risk aversion and less desire to smooth consumption, the precautionary and income-smoothing motives for saving become weaker
In [42]:
eis_vec = np.linspace(0.4, 2, 10) # our benchmark was eis=1
As = [steady_state(Pi, a_grid, y, r, beta, eis)['A'] for eis in eis_vec]
plt.plot(eis_vec, As)
plt.xlabel('Elasticity of intertemporal substitution')
plt.ylabel('Assets over quarterly income');

What remains?

  • We've learned how to compute the full steady state—including the steady state policy functions, the steady state distribution, and aggregates—of the standard incomplete markets model
  • Remaining topics to briefly cover
    • Calibration: choosing some parameters to hit steady-state calibration targets, rather than exogenously specifying all parameters
    • General equilibrium: how to think about the SIM as part of a very simple general equilibrium model
    • Expectation vectors: a tool to compute almost any desired moment of individual behavior, which will also come in handy later as part of the fake news algorithm

Calibration

  • Suppose we want to choose total assets to equal a certain level relative to income
    • say, 140% of annual household income (inspired by McKay, Nakamura, and Steinsson, who calibrate to 140% of annual GDP)
    • average quarterly income normalized to 1, so with this target A must be $4\times 1.4=5.6$
    • calibrate all other parameters directly, choose discount factor $\beta$ to hit this target
    • for one-dimensional root-finding use the robust and quick Brent's method
In [43]:
beta_calib = optimize.brentq(
                lambda beta: steady_state(Pi, a_grid, y, r, beta, eis)['A'] - 5.6,
                0.98, 0.995)
beta_calib
Out[43]:
0.987703940322874

Get entire steady state and verify we have right assets:

In [44]:
ss_calib = steady_state(Pi, a_grid, y, r, beta_calib, eis)
ss_calib['A'], ss_calib['C']
Out[44]:
(5.5999999999694, 1.01400000488598)

Check aggregate steady-state budget balance $C=1+rA$ (up to numerical error from our steady-state tolerance):

In [45]:
ss_calib['C'] - (1 + ss_calib['r']*ss_calib['A'])
Out[45]:
4.886056670372341e-09

General equilibrium

  • So far we just discussed the household problem in isolation ("partial equilibrium")
  • How to think about GE?
    • Suppose households have stochastic endowments $e(s)$ of effective labor (average of 1), following same Markov process that $y(s)$ previously did
    • Each unit produces one unit of consumption goods in a perfectly competitive production sector, so real wage is 1
    • Households save in government bonds, paying real interest rate $r$ and supplied in exogenous quantity $B$
    • Government finances this with tax $\tau = rB$ on labor, households' after-tax labor income is $y(s)=(1-\tau)e(s)$
  • Only difference here from before is that $y(s)$ in households' budget constraint includes factor of $1-\tau$

General equilibrium calibration

Calibrate similar to before:

In [46]:
B = 5.6   # annual bonds/GDP is 140%, so quarterly is 560%, and quarterly GDP is 1
tau = r*B # labor tax needed to balance steady-state government budget
e = y     # use our previous y, which had mean 1, for the labor endowment process

Because of the low real interest rate, $\tau$ is quite small:

In [47]:
tau
Out[47]:
0.013999999999999999

Now calibrate $\beta$ to be consistent with steady-state asset market clearing $A=B$, using $(1-\tau)e(s)$ as income:

In [48]:
beta_ge = optimize.brentq(
                lambda beta: steady_state(Pi, a_grid, (1-tau)*e, r, beta, eis)['A'] - B,
                0.98, 0.995)
beta_ge
Out[48]:
0.9877855433151486

Compute full steady state with beta_ge, check asset market and goods market clearing:

In [49]:
ss_ge = steady_state(Pi, a_grid, (1-tau)*e, r, beta_ge, eis)
ss_ge['A'] - B, ss_ge['C'] - 1 # asset mkt clearing, goods mkt clearing
Out[49]:
(-1.651745407116323e-11, 4.785916551952596e-09)

General equilibrium counterfactuals

  • Now that we have full calibrated GE model, can do GE comparative statics
    • Vary some parameter and re-solve for GE real interest rate $r$, keeping all other parameters constant
  • Example: vary income risk like before
    • Now total assets are fixed and can't change, but $r$ falls to equilibriate the market as risk rises and households want to save more
    • (Note this computation is a little slow, since we are looping over $\sigma$ and solving separately for each $r$. Smarter ways exist to do this, which we won't cover now.)
In [50]:
sigmas = np.linspace(0.3, 1.2, 8) # our benchmark was sigma=0.7
rs = []
for sigma in sigmas:
    e_new, pi_new, Pi_new = discretize_income(0.975, sigma, 7)
    rs.append(optimize.brentq(
                lambda r: steady_state(Pi_new, a_grid, (1-r*B)*e_new, r, beta_ge, eis)['A'] - B,
                -0.02, 0.015))
plt.plot(sigmas, 4*np.array(rs))
plt.xlabel('Cross-sectional standard deviation of income')
plt.ylabel('Equilibrium real interest rate (annualized)');

Expectation vectors

  • Suppose we have some function $x(s,a)$ of the state
  • If the state for an individual household is $(s_t,a_t)$ today, what do we expect $x(s_{t+j},a_{t+j})$ to be, assuming household follows steady-state policy rule?
  • Call this the $j$th expectation vector $\mathcal{E}_j$:
$$ \mathcal{E}_j(s,a) \equiv \mathbb{E}[x(s_{t+j}, a_{t+j}) | (s_t, a_t)] $$
  • Insight: law of iterated expectations implies a simple recursion with base case $\mathcal{E}_0(s,a) = x(s,a)$:
$$ \mathcal{E}_j(s,a) = \mathbb{E}[\mathcal{E}_{j-1}(s',a') | (s,a)] \tag{15} $$
  • So as long as we know how to calculate one-period-ahead expectations to evaluate right side of (15), we can calculate expectation vectors at any horizon!
    • This will enable us to do a lot!
    • We'll turn to one-period-ahead expectations now...

Part one of one-period-ahead expectations: take expectations with respect to lottery

Recall that given a_i and a_pi describing the lottery form of policy, we updated the beginning-of-period distribution $D_t(s,a)$ to an end-of-period distribution $D_t^{end}(s, a')$ using:

In [51]:
@numba.njit
def forward_policy(D, a_i, a_pi):
    Dend = np.zeros_like(D)
    for s in range(a_i.shape[0]):
        for a in range(a_i.shape[1]):
            # send pi(s,a) of the mass to gridpoint i(s,a)
            Dend[s, a_i[s,a]] += a_pi[s,a]*D[s,a]
            
            # send 1-pi(s,a) of the mass to gridpoint i(s,a)+1
            Dend[s, a_i[s,a]+1] += (1-a_pi[s,a])*D[s,a]
            
    return Dend

Here, we need to go in the opposite direction: start with some end-of-period variable $X_t^{end}(s,a')$, and obtain its beginning-of-period expectation $X_t(s,a) = \pi(s,a)*X_t^{end}(s,i(s,a)) + (1-\pi(s,a))*X_t^{end}(s,i(s,a)+1)$. Though we don't write it explicitly as a matrix, one can verify that this is actually the transpose of the linear mapping in forward_policy:

In [52]:
@numba.njit
def expectation_policy(Xend, a_i, a_pi):
    X = np.zeros_like(Xend)
    for s in range(a_i.shape[0]):
        for a in range(a_i.shape[1]):
            # expectation is pi(s,a)*Xend(s,i(s,a)) + (1-pi(s,a))*Xend(s,i(s,a)+1)
            X[s, a] = a_pi[s, a]*Xend[s, a_i[s, a]] + (1-a_pi[s, a])*Xend[s, a_i[s, a]+1]
            
    return X

Combine with taking expectations over the state

Recall that for the distribution, we updated $D_t^{end}(s,a')$ to $D_{t+1}^{end}(s',a')$ by multiplying by the transpose of the Markov matrix, then combined this with the update using the lottery:

In [53]:
def forward_iteration(D, Pi, a_i, a_pi):
    Dend = forward_policy(D, a_i, a_pi)    
    return Pi.T @ Dend

Similarly, here we need to obtain $X_t^{end}(s,a') = \mathbb{E}_t[X_{t+1}(s',a') | s]$ by taking expectations with respect to $s'$, then combine this with the code we already wrote to obtain $X_t(s,a)$ from $X_t^{end}(s,a')$. Altogether, this gives us a function that iterates

In [54]:
def expectation_iteration(X, a_i, a_pi):
    Xend = Pi @ X
    return expectation_policy(Xend, a_i, a_pi)

Note that taking expectations goes in the opposite order (first take expectations with respect to $s'$, then $a'$) as updating the distribution, because it's going backward in time rather than forward.

Obtaining expectation vectors

  • Write a simple function that takes in any function $x(s,a)$ of the state (in the form of an $n_s\times n_a$ array X) and calculates expectation vectors $\mathcal{E}_0,\ldots,\mathcal{E}_{T-1}$ up to some horizon $T$.
  • Use recursion based on law of iterated expectations:
$$ \mathcal{E}_j(s,a) = \mathbb{E}[\mathcal{E}_{j-1}(s',a') | (s,a)] \tag{15} $$
  • Combined with base case $\mathcal{E}_0(s,a) \equiv x(s,a)$, just repeatedly apply expectation_iteration to evaluate right side of (15):
In [55]:
def expectation_vectors(X, a_i, a_pi, T):
    # set up array of curlyEs and fill in first row with base case
    curlyE = np.empty((T, ) + X.shape)
    curlyE[0] = X
    
    # recursively apply law of iterated expectations
    for j in range(1, T):
        curlyE[j] = expectation_iteration(curlyE[j-1], a_i, a_pi)
        
    return curlyE

Application: autocorrelation of consumption

  • Want to calculate correlation between individual consumption at dates $t$ and $t+j$ in steady state
  • To do this first get covariances (where $\tilde{c}(s,a)\equiv c(s,a) - C$ is demeaned consumption and $\mathcal{E}_j^{\tilde{c}}$ are expectation vectors of $\tilde{c}$):
$$ \begin{align} \text{Cov}(c(s_t,a_t), c(s_{t+j},a_{t+j})) &= \mathbb{E}[\tilde{c}(s_t,a_t)\cdot \tilde{c}(s_{t+j},a_{t+j})] \\ &= \mathbb{E}[\tilde{c}(s_t,a_t) \cdot \mathbb{E}[\tilde{c}(s_{t+j},a_{t+j}) | (s_t,a_t)]] \\ &= \mathbb{E}[\tilde{c}(s_t,a_t) \cdot \mathcal{E}_j^{\tilde{c}}(s_t,a_t)] \end{align} $$
  • In other words, we just need to evaluate the product of $\tilde{c}$ and its $j$-period-ahead expectation $\mathcal{E}_j^{\tilde{c}}$ at all $(s,a)$, then aggregate across distribution
    • Implement on next slide!

Application: autocorrelation of consumption, continued

Let's try this for all $j$ up to $T=40$:

In [66]:
T = 40
ctilde = c - np.vdot(D, c) # demeaned consumption
E_ctilde = expectation_vectors(ctilde, a_i, a_pi, T) # expectation vectors for 0, ... , T-1

Need to take product of ctilde with the $j$th row of E_ctilde, then aggregate by steady-state distribution D.

In [72]:
Autocov_c = np.empty(T)
for j in range(T):
    Autocov_c[j] = np.vdot(D, ctilde * E_ctilde[j])

Autocorrelation is autocovariance divided by variance, and variance is autocovariance at a lag of zero, so can write:

In [73]:
Autocorr_c = Autocov_c / Autocov_c[0]

Python trickery side note: not necessary here for speed, but can calculate Autocov_c above without a for loop by flattening state space $(s,a)$ into one dimension using ravel (for ctilde*D) and reshape (for E_ctilde) and then doing matrix multiplication:

In [77]:
Autocov_c_alt = E_ctilde.reshape((T, -1)) @ (ctilde*D).ravel()
assert np.max(np.abs(Autocov_c_alt - Autocov_c)) < 1E-15 # verify this is right

Application: autocorrelation of consumption, visualized

  • Applaud ourselves for calculating this exactly, and almost instantly, without any Monte Carlo simulation!
  • Consumption pretty persistent here, reflecting both persistence in underlying income process and smoothing via assets
  • But maybe not persistent enough? Empirically, population-level correlation between $c_t$ and $c_{t+j}$ probably above .50 even at 10-year horizon (unmodeled permanent component of income?)
In [80]:
plt.plot(Autocorr_c)
plt.title('Autocorrelation of consumption')
plt.xlabel('Horizon in quarters')
plt.ylabel('Correlation');