Chapter 10.1: Overlapping Generations Models with Idiosyncratic Income Uncertainty --- Python Code

book_cover

By BURKHARD HEER and ALFRED MAUSSNER, University of Augsburg

In the following, we study the computation of the steady state of the 70-Period OLG model described in Section 10.1 of our book on 'Dynamic General Equilibrium Modeling: Computational Methods and Applications' (3rd edition, in progress). We consider a medium-scale overlapping generations model where workers are subject to idiosyncratic productivity shocks so that income is stochastic at the individual level (but not at the aggregate level).

In the following, we present a detailed discussion of the PYTHON program that computes the stationary solution to this heterogeneous-agent economy. We apply value function iteration to solve for the individual policy functions. In order to maximize the right-hand side of the Bellman equation, we use Golden Section Search as described in Section 11.6.1 of our book. In order to interpolate between grid points of the value function, we use linear or cubic interpolation as described in Section 11.2.3 of our book. In addition, we propose an algorithm for the computation of the stationary distribution.

The PYTHON program 'AK70_stochastic_income.py' can be downloaded from our web page

https://www.uni-augsburg.de/de/fakultaet/wiwi/prof/vwl/maussner/dgebook/

If you have any questions, suggestions or comments, please send an email to: Burkhard.Heer@wiwi.uni-augsburg.de

Description of the model

We study a medium-sized OLG model with respect to its ability to model inequality. It is characterized by the following features:

  1. life-cycle savings,
  2. uncertain lifetime,
  3. uncertain earnings,
  4. pay-as-you-go pensions.

Huggett (1996) analyses this kind of model with respect to the income and wealth distribution; in particular, he wonders if the OLG model is able to explain the fact that the (endogenous) wealth inequality is larger than the earnings and income inequality. We modify his model slightly.

Demographics

A period, $t$, corresponds to one year. At each period $t$, a new generation of households is born. Newborns have a real-life age of 21 denoted by $s=1$. All generations retire at the end of age $s=T^W=45$ (corresponding to a real-life age of 65) and live up to a maximum age of $s=T=70$ (real-life age 90). The number of periods during retirement is equal to $T^R=T-T^W=25$.

Let $N_t(s)$ denote the number of agents of age $s$ at $t$. We denote the total population at $t$ by $N_t$. At $t$, all agents of age $s$ survive until age $s+1$ with probability $\phi_{t}^s$, where $\phi_{t}^0=1$ and $\phi_{t}^T=0$. In period $t$, the newborn cohort grows at rate $n_t$: \begin{equation} N_{t}(1) = (1+n_t) N_{t-1}(1).\end{equation} We consider the steady state where the survival probabilities and the population growth rate are constant, $\phi^s_t=\phi^s$ and $n_t=n$.

Households

Each household comprises one (possibly retired) worker. Households maximize expected intertemporal utility at the beginning of age 1 in period $t$ \begin{equation} \label{Ch93lifeutil2} \max \sum_{s=1}^J \beta^{s-1} \left(\Pi_{j=1}^s \phi_{t+j-1}^{j-1}\right)\; \mathbb{E}_t \; \left[ u(c_{t+s-1}^{s}, l_{t+s-1}^{s} ) + v(g_{t+s-1})\right],\end{equation} where $ \beta>0$ denotes the discount factor. Instantaneous utility $u(c,1-l)$ is specified as a function of consumption $c$ and and leisure $1-l$: \begin{equation} \label{util_instantaneous} u(c,1-l)= \frac{\left(c^\gamma (1-l)^{1-\gamma} \right)^{1-\eta}}{1-\eta} \end{equation} where $1/\eta$ denotes the intertemporal elasticity of substitution (IES) and $\gamma$ is the share of consumption in utility.

During working life, the labor supply of the $s$-year-old household amounts to $0 \le l^s\le l^{max}$, $s=1,\ldots,T^W$, where we impose a constraint on the maximum working hours $l^{max}$. During retirement, $l^s_t\equiv 0$ for $s=T^W+1,\ldots,T$.

Utility from government consumption, $v(g_t)$, is additive, so government consumption per capita, $g_t$, does not have any direct effect on household behavior (only indirectly through its effects on transfers and taxes).

Total gross labor income of the $s$-year-old worker in period $t$, $\epsilon(s,\theta,e) A_t w_t l^s_t$, consists of the product of his idiosyncratic productivity $\epsilon(s,\theta,e)$, aggregate productivity $A_t$, the wage per efficiency unit $w_t$ and his working time $l^s_t$. Aggregate productivity $A_t$ grows at the exogenous rate $g_A$, which is also equal to the economic growth rate of per capita income in steady state. The household's idiosyncratic labor productivity $\epsilon(s,\theta,e)$ is stochastic and also depends on his age $s$ according to $\epsilon(s,\theta,e)=\theta e \bar y^s$. The stochastic component $\theta$ follows a Markov process and takes only a finite number $n_\theta$ of possible values in the set $\Theta=\{\theta_1=\underline{\theta},\ldots,\theta_{n_\theta}=\overline{\theta}\}$ with $\theta_i<\theta_{i+1}$ for $i=1,\ldots,n_\theta-1$. Let $prob(\theta'|\theta)$ denote the transition probability of the household with productivity $\theta$ in this period to become a household with productivity $\theta'$ in the next period for households aged $s=1,\ldots,T^W-1$. In old age, we assume without loss of generality that productivity remains constant. Notice furthermore that we assume that the transition matrix is independent of age $s< T^W$ and time-invariant. The individual $\theta$'s are assumed to be independent across agents and the law of large numbers holds (there is an infinite number of agents) so that there is no aggregate uncertainty. The labor productivity process is calibrated in detail below. In addition, the individual's productivity depends on his permanent productivity type $e\in\{e_1,e_2\}$ which is chosen to reflect the difference in earnings that stem from different levels of education (high school/college) and an age compenent $\bar y^s$, which is hump-shaped over the life-cycle.

All households receive transfers $tr_t$ from the government. The worker pays labor income taxes $\tau^l$ and social security contribution $\tau^p$ proportional to his labor income,
$\epsilon(s,\theta,e) A_t w_t l^s_t$. In old age, the retired worker receives a pension $pen_t$ which is independent of the individual's contributions. (In Chapter 10.1.2., we consider the case of contributions-dependent pensions.)

Accordingly, net non-capital income $y^s_t$ is presented by \begin{equation} y^s_t = \left\{ \begin{array}{ll} (1-\tau^l_t-\tau^p_t) \epsilon(s,\theta,e) A_t w_t l^{s}_{t} \;\; & s=1,\ldots, T^W,\\ &\\ pen_t & s= T^W+1,\ldots,T. \end{array}\right. \end{equation}

To express the budget constraint in stationary variables, we divide non-stationary individual variables by aggregate productivity $A_t$, e.g. $\tilde y^s_t \equiv y^s_t/A_t$. The budget constraint of the household at age $s=1,\ldots,T$ is then given by \begin{equation} (1+\tau^c_t) \tilde c^{s}_{t} = \tilde y^s_t+ (1-\tau^k_t) (r_{t} -\delta) \tilde k^{s}_{t} +\tilde k^{s}_{t} + R_{t}^b \tilde b^{s}_{t} + \tilde tr_{t} - (1+g_A) \tilde k_{t+1}^{s+1}-(1+g_A) \tilde b_{t+1}^{s+1},\end{equation} where $k^{s}_t$ and $b^{s}_t$ denote the capital stock and government bonds of the $s$-year-old agent at the beginning of period $t$. The household is born without assets and leaves no bequests at the end of its life, implying $k^{1}_t=k^{T+1}_t=0$ and $b^{1}_t=b^{T+1}_t=0$. She receives interest income $r_t$ and $r^b_t = R^b_t-1$ on capital and government bonds and pays income taxes on labor and capital income at the rates of $\tau^l_t$ and $\tau^k_t$, respectively. Capital depreciation $\delta k^{s}_t$ is tax exempt. Consumption is taxed at the rate $\tau^c_t$.

In addition, we impose a non-negativity constraint on assets, $b^s_t\ge 0$ and $k^s_t\ge 0$.

Technology

Output is produced with the help of capital $K_t$ and effective labor $L_t$ according to the standard Cobb-Douglas function: \begin{equation} Y_{t} = K_{t}^\alpha \left(A_t L_{t}\right)^{1-\alpha}.\end{equation}

Firms are competitive and maximize profits $\Pi_t=Y_t-r_t K_{t}-w_t A_t L_{t}$ such that factor prices are given by:

\begin{align} \label{firm1} w_t&=(1-\alpha) K_{t}^\alpha \left( A_t L_{t}\right)^{-\alpha},\\ \label{firm2} r_t&=\alpha K_{t}^{\alpha-1} \left( A_t L_{t}\right)^{1-\alpha}. \end{align}

Government and Social Security

The government levies income taxes $\tau^l$ and $\tau^k$ on labor and capital income and taxes $\tau^c$ on consumption. In addition, the government confiscates all accidental bequests $Beq_t$. It pays aggregate transfers $Tr_t$, provides a certain level $G_t$ of total public expenditures, and pays interest $r^b_t$ on the accumulated debt $B_t$. In each period, the government budget is financed by issuing government debt: \begin{equation} \label{Ch10govbudget} Tr_t+ G_t +r^b_t B_t - Tax_t - Beq_t = B_{t+1}-B_t,\end{equation} where taxes $Tax_t$ are given by \begin{equation} Tax_t = \tau^l_t A_t L_t w_t + \tau^k_t (r_t-\delta) K_t + \tau^c_t C_t, \end{equation} and $C_t$ denotes aggregate consumption.

The government provides pay-as-you-go pensions to the retirees which it finances with the contributions of the workers. Let $Pen_t$ denote aggregate pension payments. The social security budget is assumed to balance: \begin{equation} Pen_t = \tau^p_t A_t L_t w_t. \end{equation}

Households' Optimization Problem

In the following, we formulate the household optimization problem in recursive form. Therefore, we first make an assumption with respect to the portfolio choice of the individuals. In equilibrium, the allocation on the two assets is indeterminate as we have many agents and the after-tax returns on both assets, $K$ and $B$, are equal. Therefore, we make the innocuous assumption that all households hold both assets in the same proportion, which is equal to $K/(K+B)$ and $B/(K+B)$, respectively. In addition, we define households assets $a^s_t$ as the sum of the two individual assets, $k^s_t+b^s_t$.

Let the state vectors be given by $z=(s,\theta,e,a)$. In the following, we will use the abbreviation $z_s=(\theta,e,a)$ for the individual state vector of the $s$-year-old agent. For notational convenience, we drop the time-period index $t$ in the following whenever appropriate. For example, $c$, $l$, $y$ and $a$ will denote individual consumption, labor, net income and wealth at age $s$ in period $t$.

$V_t(z_s)$ denotes the value function of the $s$-year old household in period $t$. The optimization problem of the household is given by

\begin{equation} V_t(z_s) = \max_{c,l,a'} \left\{ u(c,l)+v(g) + \beta \phi^s_t \sum_{\theta'} \; prob(\theta'|\theta)\; V_{t+1}(z_{s+1})\right\} ,\end{equation}

subject to:

\begin{align} (1+\tau^c_t) c &= y +\left[ 1+ (1-\tau^k_t) (r_{t} -\delta) \right] a +tr - a',\\ a&\ge 0, \end{align}

with the terminal condition $V_t\left(z_{T+1}\right)=0$.

To formulate this optimization problem in stationary variables, define

$\tilde V_t \equiv \frac{V_t}{A_t^{\gamma(1-\eta)}}$, $\tilde z_s = (\theta,e,\tilde a) $, $\tilde v(\tilde g) = \frac{v(g)}{A_t^{\gamma(1-\eta)}}$, and

\begin{equation} \tilde u(\tilde c,1-l) = \frac{u(c,l)}{A_t^{\gamma(1-\eta)}} = \frac{\left(\tilde c^{\gamma} (1-l)^{1-\gamma}\right)^{1-\eta}}{1-\eta} \end{equation}

The Bellman equation can be rewritten in stationary form \begin{equation} \tilde V_t(\tilde z_s) = \max_{\tilde c,l,\tilde a'} \left\{ \tilde u(\tilde c,l)+\tilde v(\tilde g) +(1+g_A)^{\gamma (1-\eta)} \beta \phi^s_t \sum_{\theta'} \; prob(\theta'|\theta)\; \tilde V_{t+1}(\tilde z_{s+1})\right\},\end{equation} with the terminal condition $\tilde V_t (\tilde z_{T+1}) = 0$.

The budget constraint of the household in stationary variables is given by \begin{equation} (1+\tau^c_t) \tilde c = \tilde y +\left[ 1+ (1-\tau^k_t) (r_{t} -\delta) \right] \tilde a +\tilde tr - (1+g_A) \tilde a'\\ \end{equation} with \begin{equation} \tilde y = \left\{ \begin{array}{ll} (1-\tau^l_t-\tau^p_t) \epsilon(s,\theta,e)\, w_{t}\, l\;\; & s=1,\ldots, T^W,\\ &\\ \widetilde{pen}_t & s= T^W+1,\ldots,T. \end{array}\right. \end{equation}

Stationary Equilibrium

In the following, we consider the equilibrium in stationary variables. For the sake of simplicity, we assume that we have already discretized the state space. In particular, we have chosen the grid over the asset space $\tilde a\in {\cal A} = \{\tilde a_1,\tilde a_2,\ldots,\tilde a_{n_a}\}$ with $n_a$ grid points. Let $f(s,e,\theta,\tilde a)$ denote the (discretized) density function associated with the state $\tilde z=(s,e,\theta,\tilde a)$.

In order to express the equilibrium in terms of stationary variables, we have to divide aggregate variables by the product of aggregate productivity $A_t$ and
the mass of the total population $N_t$. Therefore, we define the following stationary variables, $\tilde X_t \equiv \frac{X_t}{A_t N_t}$, for the aggregate variables $X\in\{Pen, Tr, G, B, Beq, Tax, Y, K, C\}$. In particular, $\widetilde{Tr}_t = \widetilde{tr}_t$. Aggregate stationary labor is defined as $\tilde L_t = L_t/N_t$.

In equilibrium, aggregate effective labor supply is equal to the sum of the individual effective labor supplies: \begin{equation} \tilde L_t = \sum_{s=1}^{T^w} \sum_{i_\theta=1}^{n_\theta} \sum_{j=1}^2 \sum_{i_a=1}^{n_a} \,\epsilon(s,\theta_{i_\theta},e_j)\; l(s,\theta_{i_\theta},e_j,\tilde a_{i_a})\; f(s,\theta_{i_\theta},e_j,\tilde a_{i_a}). \end{equation} Aggregate wealth $\tilde \Omega_t$ is equal to the sum of the individual wealth levels: \begin{equation} \tilde \Omega_t = \sum_{s=1}^{T} \sum_{i_\theta=1}^{n_\theta} \sum_{j=1}^2 \sum_{i_a=1}^{n_a} \tilde a_{i_a} \; f(s,\theta_{i_\theta},e_j,\tilde a_{i_a}) . \end{equation} In capital market equilibrium, \begin{equation} \Omega_t = B_t + K_t.\end{equation} At the beginning of period $t+1$, the government collects accidental bequests from the $s$-year old households who do not survive from period $t$ till period $t+1$: \begin{equation} Beq_{t+1} = \sum_{s=2}^{T} \sum_{i_\theta=1}^{n_\theta} \sum_{j=1}^2 \sum_{i_a=1}^{n_a} (1-\phi^s_t) \left[1+(1-\tau^k_{t+1}) (r_{t+1}-\delta)\right]\, \tilde a'(s,\theta_{i_\theta},e_j,\tilde a_{i_a}) \; f(s,\theta_{i_\theta},e_j,\tilde a_{i_a},), \end{equation} where we have used the condition that the after-tax returns on the two assets are equal: \begin{equation} r^b_t = (1-\tau^k_{r}) (r_{t}-\delta) \end{equation}

In equilibrium, the goods markets clear: \begin{equation} \tilde Y_t = \tilde C_t + \tilde G_t + (1+g_A) (1+n) \tilde K_{t+1} -(1-\delta) \tilde K_{t},\end{equation} where aggregate consumption is the sum of individual consumptions: \begin{equation} C_t = \sum_{s=1}^{T} \sum_{i_\theta=1}^{n_\theta} \sum_{j=1}^2 \sum_{i_a=1}^{n_a} \tilde c(s,\theta_{i_\theta},e_j,\tilde a_{i_a}) \; f(s,\theta_{i_\theta},e_j,\tilde a_{i_a}). \end{equation}

In the following, we calibrate the steady state before we compute it.

Calibration

Our model is calibrated for the US economy for the year 2015. We assume that the US economy is in steady state in this period. Periods correspond to years. Agents are born at real lifetime age 21 which corresponds to age $s=1$. As stated above, they work $T^w=45$ years corresponding to a real lifetime age of 65 and live a maximum life of 70 years $(T^R=25)$ so that agents do not become older than real lifetime age 90. Our survival probabilities $\phi^s_t$ and population growth rates $n_t$ are taken from the United Nations (2015) which provides 5-year forecasts until the year 2100. We interpolate population data using cubic splines and assume that survival probabilities and the population growth rate are constant after 2100. The population date is stored to the Excel file 'Survival_probs.xlsx' which will be read into the PYHTON program 'AK70_stochastic_income.py'.

For the discount factor, we choose the parameter values $\beta=1.011$ in accordance with the empirical estimates of Hurd (1989), who explicitly accounts for mortality risk. In addition, we choose the intertemporal elasticity of substitution $1/\eta=1/2.0$. The preference parameter $\gamma=0.33$ is calibrated so that the average labor supply $\bar l$ is approximately equal to 0.30. We find this parameter by trying different values for $\gamma$. As a starting value, we fall back on the neoclassical growth model. For example, Heer (2019, Public Economics, p. 28) finds that $\gamma=0.338$ implies $\bar l=0.30$ in the neoclassical growth model. The model parameters are presented in following table:

Parameter Value Description
$\alpha$ 0.35 production elasticity of capital
$\delta$ 8.3% depreciation rate of capital
$g_A$ 2.0% growth rate of output
$n$ 0.75% population growth rate
$1/\eta$ 1/2 intertemporal elasticity of substitution
$\gamma$ 0.33 preference parameter for weight of utility from consumption
$\beta$ 1.011 discount factor
$\tau^l+\tau^p$ 28% tax on labor income
$\tau^k$ 36% tax on capital income plus social security tax
$\tau^c$ 5% tax on consumption
$G/Y$ 18% share of government spending in steady-state production
$B/Y$ 63% debt-output ratio
$repl$ 35.2% gross pension replacement rate
$\{$ $e_1,e_2$ $\}$ $\{$ 0.57, 1.43 $\}$ permanent productivity types

The tax parameters are set to $\tau^l+\tau^p=28\%$, $\tau^k=36\%$, and $\tau^c=5\%$, and the government consumption share in GDP is set equal to $G/Y=18\%$ following Trabandt and Uhlig (2011). The debt-output level amounts to $B/Y=63\%$. The gross replacement rate of pensions with respect to wage income (of the household with unitary individual efficiency and average labor supply $\bar l$), $repl=pen/(w\bar l)=35.2\%$, is taken from the OECD (2015).

Next, we calibrate the labor efficiency of the $s$-year old household, $\epsilon(s,\theta,e) = e \bar y^s e^\theta$. We choose the permanent efficiency types of the workers $(\epsilon_1,\epsilon_2)=(0.57,1.43)$ which reflects the wage difference between high school and university graduates. The mean efficiency index $\bar y^s$ of the s-year-old worker is taken from Hansen (1993) and is hump-shaped (with a peak around age 52). The idiosyncratic productivity shock $\theta$ follows a Markov process. The Markov process is given by: \begin{equation} \theta'=\rho \theta+\xi, \end{equation} where $\xi\sim N(0,\sigma_\xi)$ is distributed independently of age $s$. Huggett uses $\rho=0.96$ and $\sigma_\xi^2=0.045$. Furthermore, we follow Huggett (1996) and choose a log-normal distribution of earnings for the 20-year old with $\sigma_{y_1}=0.38$ and mean $\overline{y}^1$. As the log endowment of the initial generation of agents is normally distributed, the log efficiency of subsequent agents will continue to be normally distributed. This is a useful property of the earnings process, which has often been described as log-normal in the literature.

We discretize the state space $\Theta=(\theta_1,\ldots,\theta_{n_\theta})$ using $n_\theta=5$ values. The states $\theta_{i_\theta}$, $i_\theta=1,\ldots,n_\theta$, are equally spaced and range from $-m \sigma_{y_1}$ to $m\sigma_{y_1}$. We choose $m=1$ so that the Gini coefficient of hourly wages amounts to 0.374 in good accordance with empirical observations for the US economy. Our grid $\Theta$, therefore, is presented by: \begin{equation}\Theta = ( -0.7576, -0.3788, 0.0000, 0.3788, 0.7576 )\end{equation} The probability of having productivity shock $\theta_{i_\theta}$ in the first period of life is computed by integrating the area under the normal distribution implying the initial distribution among the 21-year-old agents for each permanent productivity type $e_i$, $i=1,2$: \begin{equation} ( 0.1783, 0.2010, 0.2413, 0.2010, 0.1783 ) \end{equation} Each permanent efficiency type $e_i$, $i=1,2$ has a share of 50% in each cohort. In our model, the distribution of productivity types $\theta$ in each cohort is constant and we will denote it by \begin{equation} \nu(\theta) = \left( \begin{array}{c} 0.1783\\ 0.2010\\ 0.2413\\ 0.2010\\ 0.1783 \end{array}\right) \end{equation}

The transition probabilities are computed using Tauchen's method as described in Algorithm 12.2.1 in Heer and Maussner (2009, DGE Modeling: Computational methods and applications, p. 658). As a consequence, the efficiency index $\theta$ follows a finite Markov-chain with transition matrix: \begin{equation} Prob(\theta'|\theta) = \left( \begin{array}{lllll} 0.7734 & 0.2210 & 0.0056 & 0.0000 & 0.0000\\ 0.1675 & 0.6268 & 0.2011 & 0.0046 & 0.0000\\ 0.0037 & 0.1823 & 0.6281 & 0.1823 & 0.0033\\ 0.0000 & 0.0046 & 0.2011 & 0.6268 & 0.1675\\ 0.0000 & 0.0000 & 0.0056 & 0.2210 & 0.7734 \end{array}\right) \end{equation}

results

Results

The above figure illustrates the Lorenz curves of earnings (green line) and wealth (orange line) which will result from our computation (final output). In steady state of our model, the (earnings-)poorest quintile of the workers receives only 3% of total earnings, while the top quintile earns 54% of total earnings. The Gini coefficient of earnings amounts to 0.50. Earnings are more concentrated than hourly wages as the more productive workers supply higher labor than the less productive workers in our model (as observed empirically).

Wealth is more concentrated than earnings. The (wealth-)poorest quintile holds no wealth at all (20% of the population are credit-constrained), while the top quintile of the wealth distribution holds 67% of total wealth. The Gini coefficient amounts to 0.66. Therefore, our model is able to replicate the fact that wealth is distributed more unequally than earnings. Our inequality of wealth falls a little short of empirical values because we neglect 1) bequests (Rockefeller's progeny), 2) entrepreneurship (Bill Gates), or 3) uncertain medical expenditures in old age in our model, among others.

Computation of the Steady State

The algorithm for the computation of the steady state consists of the following steps:

  1. Parameterize the model and choose asset grids for the individual state space.

  2. Make initial guesses of the steady state values of the aggregate capital stock $\tilde K$, labor $\tilde L$, mean working hours $\bar l$, labor income taxes $\tau^l$, the social security contribution rate $\tau^p$ and government transfers $\widetilde{tr}$.

  3. Compute the values $w$, $r$, and $\widetilde{pen}$, which solve the firm's Euler equation and the budget of the social security authority.

  4. Compute the household's decision functions by backward induction using value function iteration.

  5. Compute the optimal path for consumption, savings and labor supply for the new-born generation by forward induction given the initial asset level $\tilde a^1=0$ and distribution of idiosyncratic productivities $e$ and $\theta$.

  6. Compute the aggregate savings $\tilde \Omega$, labor supply $\tilde L$, mean working hours $\bar l$, aggregate taxes and transfers $\widetilde{tr}$.

  7. Update the aggregate variables and return to step 3 until convergence.

  8. Update the asset grid of the individual state space if necessary and return to step 3 until convergence.

Let us look at the program and anaylse the individual steps in detail.

Step 0: Setting up the PYTHON computer program

Please download the PYTHON code 'AK70_stochastic_income.py' and run it. In line 263, you specify the interpolation method. We set kind_tol = 'linear' for linear interpolation. If you would like to try cubic interpolation instead, just set kind_tol = 'cubic'. Notice that the computational time almost doubles with cubic interpolation.

The computational time amounts to approximately 29 hours on my computer (the GAUSS code is about 60 times faster). The solution for the aggregate variables are $\tilde K=1.596$, $\tilde L=0.310$, $\bar l = 0.305$, $\widetilde{tr}=0.0266$, $\widetilde{pen}=0.116$ and $\tau^p=11.6$%.

The program consists of an outer loop over the aggregate variables $(\tilde K,\tilde L,\widetilde{tr},\widetilde{pen}, \tau^p)$ and an inner loop that computes the policy functions as a solution to the value function iteration problem. After 29 iterations of the outer loop, the program execution stops when the aggregate capital stock $\tilde K$ and efficienct labor $\tilde L$ in two consecutive iterations diverge by less than 0.0001. In the inner loop, we just need one iteration since we can solve the finite lifetime value function problem backwards starting in the last period of life, $T=70$ and initialize the value function at this age as described above in the economic model.

In the following, I will present functions as they appear in the main part of the program. In the download version, we specified the functions in the beginning of the program instead so that it allows for easy editing.

In the first part of the program, we load the necessary libraries numpy, scipy, time, matplotlib and math. If you have not installed them or have difficulties to load them, please consult the above referenced manuscript of THOMAS J. SARGENT and JOHN STACHURSKI on 'Python Programming for Economics and Finance':

https://python-programming.quantecon.org/intro.html

In addition, we set the interpolation mode to 'linear'.

Step 1: Parameterization

Next, we set some numerical parameters. In particular, we set up an equi-spaced grid with $𝑛𝑎=501$ grid points over the individual asset space $\tilde \Omega = \{\tilde a_0=\tilde a^{min},\ldots,\tilde a_{n_a}=\tilde a^{max}\}$ with the minimum and maximum asset levels equal to $\tilde a^{min}=0$ and $\tilde a^{max}=20.0$. The upper boundary point $\tilde a^{max}$ is found with some trial and error so that the measure of households with wealth level equal and close to $\tilde a^{max}$ is zero. Since our algorithm restricts the optimal policy function $\tilde a'(.)$ to lie on the interval $[\tilde a^{min},\tilde a^{max}]$, the behavior of the policy functions displays abrupt changes at the upper boundary of the asset space interval and we want to restrict the evaluations of policy functions to $\tilde a \ll \tilde a^{max}$.

For the distribution, we choose a finer gird with twice as many grid points $nag=1002$ over the same interval (following the recommendation by Rios-Rull (1997)). In order to compute the Lorenz curves for earnings and income, we also specify two equispaced grids of $nearn=500$ points on the interval [0,5].

We will also define some constants for the tolerance of the final solution for the aggregate capital stock, $tol=0.0001$, and the tolerance of our optimization routine, 𝑡𝑜𝑙1=1𝑒−10 and set the maximum working hours equal to $lmax=0.60$. The number of iterations in the outer loop over aggregate variables is set equal to $nq=30$. In the following, however, we will just consider one loop for the ease of exposition.

The other parameters will be explained as we describe the individual parts of the program in more detail and can be ignored by the reader for now.

We calibrate the parameters of the model as presented in the table above. For this reason, we import the surival probabilities of the $s$-year old, $\phi^s$, $s=1,\ldots,70$, and the age-efficiency component, $\bar y^s$, $s=1,\ldots,45$, from the Excel file 'survival_probs.xlsx'. Be careful to download the Excel file 'survival_probs.xlsx' and save it in the same directory as the Python code 'AK70_stochastic_income.py'.

With the help of the stationary survival probabilities $\phi^s$ and the population growth rate $n$, we can compute the stationary age distribution in the population. Let $N_t$ denote total population and $N_t(s)$ the number of $s$-year old households. We define the measure of the $s$-year old in period $t$, $\mu^s_t$, as

\begin{equation} \mu^s_t \equiv \frac{N_t(s)}{N_t} \end{equation}

The sum of all measures $\mu^s$ is equal to one by definition: \begin{equation} \sum_{s=1}^{70} \mu^s_t = 1.0 \end{equation}

To compute the stationary measure $\{\mu^s\}_{s=1}^{70}$, we simply set up a vector with $nage=70$ entries, initialize the mass of the $1$-year old, $\mu^1$ equal to one (corresponding to the vector 'mass[0]'' in the program) and iterate over age $s=1,\ldots,nage-1$ as follows:

\begin{equation} \mu^{s+1} = \frac{\phi^s}{1+n} \mu^s \end{equation}

This formula is derived from the following two equations: \begin{equation} N_t(1) = (1+n) N_{t-1}(1) \end{equation} and \begin{equation} N_t(s) = \phi^{s-1} N_{t-1}(s-1) \end{equation} Division of both equation by population size $N_t$ and noticing that $1+n=\frac{N_t}{N_{t-1}}$ implies our formula.

After we have computed all $\mu^s$, $s=1,\ldots,70$, we normalize the measures so that their sum is equal to one. The measure of the cohorts in our calibration declines monotonically with age $s$:

Next, we also choose the discretization of the individual idiosycnratic productivity space $\Theta$ and the transition matrix of the Markov process, $prob(\theta'|\theta)$ as described in the previous section on the calibration.

In the following loop, we compute the Gini coefficient of the hourly wages. Notice that the total measure of the workers is equal to 'mass[0:nw]' and, therefore, the $s$-year-old worker has measure 'mass[s-1]/sum(mass[0:nw])' in the population of workers. Each permanent type has measure 1/2, while the worker with the idiosycnratic stochastic productivity component 'ye1[iy]' has measure 'muy[iy]'.

In order to compute the Gini coefficient of the wage distribution, we also have to order the wages in ascending order:

The Gini coefficient of wages amounts to 0.37 in good accordance with empirical data from the US economy. The Gini coefficient of wages is independent of the endogenous aggregate variables $\tilde K$ or $\tilde L$ since individual wages are simply proportional to individual productivity $\epsilon(s,\theta,e)$.

Step 2: Initialization

In this step, we need to provide initial guesses of the steady state values of the aggregate capital stock $\tilde K$, labor $\tilde L$, mean working hours $\bar l$, labor income taxes $\tau^l$, the social security contribution rate $\tau^p$ and government transfers $\widetilde{tr}$. We assume that households work 30% of their time, $\bar l=0.30$. We also simply assume that aggregate efficient labor is also equal to $\tilde L=0.30$. Of course, this is only an approximation in our economy. On the one hand, only 78\% of the households are working. Since we normalize the total measure of households to one, $\tilde L$ should be lower for this reason. On the other hand, we observe that workers with higher productivity also work longer hours. Therefore, $\tilde L$ should be higher than the number of working hours. As it turns out, the result of our computation is insensitive with respect to the initial guesses; however, computational time until convergence to the final values may depend on our initial choice.

Using a real interest rate equal to $r=3\%$ and efficient labor supply $\tilde L=0.30$, we can compute the implied steady-state value of the capital stock from the first-order condition of the firm. Therefore, we have to define the functions 'wage(k,l)' and 'interest(k,l)' which compute the two factor prices of capital 'k' and aggregate efficient labor 'l'. For later use, we also set up the Cobb-Douglas production function 'production(k,l)'. We can solve the first-order condition of the firm with respect to capital for capital 'kbar' (corresponding to $\tilde K$ in the model) implying 'kbar=1.708'.

In our model, aggregate capital ('kbar'=1.708) is lower than aggregate savings, $\tilde K< \tilde \Omega$, as part of the savings is invested into government bonds $\tilde B$. In equilibrium, government bonds amount to 63% of GDP. Since the wealth-GDP ratio is approximately 3.0 in the United States, government bonds amount to approximately 21% of wealth and, hence, physical capital is approximately 79% of wealth. Therefore, we can get an initial guess for our wealth $\tilde \Omega$ by using the approximation $\tilde \Omega\approx 1.2 \tilde K$.

As for the final remaining aggregate quantities, we need to provide a guess for the pension contribution rate $\tau^p$ and government transfers $\tilde tr$. Given $\tilde K$ and $\tilde L$, we can compute the wage rate $w$ and, hence, pensions with the help of the replacement rate, $\widetilde{pen} = repl \times w \, \bar l$. With the help of the social security budget, we can provide a guess for $\tau^p$ (noticing that the measure of the retirees is equal to 'sum(mass[nw:nage]'. Given that $\tau^l+\tau^p=28\%$ in the US economy, we can also compute the labor income tax rate $\tau^l$, which is approximately equal to 20% in our initialization. For government transfers, we simply use a small number as initial guess, 'trbar=0.01'.

Step 3: Computation of the Optimal Policy Functions with (Finite) Value Function Iteration

Next, we introduce the optimization function which we apply to the right-hand side of the Bellman equation in the value function iteration. We will use Golden Section Search (as described in Section 11.6.1 in our book). The function is called using the command 'GoldenSectionMax(f,ay,by,cy,tol)' where f is a one-dimensional function, the maximum is bracketed by the points ay and cy, by is an intermediate point and tol describes the maximum difference between the upper and lower boundary, cy-ay. In the above parameterization of the program, we have chosen the tolerance 'tol1=1e-5'.

We use a simple function 'testfunc(x)' with the maximum at point x=2.0 to test our function GoldenSectionMax(.).

Value Function Iteration: Retirees

In the inner loop, we compute the policy functions of the worker and retiree using value function iteration. First, we initialize the value function of the retirees using the command 'vr = np.zeros((na,tr))', where the row index refers to the asset level $\tilde a$ and the column index to the 25 retirement ages $iage=s-nw=0,\ldots,24$.

To compute the value functions at retirement age $iage=0,\ldots,24$, we start in the last period and iterate backwards in time (or age). The value function in the last period of the life, $T=T^W+T^R=70$ with given exogenous pensions $\widetilde{pen}$ and interest rate $r$ is presented by

\begin{equation}\tilde V^{T}(\tilde a^{T}) = u(\tilde c^{T},1) \end{equation}

with

\begin{equation} \tilde c^T_t = \frac{\widetilde{pen} + [1+(1-\tau^k)(r-\delta)] \tilde a^T-\widetilde{tr}}{1+\tau^c}. \end{equation}

Retired agents at age $T=70$ (corresponding to the iteration index $nr-1=24$ -- remember that Python starts indexing a vector with the number '0' so that a vector with 'nr' rows has the index 'nr-1' for the last row) consume their total income consisting of pensions and interest income. We compute this value function at all asset grid points $a_{i_a} \in \{a_0=0,a_1,\ldots,a_{501}=30\}$ and store them in the value function $vr[ia,24]$ for $ia=0,\ldots,501$. In order to compute the utility $u(c,1-0)$ for given consumption $c$ and leisure $1-l$, we use the function 'u(c,l)' with the arguments 'c' and 'l' for consumption and labor, respectively. Remember that in old age, labor supply is equal to zero so that leisure amounts to $1-l=1$. We also store the optimal consumption policy in the matrix 'cropt' and the optimal next-period wealth $a'$ (which is equal to zero in the final year of the life) in the matrix 'aropt'.

In the next step, we consider the retiree in his second-to-last period at retirement age $iage=23$. We use Golden section search to compute the optimal solution to the right-hand side of the Bellman equation:

\begin{equation} \tilde v^{23}(\tilde a^{23}) = \max_{\tilde a^{24}} \left\{ u\left(\frac{\widetilde{pen} +[1+(1-\tau^k)(r-\delta)] \tilde a^{23}+ \widetilde{tr}- (1+g_A) \tilde a^{24}}{1+\tau^c},1 \right) + (1+g_A)^{\gamma (1-\eta)} \beta \; \tilde v^{24} (\tilde a^{24}) \right\}. \end{equation}

This basic step is performed in the program part below in the line

 aropt[ia,iage-1]=GoldenSectionMax(value1,ax,bx,cx,tol1)

The right-hand side of the Bellman equation which we are maximizing with respect to the choice of next-period assets $\tilde a'$ is specified in the function 'value1(x)'. In order to compute the right-hand side of the Bellman equation at age $s$, we need to evaluate the value function at age $s+1$ for the asset level $\tilde a'=x$. Since x may not be a grid point that we stored in the matrix 'vropt', we need to use interpolation to find the value. We use linear inteprolation, but you may change the method into 'cubic' as described above. We, therefore, apply the command 'interpolate.interp1d(a,vr[:,iage], kind='linear')'.

As a second input into the Golden Section Search Routine, we need to bracket the maximum by the interval $[ax,cx]$ with an intermediate value $bx$. We, therefore, start to search over the grid points $\tilde a_{ia}$ until the value of the rhs of the Bellman equation starts to decline. Therefore, we first initialize the value 'v0' with the value 'neg=-1e10' and adjust it whenever we find a higher value for the rhs of the Bellman equation with the help of $\tilde a_{ia}$. Once the rhs of the Bellman equation starts to decline, we know that we have bracket the maximum. Remember that the rhs of the Bellman equation is a concave function of $\tilde a'$.

We also make use of the monotonocity of the value function and the policy function $\tilde a'(.)$. Since both functions increase monotonically with $\tilde a$, we also speed up our algorithm when we search for $[ax,cx]$. Assume that we have found the lower boundary $ax$ for the grid point $\tilde a_{i_a}$. When we start searching for the lower boundary $ax$ for the next grid point, $\tilde a_{i_a+1}>\tilde a_{i_a}$, we simply start at the value of $ax$ found in the last iteration.

In the following computation, we also need to take special care if $ax$ ($cx$) happens to be the boundary point $\tilde a^{min}$ ($\tilde a^{max}$). In this case, we have to evaluate if we have a corner solution $\tilde a'=\tilde a^{min}$ ($\tilde a'=\tilde a^{max}$). For example, we compare the rhs of the Bellman equation at $\tilde a^{min}$ with the one evaluated at $\tilde a^{min}+eps$, where 'eps' is a small constant. If the value of 'value1()' at $\tilde a^{min}$ is larger than at $\tilde a^{min}+eps$, we have a corner solution. Otherwise, we may apply Golden Section search again.

The value and policy functions of the 46-year old (first period of retirement) are displayed above. Notice that all three functions $\tilde a'(.)$, $\tilde c(.)$ and $\tilde v(.)$ are monotone functions and that $\tilde v(.)$ is concave.

Value Function Iteration of the Worker

Next, we turn to the value function iteration of the worker. The procedure is almost the same as in the case of the retiree. The worker also maximizes the right-hand side of the Bellman equation. However, different from the retiree, he also chooses optimal labor supply.

There are different ways how we can implement this two-dimensional optimization problem. We have choosen to break it up in two nested optimization problems. In the outer loop, we optimize over the next-period wealth $a'$ just like in the case of the retiree. The Golden Section Search algorithm maximizes the value function 'value2(x)' of the worker where the argument $x$ of the function is the next-period wealth level $\tilde a'$. In the inner loop (inside the function 'value2(x)'), we compute the optimal labor supply given this and next-period wealth, $\tilde a$ and $\tilde a'$ using the first-order condition of the worker with respect to his labor supply:

\begin{equation} \frac{(1-\tau^l-\tau^p) \epsilon(s,\theta,e) w}{1+\tau^c} = \frac{1-\gamma}{\gamma} \frac{\tilde c}{1-l}. \end{equation}

After the substition of $\tilde c$ from the budget constraint of the worker, we can solve for the labor supply $l$ of the worker:

\begin{equation} l = \gamma - \frac{1-\gamma}{(1-\tau^l-\tau^p) \epsilon(s,\theta,e)w} \left( \left[ 1+(1-\tau^k)(r-\delta) \right]\tilde a +\widetilde{tr}-(1+g_A) \tilde a' \right) \end{equation}

We also have to impose the constraint $l\ge 0$ and $l \le 0.60$ if $l$ happens to lie outside the admissible space.

Notice that, in the present case, it is easy to compute the optimal labor supply since we can solve for $l$ explicitly. If we use a different utility function or consider contributions-based pensions, labor supply may only be computed implicitly and we have to use more sophisticated methods. In the case of a different utility function, we may have to solve a non-linear equation problem. If we also consider contributions-based pensions, the computation of the optimal labor supply is more complicated and we will discuss it in the next section of our book.

The implementation of this nested optimization is straightforward. The function 'value2(x)' computes the right-hand side of the Bellman equation for next-period wealth $\tilde a'=x$. The present period wealth is given by the global variable 'k0'. We also use the productivity type $\epsilon(s,\theta,e)$ as given by the global variables 'perm[iperm]', 'ye1[iy]' and 'ef[s]'. With the help of these variables, $\tilde a'$ and $\tilde a$, and the factor prices $w$ and $r$, we are able to compute the optimal labor supply $l$. We also have to check if $l\ge 0$ and $l\le 0.60$. If not, we set it either equal to zero, $l=0$, or $l=0.60$. In order to compute the Bellman equation, we also have to determine the value of the next-period value function at age $i+1$ and wealth $x=\tilde a'$. We again interpolate either linearly or with a cubic spline. Notice that at age $s=45$ (corresponding to $iage=44$ in the iteration over the age of the worker), the next-period value function is the value function of the retiree during his first period of retirement. In this case, we have to interpolate 'vw0=vr[:,0]''; otherwise, the next-period value function is given by 'vw0=vw[:,iage+1]''.

In the outer loop, we just use golden section search with respect to next-period wealth $\tilde a'$. This is implemented in the line

'k1 = GoldenSectionMax(value2,ax,bx,cx,tol1)'

Again, we have to bracket the maximum by the interval $[ax,cx]$ and have to take special care with the boundary points $\tilde a^{min}$ and $\tilde a^{max}$ just like in the case with the retired agent.

The value function of the worker at ages $iage=44,\ldots,0$ is computed in the next part of the program:

It is always good idea to plot the policy function for some different productivity types and ages (remember that the number of all productivities and ages of the workers alone amounts to $nperm \times ny \times nw = 450$). Therefore, we need to restrict our attention to the study of a random choice. Above, we graph the labor supply function of the skilled and unskilled worker at age $s=21$ with mean idosyncratic productivity $iy=3$ (corresponding to $\theta=0.378$). Labor supply falls with higher wealth $\tilde a$. The high-skilled (orange curve) has a much higher labor supply than the low-skilled (blue curve) for given wealth $\tilde a$ (but since the high-skilled holds higher wealth, the difference in observed working hours is smaller). In addition, the lower boundary $l\ge 0$ starts binding for the low-skilled worker for a wealth level in excess of 12.0.

Euler equation residual

To evaluate the accuracy of our policy function approximation, we study the so-called 'Euler residual' which is defined as the percentage deviation of the Euler equation. (The Euler residual is considered in many chapters of our book on 'Dynamic General Equilibrium Modeling', for example in Chapter 6 where we consider projection methods.) For the $s$-year old worker with wealth level $\tilde a$ and productivity type $\epsilon(s,\theta,e)$, the Euler equation residual is defined by:

\begin{equation} R(\tilde a) = 1 - \frac{\tilde u_c(\tilde c,l)}{\beta \phi^s \mathbb{E} \left\{\tilde u_c(\tilde c',l')\right\}}, \end{equation}

where $\tilde c'$ and $l'$ are next-period consumption level and labor supply and expectations $\mathbb{E}$ are conditional on information at age $s$ in period $t$ (after observing the productivity shock in this period).

If $R(\tilde a)=0$, the Euler equation error is zero and we have a perfect fit. The Euler equation residual is a popular measure of accuracy. In particular, Watson (2000, 'Accuracy of Numerical Solutions using the Euler Equation Residuals', Econometrica) has shown that the accuracy of the numerical approximation of the policy function is of the same order of magnitude as the numerical error of the Euler equation residual.

Therefore, we evaluate the Euler equation error in the following to get an idea about the goodness of fit. We use a finer grid of the asset level for the Euler equation residual than that for the policy function. We want to consider points off the policy function grid where we have to interpolate between grid points and accuracy might be smaller. As our measure of fit we use the average absolute Euler equation residual.

In order to compute the Euler equation error, we also need to define functions for the marginal utility of consumption for the present and the next period, $uc(.)$ and $uc1(.)$ in our Python code. The evaluation of next-period marginal utility is very time-consuming because we need to interpolate the next-period policy functions (stored in 'cwopt' and 'lopt') to derive $c'$ and $l'$. Therefore, you should compute the Euler residual only when you develop the program to check for the accuracy of your computer code and after the final outer loop over the aggregate variables when you want to determine if the number of grid points and the mode of interpolation is sufficient.

The Euler equation residual is equal to 0.11% and 0.26% for the workers and retirees, respectively. Given the accuracy of our Golden Section Search routine (1e-5), this accuracy is admissible and we are pretty certain that we have programmed the value function iteration problem correctly. The solution to the dynamic programming problem and from the Euler equation coincide. If you want to increase accuracy, you need to 1) consider more grid points $na$ and/or 2) use cubic interpolation instead of linear interpolation. As often in computational economics, you are faced with the trade-off between speed and accuracy.

Step 5: Distribution Function

In Step 5, we compute the endogenous wealth distribution in each cohort over an equispaced grid of the asset space $[\tilde a^{min},\tilde a^{max}]\times \Theta$ with the number of grid points \begin{equation} n_{ag} \times n_{e}\times n_\theta \times nw + n_{ag} \times nr = 1202 \times 2 \times 5 \times 45 + 1202 \times 25 = 570,950. \end{equation}

The distribution is stored in the variables 'gkw' and 'gk' for the worker and retiree, respectively. In the case of the worker, 'gkw' is a four-dimensional object, while, in the case of the retired household, 'gkr' is only two-dimensional because the productivity type does not affect the behavior in old age.

We start with the newborn generation at age $s=1$ with zero wealth. Furthermore, we know the distribution of the idiosyncratic productivity at age 1 which is given by $\nu(\theta)$ (see calibration section). Each permanent productivity has measure 1/2 in each cohort of the workers. The distribution is initialized in the following part of the program:

We also initialize the variables 'kgen', 'gwealth', 'gearn' and 'gincome' which store the wealth of the $s$-year-old and the distribution of wealth, earnings and income.

Given the distribution of the wealth level $\tilde a$ and the productivity $\epsilon(s,\theta,e)$ at age $s=1$, we can compute the distribution of $(\tilde a, \theta, e)$ at age $s=2$ by using the optimal decision functions of the agents with respect to labor supply $l$, consumption $\tilde c$ and next-period wealth $\tilde a'$ and the transition probabilities for the idiosyncratic productivities. The policy functions at age $s=1$ are stored in the variables 'lopt[:,:,:,0]', 'cwopt[:,:,:,0]' and 'awopt[:,:,:,0]'.

For expositional reasons, let us look at a specific example, e.g. the low-skilled worker with the idiosyncratic productivity $\theta_4=0.3788$. His measure is equal to 0.00213 and given by the product of 1) the measure of the 1-year old, $\mu^1$, the share of workers with idiosyncratic productivity $\theta_4$, $\nu(\theta_4)$, and 3) the share of low-skilled workers among the workers (equal to 1/2). His optimal next-period wealth is presented by 'awopt[iag,iy,iperm,iage]' with 'iag=0', 'iy=3', 'iperm=0' and 'iage=0':

His measure is given by

Therefore, we know that the households at age $s=1$ in period $t$ with measure equal to 0.00213 will choose to have next-period wealth $\tilde a'=0.008365$. Some households will die so that only $\phi^1$ of the 1-year-old households survive. In addition, the share of each cohort decreases between period $t$ and period $t+1$ by the factor $1/(1+n)$ because the population grows. As a consequence, the measure of household in the next period that is implied by this type of worker amounts to \begin{equation} 0.002129 \frac{\phi^1}{1+n} = 0.002111. \end{equation} Similarly, we compute the dynamics of the households over the complete productivity and asset space. In the next iteration, we increase the age by one.

We also have to discuss how we handle the case when the next-period asset $\tilde a'$ is not a grid point on the grid 'ag' for the distribution function. For example, in our example above, the next-period wealth amounts to $\tilde a'=0.008365$ and lies between the first two grid points $ag_1=0$ and $ag_2=0.01998$. In this case, we proceed as described in Chapter 7.2 of our book on 'Dynamic General Equilibrium Modeling'. We introduce a simple lottery. If the optimal next-period wealth $\tilde a'$ happens to lie between $ag_{ia-1}$ and $ag_{ia}$, $ag_{ia-1}<\tilde a'<ag_{ia}$, we simply assume that the next-period wealth level will be $ag_{ia}$ with probability $(\tilde a'-ag_{ia-1})/(ag_{ia}-ag_{ia-1})$ and $ag_{ia-1}$ with the complementary probability $(ag_{ia}-\tilde a')/(ag_{ia}-ag_{ia-1})$. This lottery allocation is computed in the last 5 lines of the following program part:

In the above program part, we have also compute the distribution of earnings and income and the wealth held by the $s$-year old cohort with the help of the measure 'gkw[iag,iy,iperm,iage]' where 'iag', $iag = 1,\ldots,nag$, denotes the index of the asset grid point, 'ag[ig]', 'iy', $iy = 1,\ldots,ny$, denotes the index of the idiosyncratic productivity type, 'ye[iy]', 'iperm', $iperm=1,2$, denotes the index of the permanent productivity type, 'perm[iperm]', and 'iage' denotes the age of the household. (Gross) income is defined as the sum of earnings and interest income (but excluding transfers).

In the last period of working age, ($s=45$ in the model, 'iage=44' in the program code), we need to take care of the fact, that the distribution of wealth among the retired is not stored for the different productivity types but only depending on age $s$. Therefore, the measure of all households with optimal next-period wealth $\tilde a'=a_{iag}$ at age 'iage=44' is added to the measure 'gkr[iag,0]' of the retired households at their first year of retirement for all productivity types, 'iperm' and 'iy' of the worker.

In the final year as a worker and during retirement, we continue to iterate over the dynamics of the distribution as follows:

In our computation, the aggregate wealth of the households during their first year of retirement amounts to 0.0361. The average wealth of the 46-year old is presented by

In the following, we compute the Gini coefficients of the wealth, earnings and income distribution:

Notice that we can already confirm the result of Huggett (1996) in this first iteration where the aggregate variables are not in steady state. Given the empirival distribution of wages, the OLG model implies a more unequal distribution of income and earnings than that of wages, and wealth is distributed even more unequally than earnings and wealth.

Next, we update the aggregate variables $\tilde \Omega$ which is simply equal to the sum of all savings. To derive the new value of total wealth for the next loop over the aggregate variables, we simply use an average of the old and new value with weights 'phi=0.8' and '1-phi=0.2'. To compute the new aggregate capital stock $\tilde K$, we need to compute production $\tilde Y$ first (using the old values of the capital stock and aggregate labor) and use the calibration that debt $\tilde B$ is equal to $0.63 \tilde Y$ to compute $\tilde B$. Aggregate capital $\tilde K$ is implied by the capital market equilibrium: \begin{equation} \tilde K = \tilde \Omega -\tilde B \end{equation} Using this new value for the capital stock ('knew'), we update 'kbar' (again, using a weighted average of the new and the old value of the capital stock) to derive 'kbar=1.618'.

We also need to compute the aggregate values of labor supply $\tilde L$, bequests $\widetilde{Beq}$, consumption $\tilde C$ and the mean labor supply $\bar l$. For this reason, we simply add up the individual variables given the distribution 'gkw' and 'gkr'. We have stored the optimal policy functions of the worker in the variables 'lopt[ia,iy,iperm,iage]' and 'cwopt[ia,iy,iperm,iage]' for labor supply and consumption, respectively. Given that the grid points 'ag[iag]' of the distribution are finer than the grid points of the policy function 'a[ia]', we need to interpolate to find the optimal policy function at grid points 'ag[iag]'. (The reason we choose a finer grid for the distribution than for the policy function is motivated in our book and follows a recommendation of Rios-Rull (1999).)

One reason to compute the distribution of individual wealth using the variable 'gwealth' consists of its latter use in the Step 8 where we decide if the asset grid is chosen correctly. In particular, we would like to choose an asset grid that is wide enough so that the upper boundary $\tilde a^{max}$ does not bind and households do not hold maximum wealth (where the approximation of the policy function is poor). When we study the measure of the upper grid points in 'gwealth' we find that they are all equal to zero. Similarly, the grid is not too wide so that we also observe agents who hold wealth level in the upper quintile of the asset space. Therefore, the grid is chosen wisely. Otherwise, you need to adjust the grid.

Finally, we have to update the aggregate variables, again using a weighted average from the starting value and the present iteration. We also compute aggregate taxes which are the sum of labor income taxes, capital income taxes and consumption taxes (for this reason, we have computed aggregate consumption $\tilde C$). With the help of the steady-state government budget constraint, we can compute the new values of transfers: \begin{equation} \widetilde{tr} = \widetilde{Tr} = \widetilde{Tax} + \widetilde{Beq}+ [(1+g_A)(1+n)- (1+r^b)] \tilde B- \tilde G, \end{equation}

where aggregate government consumption $\tilde G$ is computed with the help of the calibration parameters and amounts to $0.18 \tilde Y$.

The new pensions are computed with the help of average working hours and the replacement rate \begin{equation} \widetilde{pen} = repl \; w \bar l \end{equation}

With the help of the social security budget, we can compute the new value of the social security tax (pension contribution rate) $\tau^p$: \begin{equation} \tau^p = \frac{\widetilde{Pen}}{w \tilde L} \end{equation}

We save and plot the policy functions and distributions:

The three age profiles behave as expected. The (average-)wealth-age profile is hump-shaped and peaks prior to retirement. Households save for old age when their non-wealth income shrinks (pensions are below earnings). Consumption peaks in the last period of the working life. At the first year of retirement, leisure jumps to 1.0 (retirees do not work) and, in order to smooth utility over time, consumption falls. Working hours (averages among the $s$-year-old) peak at age 32 and decline thereafter. The initial increase of the labor supply is caused by the increase in the age productivity $\bar y^s$; however, $\bar y^s$ peaks at age 52 while labor starts to decline prior to this age because of the wealth effect. Workers get wealth-richer with increasing age and, for this reason, reduce their labor supply ceteris paribus.

This terminate the inner loop and we continue to iterate over the aggregate variable in the outer loop of the program 'AK70_stochastic_income.py' until convergence, which takes place after 29 iterations.

This version: May 12, 2021