Simulating the Solow model

We consider the standard Solow model with a Cobb-Douglas production function, a population and a labour efficiency growing at constant rates, and a constant saving rate. The law of motion of the aggregate stock of physical capital, \(K\), is therefore:

\[ \dot K(t) = s K(t)^{\alpha}\bigl(A(t)L(t)\bigr)^{1-\alpha} - \delta K(t) \]

where \(s\in]0,1[\) is the saving rate, \(\delta\in]0,1[\) the depreciation rate, \(L\) the population, \(A\) the efficiency of labour, and we denote by \(Y=K^{\alpha}\bigl(AL\bigr)^{1-\alpha}\) output (identical to income here). This equation tells us that the aggregate stock of physical capital of the economy increases if and only if investment exceeds the depreciation of capital.

For an arbitrary initial condition, we would like to determine the path of capital or of output from this equation. It is a nonlinear differential equation which we could solve, at least numerically, but as it stands the analysis of the dynamics would be awkward because this equation is not autonomous: the link between the change in \(K\) and its level is not stable over time, since the population and the efficiency of labour grow at every instant. We first transform the model by eliminating these two exogenous sources of growth. Denoting by \(n\) and \(x\) the (assumed constant) growth rates of the population \(L\) and of the efficiency of labour \(A\), one easily shows that the dynamics of the capital stock per efficient worker \(\hat k(t) = \frac{K(t)}{A(t)L(t)}\) is given by:

\[ \dot {\hat k}(t) = s \hat k(t)^{\alpha} - (n+x+\delta)\hat k(t) \]

The capital stock per efficient worker increases if and only if investment per efficient worker exceeds the depreciation of capital per efficient worker. We now have a time-invariant relation between the change and the level, and it will be simpler to study the dynamics from this equation. We can always go back to the aggregate variables, noting that by definition \(K(t)=\hat k(t)A(t)L(t)\).

This differential equation has a unique strictly positive steady state:

\[ \hat k^{\star} = \left(\frac{s}{n+x+\delta}\right)^{\frac{1}{1-\alpha}} \]

and since output per efficient worker is given by \(\hat y = \hat k^{\alpha}\), we also have:

\[ \hat y^{\star} = \left(\frac{s}{n+x+\delta}\right)^{\frac{\alpha}{1-\alpha}} \]

One can show that this steady state is globally asymptotically stable, that is, the capital stock per efficient worker converges to \(\hat k^{\star}\) for any initial condition \(\hat k_0>0\). One can also show that this convergence to the steady state is monotonic. We shall come back to the asymptotic behaviour of the dynamics below, by solving the differential equation analytically, but this result could be established more generally (this will be the subject of another note) without identifying explicitly the solution of the differential equation.

We use Python, together with the numpy and scipy libraries, to solve the differential equation numerically. We first have to define the problem to be solved. The following function returns the change in the capital stock per efficient worker associated with a level of capital per efficient worker:

def dotk(k, t, alpha, s, n, x, delta):
    return s*k**alpha-(n+x+delta)*k

Note that the second argument is useless in our case; it allows one to handle non-autonomous problems (our problem is autonomous since we have eliminated the two exogenous sources of growth). The following function returns the steady state:

def kstar(alpha, s, n, x, delta):
    return (s/(n+x+delta))**(1/(1-alpha))

we do not need to know the steady state to solve the differential equation, but it will allow us to choose an initial condition with a clear idea of the initial distance to the steady state (which in principle must also be the long-run level of the capital stock per efficient worker).

Finally, the following block of code uses the odeint function of the scipy library to solve the differential equation:

import numpy as np
from scipy.integrate import odeint
# Time discretisation
t = np.linspace(0, 200, 10000)
# Calibration of the model
alpha = .33
s = .20;
n = .02
x = .02
delta = .02
# Initial conditions (relative to the steady state)
kinit0 = 0.05*kstar(alpha, s, n, x, delta)
kinit1 = 1.95*kstar(alpha, s, n, x, delta)
# Solution of the differential equation (transition to the steady state)
kpath0 = odeint(dotk, kinit0, t, args=(alpha, s, n, x, delta))
kpath1 = odeint(dotk, kinit1, t, args=(alpha, s, n, x, delta))

The following figure shows two paths for the capital stock per efficient worker. The red path corresponds to the case where the initial condition is above the steady state. This path is monotonically decreasing, as expected. The blue path corresponds to the case where the initial condition is below the steady state; it is monotonically increasing. The two paths merge (the difference is no longer visible after 130), since they both converge to the same steady state.

Another way of representing the dynamics is to illustrate the relation between the growth rate (on the vertical axis) and the level (on the horizontal axis). The growth rate is defined as the ratio of the change in \(\hat k\) to its level. We therefore have:

Clearly, because of the assumption of a decreasing marginal return to capital (\(\alpha < 1\)), the growth rate is a monotonically decreasing function of the level. The growth rate can be rewritten more elegantly by bringing out the steady state. Factoring out the depreciation rate of capital per efficient worker:

then expressing the ratio of the saving rate to the depreciation rate of capital per efficient worker in terms of the steady state, we finally obtain the following expression of the growth rate:

One easily checks that the growth rate is positive if and only if \(\hat k < \hat k^{\star}\) and that the growth rate is a monotonically decreasing function of the level of the capital stock per efficient worker (both properties are consequences of the assumption \(\alpha < 1\)).

The following figure shows the function \(g(\hat k)\) (the blue curve); it does go through 0 at the steady state (intersection of the red lines) and, in line with what we saw above, it is monotonically decreasing (one can guess the vertical asymptote at 0 and the horizontal asymptote at \(-(n+x+\delta)\)).

All in all, these numerical simulations are in line with what we expected. Simulating this model is not necessarily very interesting, but the numerical approach can prove useful, even indispensable, once the model is made more complex.

One should nevertheless question the reliability of this numerical solution. From a qualitative point of view we have just seen that there was a priori no problem here (we recover the expected properties). In what follows we address this issue from a quantitative point of view. Without going into a detailed description of the algorithm behind the odeint function, it is important to keep in mind that a numerical approach always relies on approximations. One must therefore ask whether the approximation errors are large or negligible. It is also useful to know where these errors are possibly larger.

Fortunately, this model has an analytical solution. In the next section we present this solution and compare the exact solution with a solution obtained by numerical simulation.

We know that the dynamics of the capital stock per efficient worker obeys the following differential equation:

\[ \dot {\hat k}(t) = s \hat k(t)^{\alpha} - (n+x+\delta)\hat k(t) \]

This nonlinear equation cannot be solved directly, but a simple change of variable reduces it to a linear differential equation (which we know how to solve). Rather than studying the dynamics of the capital stock per efficient worker, we shall look at the dynamics of the inverse of the average product of capital, \(z(t)=\hat k(t)^{1-\alpha}\). Excluding the case \(\hat k(t)=0\), the differential equation can be rewritten as:

\[ \dot {\hat k}(t) {\hat k}(t)^{-\alpha} = s - (n+x+\delta)\hat k(t)^{1-\alpha} \]

One recognises \(z(t)\) on the right-hand side. Noting that the derivative of the inverse of the average product of capital with respect to time is \(\dot z(t) = (1-\alpha) {\hat k}(t)^{-\alpha} \dot {\hat k}(t)\), the differential equation can finally be written as:

\[ \dot z(t) = (1-\alpha)s - (1-\alpha)(n+x+\delta) z(t) \]

The dynamics of the inverse of the average product of capital is linear! We know how to solve this differential equation (a linear ODE with a constant right-hand side and constant coefficients). The general solution of the complete equation is the sum of a particular solution (the steady state \(z^{\star}=\frac{s}{n+x+\delta}\)) and of the general solution of the homogeneous equation (\(Ae^{-(1-\alpha)(n+x+\delta)t}\)):

\[ z(t) = \frac{s}{n+x+\delta} + A e^{-(1-\alpha)(n+x+\delta)t} \]

Given an initial condition \(z(0)\), the solution satisfying this condition is:

\[ z(t) = \frac{s}{n+x+\delta} + \left(z(0)-\frac{s}{n+x+\delta}\right) e^{-(1-\alpha)(n+x+\delta)t} \]

which we can rewrite in terms of \(\hat k\):

\[ \hat k(t) = \left(\frac{s}{n+x+\delta} + \left(\hat k(0)^{1-\alpha}-\frac{s}{n+x+\delta}\right) e^{-(1-\alpha)(n+x+\delta)t}\right)^{\frac{1}{1-\alpha}} \]

or equivalently:

\[ \hat k(t) = \hat k^{\star} \left(1 + \left(\left(\frac{\hat k(0)}{\hat k^{\star}}\right)^{1-\alpha}-1\right) e^{-(1-\alpha)(n+x+\delta)t}\right)^{\frac{1}{1-\alpha}} \]

One easily checks that \(\lim_{t\rightarrow\infty} \hat k(t) = \hat k^{\star}\) for any positive initial condition. One checks just as easily that the path of capital per efficient worker is monotonically increasing if \(\hat k(0) < \hat k^{\star}\) and monotonically decreasing if \(\hat k(0)>\hat k^{\star}\). We thus have that, in this model, the steady state is globally asymptotically stable. In passing, one recognises under the exponential the speed of convergence \(\beta^{\star} = (1-\alpha)(n+x+\delta)\).

We now wish to compare this solution, exact by construction, with the numerical solution. The following function computes \(\hat k\) at any instant \(t\) given an initial condition \(\hat k(0)\).

def exact(t, k0, kinf, alpha, beta):
    from numpy import exp
    return kinf*(1+((k0/kinf)**(1-alpha)-1)*exp(-beta*t))**(1/(1-alpha))

The following chart shows the numerical error, that is the difference between the exact solution and the numerical solution, against the level of the capital stock per efficient worker. Clearly the numerical error is very small (of the order of \(10^{-7}\) for a variable taking values between 0 and 12). We could probably make it even smaller by changing the defaults of the odeint routine.

The conclusion is therefore that the odeint routine does a very good job here. One notes, though this may not be a general result, that the errors are smaller (in absolute value) close to the steady state. In the lectures, we obtained the quantitative predictions of the Solow model (for instance the speed of convergence) by resorting to a (log-)linear approximation of the Solow model. In the next section we use the exact solution to compute the approximation errors in that case.

The linear approximation is built from a first-order Taylor expansion of the function \(g(\hat k)\) in a neighbourhood of the steady state \(\hat k^{\star}\):

\[ g(\hat k) \approx g(\hat k^{\star}) + g'(\hat k^{\star})\left(\hat k - \hat k^{\star}\right) \]

or, since the growth rate is zero at the steady state:

\[ g(\hat k) \approx g'(\hat k^{\star})\left(\hat k - \hat k^{\star}\right) \]

Graphically, this amounts to replacing, in figure 2, the blue curve by its tangent at \(\hat k^{\star}\) (the green line in figure 4):

One can see on this chart that the approximation errors (the difference between the blue curve and the green line) are here much larger than the numerical errors computed in the previous section. Note that the errors:

• are larger the further away one is from the steady state,
• are asymmetric, in the sense that they are larger on the left than on the right, which is explained by the greater curvature of the function \(g(\hat k)\) on the left,
• always have the same sign: since \(g(\hat k)\) is a convex function, the approximate growth rate is always below the theoretical growth rate.

These differences, which we observe here graphically, obviously have consequences for the paths we can compute. With the first-order approximation, we have:

\[ g(\hat k) \approx -(1-\alpha)s \left. \hat k^{\star}\right. ^{\alpha-1}\frac{\hat k - \hat k^{\star}}{\hat k^{\star}} \]

Substituting the definition of the steady state, we get:

\[ g(\hat k) \approx -(1-\alpha)(n+x+\delta)\frac{\hat k - \hat k^{\star}}{\hat k^{\star}} \]

Recalling that the growth rate of \(\hat k\) can be written as the change in \(\log \hat k\), we have equivalently:

\[ \dot{\log \hat k} \approx (1-\alpha)(n+x+\delta)\left(1-\frac{\hat k}{\hat k^{\star}}\right) \]

Finally, the term \(1-\frac{\hat k}{\hat k^{\star}}\), which measures the distance to the steady state, can be approximated by a \(\log\) function in a neighbourhood of the steady state (we know that \(\log(1-x)\approx -x\)):

\[ \dot{\log \hat k} \approx -(1-\alpha)(n+x+\delta)\log \left(\frac{\hat k}{\hat k^{\star}}\right) \]

or equivalently:

\[ \dot{\log \frac{\hat k}{\hat k^{\star}}} \approx -(1-\alpha)(n+x+\delta)\log \frac{\hat k}{\hat k^{\star}} \]

This equation tells us that the growth rate of the distance to the steady state, measured by \(\log \frac{\hat k}{\hat k^{\star}}\), is constant and equal to \(-(1-\alpha)(n+x+\delta)\). We denote by \(\beta^{\star}=(1-\alpha)(n+x+\delta)\) the speed of convergence; it is the rate at which the distance to the steady state decays. This is a homogeneous linear differential equation with constant coefficients, which we know how to solve. Hence, the distance to the steady state at any instant \(t\), given an initial distance to the steady state, is given by:

\[ \log\frac{\hat k(t)}{\hat k^{\star}} \approx \log\frac{\hat k(0)}{\hat k^{\star}} e^{-\beta^{\star} t} \]

\[ \Leftrightarrow \log \hat k(t) \approx \left(1-e^{-\beta^{\star} t}\right)\log \hat k^{\star} + e^{-\beta^{\star} t}\log\hat k(0) \]

\[ \Leftrightarrow \hat k(t) \approx {\hat k^{\star}}^{1-e^{-\beta^{\star} t}} {\hat k(0)}^{e^{-\beta^{\star} t}} \]

The last equation, derived from the log-linearisation of the dynamics, must be compared with the exact equation obtained above for \(\hat k(t)\). Figure 5 compares the approximate (dashed) and exact (solid) dynamics. As we could anticipate, the difference is larger when the initial condition is below the steady state, because it is when the capital stock approaches zero that the curvature of the model becomes larger.

Note that the approximate transition, when the initial condition is close to zero, does not have the same curvature as the exact transition. It also seems significantly slower to reach the steady state. This suggests that the measure of the speed of adjustment we are used to (\(\beta^{\star} = (1-\alpha)(n+x+\delta)\), obtained by log-linearising the dynamics) underestimates the speed of adjustment towards the steady state.

The speed of adjustment towards the steady state can be defined, more generally, as the rate of decay of the distance to the steady state. Without resorting to approximations, as in the previous section, we set:

\[ \beta(t) = - \frac{\frac{\mathrm d}{\mathrm dt}\log\frac{\hat k(t)}{\hat k^{\star}}}{\log\frac{\hat k(t)}{\hat k^{\star}}} \]

where the distance to the steady state is measured by \(\log \hat k(t)/\hat k^{\star}\) ¹.

Substituting the law of motion of the stock of physical capital, one easily obtains the following expression for the speed of convergence:

\[ \beta (t) = -(n+x+\delta) \frac{\zeta(t)^{\alpha-1}-1}{\log \zeta(t)} \]

with \(\zeta(t) = \hat k(t) / \hat k^{\star}\). We know that along the transition \(\zeta(t)\) approaches 1 monotonically [fn:2: If the economy starts above the steady state, \(\zeta(t)\geq 1\) for all \(t\) and decreases monotonically towards

1. If the economy starts below the steady state, \(\zeta(t)\leq 1\) for all \(t\)

and increases monotonically towards 1.], so the speed of adjustment will vary during the transition. As \(t\) tends to infinity (or \(\zeta\) tends to 1), the numerator and the denominator of the ratio in the last equation both tend to 0; we thus have an indeterminate form, but using l'Hôpital's rule we can show that this ratio tends to \(\alpha-1\) as \(\zeta\) tends to 1, and hence that:

\[ \lim_{t\rightarrow\infty} \beta(t) = (1-\alpha)(n+x+\delta)\equiv \beta^{\star} \]

the (constant) speed of convergence we obtain when considering the log-linear approximation in a neighbourhood of the steady state (see the previous section). Substituting the exact solution into the expression of the (exact) speed of convergence² one can compute the speed of convergence at every instant. Figure 6 shows the speed of convergence as a function of the level of capital per efficient worker. To make the chart more readable, the deviations from the steady state have been reduced (we consider values of \(\hat k\) between \(50\%\) and \(150\%\) of the steady state).

One observes that the exact speed of convergence departs further from the constant speed of convergence \(\beta^{\star}\), obtained by approximating the model in a neighbourhood of the steady state (\(4\%\) for our calibration of the model), when the economy is below the steady state (the blue curve). In figure 7, \(\beta(t)\) is plotted against time, with \(\hat k(0) = .5\hat k^{\star}\) (blue curve) or \(\hat k(0) = 1.5\hat k^{\star}\) (red curve).

When the initial endowment of the economy is below the steady state, one observes that the exact speed of convergence is above \(\beta^{\star}\) for all \(t\) and decreases monotonically. The deviations from \(\beta^{\star}\) are larger when the economy is below the steady state. Indeed, when the economy is above the steady state the speed of convergence is bounded below by 0 (it must be positive), whereas it can become arbitrarily large as \(\zeta\) tends to 0. But above all, we know that the curvature of the model is larger when the capital stock approaches 0, which explains why the difference between the exact speed of convergence and the speed of convergence resulting from a local approximation is larger. In the limit, as \(\hat k\) (or \(\zeta\)) tends to 0, the Inada condition tells us that the marginal return to capital is infinite. In that case, the speed of convergence tends to infinity.