Chapter 7: Computation of Stationary Distributions --- Python Code

Section 7.2: The Stationary Equilibrium of a Heterogeneous-Agent Economy

In the following, we study the computation of the stationary distribution function of a simple Ramsey model described in Section 7.1 of our book on 'Dynamic General Equilibrium Modeling: Computational Methods and Applications'. It is based upon the simple neoclassical growth model. However, households are subject to an exogenous employment shock and are either employed or unemployed. Depending on their employment history, they accumulate different levels of wealth. We compute the stationary distribution of wealth among both the employed and unemployed households. In our computation, we follow the Algorithms 7.2.1 and 7.2.3 from the book.

Let us briefly restate the main equations of the model. Households maximize their intertemporal utility

$\nonumber E_0 \sum_{t=0}^\infty \beta^t u\left(c_t\right)$

subject to the budget constraint

$a_{t+1}=\left\{ \begin{array}{ll} \left(1+(1-\tau_t) r_t\right) a_t+ (1-\tau_t) w_t-c_t, &\mbox{ if } \epsilon_t=e,\\ \left(1+(1-\tau_t) r_t\right) a_t + b_t-c_t, &\mbox{ if } \epsilon_t=u, \end{array} \right. $

Employed workers work one unit and receive net wage $(1-\tau_t) w_t$, while unemployed workers receive unempoyment compensation $b_t$.

The employment status follows a Markov process which is characterized by the employment transition matrix $\pi$ between the state of employment ('e') and unemployment ('u')

$\pi(\epsilon'|\epsilon)= Prob\left\{\epsilon_{t+1}=\epsilon'|\epsilon_{t}=\epsilon\right\}= \left(\begin{array}{ll} p_{ee} & p_{eu}\\ p_{ue} & p_{uu} \end{array} \right)$

In equilibrium, the aggregate consistency condition holds

$K_t=\sum_{\epsilon_t\in \{e,u\}}\int_{a_{min}}^\infty a_t f_t(\epsilon_t,a_t)\; da_t,$

where $f_t$ denotes the density function corresponding to the distribution functions $F_t$. In stationary equilibrium, $f_t(.)=f(.)$ is constant and the fiscal budget of the government balances,

$ \tau_t (w_t N_t + r_t K_t) = b_t (1-N_t)$

The aggregate consistency condition with respect to the labor market and the descripton of the production sector are presented in more detail in Section 7.1 of our book.

The basic steps in the computation of the stationary distriution $f()$ consists of an outer and inner loop. In the outer loop, we compute the aggregate capital stock and, hence, the factor prices and unemployment compensation. In each iteration, we find a new value of the stationary capital stock and simply update the old value by 5%. For higher values of the updating parameter, we find that the algorithm may not converge.

In the inner loop, we first compute the optimal policy function of the household using value function iteration. The solution provides the optimal next-period asset a' for present-period asset a and employment status e. Next, we simulate the economy for a given initial distribution of assets. We use a very simple initial distribution. Everyone holds the same (=average) wealth. Given the optimal policy function, we can compute the next-period distribution of assets and employment status. We continue to simulate the economy over a long time horizon of n periods so that the initial distribution does not affect the long-run distribution. The invariant long-run distribution is called the stationary distribution.

In the following, we describe the inner loop in more detail; the code of the outer loop only contains one line for the update of the aggregate capital stock. The computer code 'RCh7_denf.py' can be downloaded from our homepage

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

It translates our GAUSS code 'RCH7_denf.g' line by line and the number of arithmetic operations are identical for these programs.

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 Sargent and Starchurski.

In the following, we will define functions that we are using in the program and select parameter values. In the program 'RCh7_denf.py', these two parts are separated so that it allows for an easier editing. In order to ease the exposition and help you to learn about the program, I will switch between these two parts in the following.

We will start with the definition of the optimization function which we apply 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.

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

In the next step, we compute the ergodic distribution of employment. The transition is described by a Markov process with transition matrix

$\pi(\epsilon'|\epsilon)= Prob\left\{\epsilon_{t+1}=\epsilon'|\epsilon_{t}=\epsilon\right\}= \left(\begin{array}{ll} 0.9565 & 0.0435\\ 0.5 & 0.5 \end{array} \right)$

For example, the probability to remain employed (become unemployed) next period if you are employed this period amounts to 95.65% (4.35%). The unemployed remains unemployed with a probability of 50%.

We compute the ergodic distribution (0.92,0.08) with the help of the function 'equivec(p)'. The computation of the ergodic distribution is described in more detail in Section 12.2 of our book.

In the next step, we define the functions for utility, the wage and the interest rate. Therefore, we also need to calibrate the parameters of these functions as described in the Section 7.1 of the book.

In addition, we set a replace rate of unemployment compensation (with respect to the net wage) equal to 25% and an inital interest rate of 0.5% (given a period length of 1/8 of a year).

We also have to set a couple of numerical parameters. One set of parameters is on the coarseness of the grid over the individual asset space. Another set of parameters specifies the accuracy of our solutions on the value function iteration and the convergence of the distribution. We also specify the grids over the asset space. The asset grid 'a' is used for the value function and optimal policy function, the coarser grid 'ag' is used for the distribution function. All parameters are described as comments in the program.

In the first reading of this tutorial, you may skip this section; we will refer to these parameters in the following individual parts of the optimization problem and simulation of the distribution dynamics whenever appropriate.

The employment rate is equal to 92% in the population which is normalized to one. Each employed household supplies one unit of labor so that aggregate employment is also equal to 0.92.

We initialize the aggregate capital stock K with the steady state capital stock of the corresponding regresentative agent model (with exogenous labor equal to 0.92) and compute the factor prices of labor and capital and unemployment compensation.

The aggregate capital stock amounts to kk0=247.62, the wage rate w is equal to 4.797, and unemployment compensation b is set equal to 1.199.

The next big step consists of computing the optimal policy functions of the households with wealth a and employment status e (in {0,1}) using value function iteration. You find a description of this method in Sections 4.1 and 4.2 of our book.

In a first step, we need to come up with an intial guess for the value function. The iteration over the value function converges faster, if we use a good initial guess. One natural guess is to use the values for the value function that are implied by a household behavior where they simply consume their current income forever. In this case, utility is constant in each period and the value function is given by the sum over a geometric series with discount factor $\beta$. The initial guess is stored in the variable 'v' which has dimension (na x 2).

For example, we initialize the value function of the employed with wealth level a=103.07 with the value -38.396.

In the value function iteration, we will use the variable 'vold' for the value function in iteration j and 'v' in the next iteration j+1. In the next iteration, 'v' will become 'vold' and we search for the new value 'v'.

In our maximization problem, we will have to evaluate the value function 'vold' at points a which do not necessarily lie on a grid point. Since we store the value function only at grid points, we need to interpolate. This is done using the function 'value(x,y)' where 'x' and 'y' specify the asset level and the employment status, respectively. For interpolation, we use the library scipy and the command 'interpolate.interpolate1d()'. The following text code sets up this part of the code. As an example, we compute the value function of the unemployed (y=1) and asset level x=100 equal to v(100,1)=-118.257.

We need to define two more functions before we can start with the solution of the value function problem.

The first function evaluates the right-hand side of the Bellman equation. Therefore, we need to provide the present and next-period asset levels of the household, 'a0' and 'a1', as well as her present employment status 'y'. The function 'bellman(a0,a1,y)' returns the sum of the present utility and expected (discounted) next-period value. Therefore, it has to evaluate the value function at the next period asset level at point a1. Since we only store the value function at grid points and a1 might not be a grid point, we have to interpolate the value function with the help of the function 'value(.)'.

Notice that, in the function 'bellman()', we have to take care of the cases where consumption is negative so that the function utility(c) would result in a breakdown of the code (exponent of a negative function). In addition, we need to take care of the case when next-period wealth is above the upper grid point (amax) for the asset grid. In our solution, these corner solutions are irrelevant.

For values of this-period wealth, a0=90, next-period wealth, a1=80, and the employed household, y=0, the right-hand side of the Bellman equation amounts to -42.861.

Finally, we have to define an auxiliary function. We want to maximize the right-hand side of the Bellman equation 'bellman(a0,a1,y)'. This function has 3 arguments. In our optimization routine 'GoldenSectionMax', however, we need to pass on a function 'f()' with one argument only. Therefore, we define a new function 'value1(a1)' in the variable 'a1' only. The variables 'a0' and 'y' are specified as the global variables 'ainit' and 'e' in the main part of the program, respectively.

In the next part, we perform the iteration over the value function for given aggregate factor prices w and r. The iteration stops if either the number of iterations j exceeds nit-1=99 or if the mean difference between the entries in two successive value function iterations, 'vold-v', is below tol=0.001.

The computational time of the value function iteration amounts to 48 seconds on my computer. The solution of the value functions are plotted in the figure above for the employed (blue line) and unemployed (red line). For large values of the individual wealth 'a' on the abscissa, the two value functions almost coincide. At low values of the asset level around a=0, the two households have noticeable different values.

Let us discuss the individual parts of the value function iteration above:

aopt[i,e] = GoldenSectionMax(value1,ax,bx,cx,tol1)

In this line, we compute the optimal next-period asset a' by maximizing the right-hand side of the Bellman equation:

$ V(a,e) = \max_{a'} \frac{c^{1-\sigma}}{1-\sigma} + \beta \left[ p_{e0} V(a',0) + p_{e1} V(a',1)\right]$

where $p_{e0}$ and $p_{e1}$ denote the transition probabilities from employment status e in period t to the employed and unemployed status in period t+1 and consumption follows from the budget constraints of the employed

$ c = (1-\tau) w + (1+(1-\tau) r) a - a'$

and unemployed

$ c = b + (1+(1-\tau) r) a - a'$

The optimal policy is stored in the variable aopt[a,e]. In order to perform the Golden Section search, we have to bracket the maximum. Therefore, we allocate the two grid points ax and cx which contain the maximum. Let 'bx' denote the midpoint. Assume we have allocated the points ax, bx, cx. We use golden section search to find the maximum of the function 'value1(a')' which evaluates the Bellman equation at any point a' for given values ainit and e, the present-period asset level and employment status. The values ainit and e are assigned at the beginning of the two inner loops.

The right-hand side of the Bellman equation is computed in the following line of the function bellman(ainit,a1,y):

'utility(c) + beta(prob[y,0]value(a1,0) + prob[y,1]*value(a1,1))'

i) We start the computation of the value function for the agent with wealth a[0] and employment status e=0 (employed). His value function is initialized with the large negative value 'v0=negvalue'. Next, we evaluate the rhs of the Bellman equation at point a[0]=amin. The value of the value function for a'=amin is computed as 'v1'. If v1 is larger than v0, we store the new value in the variable v0. Otherwise, we finish the iteration over the next-period asset level (because the value function is concave and monotone). We contiune like this until we have found the grid point a[l] that maximizes the rhs of the Bellman equation. We store the maximum grid point in the variable 'bx' and the adjacent points in the variables ax and cx. If the next period asset level on the grid a that maximizes the Bellman equation is amin, our program above implies that ax=bx. If the optimal next-period wealth is found to be amax, we have the condition bx=cx. Otherwise, we have an interior grid point, which is the standard case.

In addition, we are exploiting the monotonocity of the value function with respect to its argument a. If a increases, we know that a' also increases. Therefore, if we have found the maximum of the Bellman equation to be a[l] for the household with wealth level a[i], we know that the optimal next-period wealth level of the household with wealth a[i+1] must be at least a[l]. Consequently, we start our search for the next-period optimal wealth at this point and use the index l0 for the starting point of the search. Notice that, therefore, we iterate over a and not e in the most inner loop to economize on the computational speed and exploit concavity of V(a,.).

ii) To handle the two cases where the maximum next-period grid point on a is either a[0]=amin or a[na-1]=amax, we compare the value of the rhs of the Bellman equation for 'amin' with 'amin+eps' and for 'amax' with 'amax-eps', where eps is a small constant. If amin (amax) is the larger of the two, the optimal next-period wealth is the boundary point amin (amax). If not, we search for the maximum in the interval [ax=a[0],cx=a[1]] with bx=a[0]+eps=amin+eps (or in the interval [ax=a[na-2],cx=a[na-1]] with bx=a[na-1]-eps=amax-eps.)

Finally, we compute the difference of the old and the new value function and assign it to the variable crit. We update our value function. Our computed value function v becomes the value function 'vold' in the next iteration j+1:

'vold = v + 0'

Notice that, in Python, we have to add the '0' on the rhs of this assignment. Otherwise, vold would be adjusted each time we assign a new value v[i,e] in our iteration over the value function. This is a special feature of Python which is absent from most other computer languages. If you are used to program in other languages like Gauss or Matlab, you need to be careful! (The second particularity of Python which often causes coding errors for people like me who are used to other computer languages are the indexing of the vectors, which starts with a[0] and not a[1] and ends with a[na-1] rather than a[na].)

For given policy function 'aopt(a,e)' for the next-period wealth given present wealth 'a' and employment status 'e', we are able to simulate the dynamics of the distribution function (keeping the production factors K and N and, hence, the factor prices w and r fixed) for given initial distribution. We will use the equal distribution where all agents hold the same wealth (equal to the aggregate capital stock). This part of the computer program is described by the following code lines:

The time for the computation of the distribution dynamics over 500 periods amounts to 14 minutes 9 seconds on my computer. Notice that this time is prohibitive. In the final iteration q in the outer loop over the aggregate capital stock K, we will increase the number of periods to $n=25,000$. We will comment on this separately in the conclusion to this application. Let me mention that other computer languages including Gauss, C++, Fortran or Julia are much faster with respect to the present application. For example, the computational time in Gauss amounts to much less than a minute. Python is particularly slow in large loops.

In the upper figure, we illustrate the density function of the assets for the employed (blue line) and unemployed (orange line) at the end of the transition in period 500. In the bottom figure, the mean asset level of the households during the transition is graphed. We notice that it does not converge yet (because we are not close to the stationary capital stock and the value function still needs to get more accurate). At the end of the transition, we find the capital stock to amount to kk1=374.907 (correspondig to $K_t$ in the model). We use this value to update our new guess of the capital stock in the next iteration q+1 over the aggregate capital stock.

Notice that we started with a low number of transition periods equal to 500. If we had started with a much larger number (like 25,000 in the final iteration), we would have approached the upper limit on the asset grid a. The reason is as follows: We start with a relatively low value for K (as we will find out shortly). For this reason, the interest rate r is higher than in the stationary equilibrium and households save more. Consequently, over time, average wealth increases above the stationary value of K.

If we had started with another guess for the initial capital stock K much higher than the stationary capital stock (so that the interest rate is lower and people save less), we might have approached the lower limit on the asset grid. As a consequence, if we start with a number of transition periods that is too high we might end up jumping between amin and amax as our new values of the capital stock K after the transition. For this reason, we slowly increase the number of transition periods from 500 to 25,000 in the outer loop.

Let me describe the individual elements of the iteration of the distribution function above in turn:

We also take care to handle the case where k0 happens to lie on the boundaries of the asset space (at amin or amax).

For the employed household with ag[17]=253.17, for example, optimal next-period wealth amounts to k1=aopt[17,0]=253.72.

In these six lines, we look at each measure of households with wealth ag[i] and employment status e and look at the resulting change in next-period density. For example, the employed households with wealth ag[17]=253.17 have measure 0.92 in period q1=0. In the next period, q1=1, their capital stock k1 amounts to 253.72. In order to allocate the measure of households equal to 0.92 to the points of the grid ag in the next period, we have to find out the adjacent points of k1=253.72. In particular, k1 lies between ag[17]=253.72 and a[18]=268.18 as implied by the first line above, j=18. In order to allocate the measure of households to the two grid points ag[17] and ag[18], we allocate the share 1-n0=96.3% of the total measure to the grid point ag[17] and n0=3.7% to the grid point ag[18]. The weighted average, (1-n0)ag[17]+ n0ag[18] is equal to the next-period capital stock k1 (see Section 7.2 in our book). In addition, we have to take care of the fact, that the employed households remain employed with probability prob[0,0]=0.9565 and become unemployed with the complementary probability prob[0,1]=0.0435. After this step, we therefore end up with the new distribution

gk[17:19,:] = array([[0.8476678, 0.0385505], [0.032277 , 0.0014679]])

After we have also treated the unemployed agent with wealth ag[17] in the same way, we have finished the computation of the distribution gk in period q1=1:

gk[16:19,:]=array([[0.00759829, 0.00759829], [0.88008791, 0.07097061], [0.032277 , 0.0014679 ]])

All other entries of the distribution matrix 'gk' are zero after the first iteration.

In the computer code above, we again also take care of the cases where k1 is equal to the boundaries of the interval over the asset space a, amin or amax.

The distribution that is computed after the complete transition in period q1=ngk=500 is illustrated above.

Finally, we have to update the aggregate capital stock. In our simulation of the distribution dynamics, we have found a new capital stock kk1=374.91 in the last period of transition, $q1=500$. We update the old value kk0 using a weighted average of the new and the old capital stock. We only slowly adjusts the aggregate capital stock (by 5%) in order to ensure convergence.

In the second outer iteration, we start with a new value of the aggregate capital stock kk0=260.03. We continue to iterate over K until convergence. After $q=51$ iterations, the program converges to the stationary capital stock $K=247.3$.

In conclusion, let us comment on the computational speed. We have found that the loop over the dynamics of the distribution function is too time-consuming. Therefore, in the program 'RCh7_denf.py', we decreased the number of grid points on the asset grid 'ag' to 201 (equal to the number of points n the asset grid 'a'). In this case, the computational time still amounts to 2 days : 1 hour : 11 minutes. The corresponding GAUSS program 'RCH7_denf.g' which uses exactly the same number of arithmetic operations only requires a computational time of 1 hour : 18 minutes. If you use C++, Fortran or Julia, the computational time can be reduced further.

There are various ways to speed up the computation in Python. As mention above, these include the use of the library Numpa or parallelization. In addtion, one could write single pieces of the code, for example the iteration over the dynamics of the distribution function, as C++ code and import it into Python.

Comments? Suggestions? Send an email to: Burkhard.Heer@wiwi.uni-augsburg.de

Last update: January 11, 2021