#+OPTIONS:   H:3 num:nil toc:nil \n:nil @:t ::t |:t ^:nil -:t f:t *:t TeX:t LaTeX:t skip:t d:t tags:not-in-toc creator:t timestamp:nil author:nil title:nil
#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="../../../css/stylesheet-blog.css">
#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="../../../fontawesome/css/all.css" />
#+HTML_HEAD: <link rel="stylesheet" href="../../../highlight/styles/base16/solarized-light.min.css">
#+HTML_HEAD: <script src="../../../highlight/highlight.min.js"></script>
#+HTML_HEAD: <script>hljs.highlightAll();</script>
#+LANGUAGE: en-GB
#+TITLE: Simulating the Solow model
#+DATE: February 2022
#+AUTHOR: Stéphane Adjemian
#+EMAIL: stephane.adjemian@univ-lemans.fr
#+PROPERTY: header-args:python :python /tmp/blog-solow-simulation/bin/python

#+BEGIN_QUOTE
This note shows how to simulate the Solow model with Python, discusses the
accuracy of the numerical simulations and compares them with what we would
obtain by simulating the linearised model.
#+END_QUOTE

* Numerical simulation of 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:

#+BEGIN_SRC bash :results silent :exports none :async t
python3 -m venv /tmp/blog-solow-simulation
source /tmp/blog-solow-simulation/bin/activate
pip install numpy scipy matplotlib
#+END_SRC


#+begin_src python :session solow-simulation-en  :exports code :results none
  def dotk(k, t, alpha, s, n, x, delta):
      return s*k**alpha-(n+x+delta)*k
#+end_src

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:

#+begin_src python :session solow-simulation-en :exports code :results none
  def kstar(alpha, s, n, x, delta):
      return (s/(n+x+delta))**(1/(1-alpha))
#+end_src

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:

#+begin_src python :session solow-simulation-en :exports code :results none
  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))
#+end_src

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.

#+begin_src python :session solow-simulation-en :exports none :results none
  import matplotlib.pyplot as plt
  plt.figure(1)
  plt.plot(t, kpath0, 'b')
  plt.plot(t, kpath1, 'r')
  plt.savefig("k-transition-1.svg", transparent=True)
#+end_src

#+CAPTION: *Paths of physical capital per efficient worker.*
#+LABEL: fig:k-transition-1
[[file:k-transition-1.svg]]

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:

\begin{equation*}
g(\hat k) = s \hat k^{\alpha-1} - (n+x+\delta)
\end{equation*}

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:

\begin{equation*}
g(\hat k) =  (n+x+\delta)\left(\frac{s}{n+x+\delta}\hat k^{\alpha-1} - 1\right)
\end{equation*}

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:

\begin{equation*}
g(\hat k) =  (n+x+\delta)\left(\left(\frac{\hat k}{\hat k^{\star}}\right)^{\alpha-1} - 1\right)
\end{equation*}

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)$).
#+begin_src python :session solow-simulation-en :exports none :results none
  plt.figure(2)
  plt.plot(kpath0, dotk(kpath0, t, alpha, s, n, x, delta)/kpath0, 'b')
  plt.plot(kpath1, dotk(kpath1, t, alpha, s, n, x, delta)/kpath1, 'b')
  plt.axvline(x=kstar(alpha, s, n, x, delta), c='r', linewidth=1, linestyle='--')
  plt.axhline(y=0, c='r', linewidth=1, linestyle='--')
  plt.savefig("k-growth-1.svg", transparent=True)
#+end_src

#+CAPTION: *Growth rate of the stock of physical capital per efficient worker.*
#+LABEL: fig:k-growth-1
[[file:k-growth-1.svg]]


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.

* Comparison with the exact solution

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)$.

#+begin_src python :session solow-simulation-en :exports code :results none
  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))
#+end_src

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.

#+begin_src python :session solow-simulation-en :exports none :results none
  vexact = np.vectorize(exact)
  # Initial conditions
  kinit0 = 0.0001*kstar(alpha, s, n, x, delta)
  kinit1 = 2.0000*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))
  aerror0 = vexact(t, kinit0, kstar(alpha, s, n, x, delta), alpha, (1-alpha)*(n+x+delta))-np.transpose(kpath0)
  aerror1 = vexact(t, kinit1, kstar(alpha, s, n, x, delta), alpha, (1-alpha)*(n+x+delta))-np.transpose(kpath1)
  plt.figure(3)
  plt.plot(kpath0, np.transpose(aerror0), 'b')
  plt.plot(kpath1, np.transpose(aerror1), 'b')
  plt.savefig("k-error-1.svg", transparent=True)
#+end_src

#+CAPTION: *Errors of the numerical solution*
#+LABEL: fig:accuracy-errors-1
[[file:k-error-1.svg]]


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.

* Comparison with the linear approximation

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 [[fig:k-growth-1]], the blue
curve by its tangent at $\hat k^{\star}$ (the green line in figure
[[fig:k-growth-2]]):

#+begin_src python :session solow-simulation-en :exports none :results none
  def tangent(k, alpha, s, n, x, delta):
      kinf = kstar(alpha, s, n, x, delta)
      slope = -(1-alpha)*s*kinf**(alpha-2)
      return slope*(k-kinf)
#+end_src

#+begin_src python :session solow-simulation-en :results none :exports none
  plt.figure(4)
  # Initial conditions
  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))
  plt.plot(kpath0, dotk(kpath0, t, alpha, s, n, x, delta)/kpath0, 'b')
  plt.plot(kpath1, dotk(kpath1, t, alpha, s, n, x, delta)/kpath1, 'b')
  plt.axvline(x=kstar(alpha, s, n, x, delta), c='r', linewidth=1, linestyle='--')
  plt.axhline(y=0, c='r', linewidth=1, linestyle='--')
  plt.plot(kpath0, tangent(kpath0, alpha, s, n, x, delta), 'g')
  plt.plot(kpath1, tangent(kpath1, alpha, s, n, x, delta), 'g')
  plt.savefig("k-growth-2.svg", transparent=True)
#+end_src

#+CAPTION: *Approximate growth rate of the stock of physical capital per efficient worker.*
#+LABEL: fig:k-growth-2
[[file:k-growth-2.svg]]

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
[[fig:k-approx-1]] 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.

#+begin_src python :session solow-simulation-en :exports none :results none
  def kapprox(k0, t, kinf, beta):
      from numpy import exp
      return kinf**(1-exp(-beta*t))*k0**exp(-beta*t)
#+end_src

#+begin_src python :session solow-simulation-en :exports none :results none
  plt.figure(5)
  # Initial conditions
  kinit0 = 0.001*kstar(alpha, s, n, x, delta)
  kinit1 = 2.000*kstar(alpha, s, n, x, delta)
  # Exact transitions
  kpath0 = vexact(t, kinit0, kstar(alpha, s, n, x, delta), alpha, (1-alpha)*(n+x+delta))
  kpath1 = vexact(t, kinit1, kstar(alpha, s, n, x, delta), alpha, (1-alpha)*(n+x+delta))
  # Approximate transitions
  apath0 = kapprox(kinit0, t, kstar(alpha, s, n, x, delta), (1-alpha)*(n+x+delta))
  apath1 = kapprox(kinit1, t, kstar(alpha, s, n, x, delta), (1-alpha)*(n+x+delta))
  plt.plot(t, kpath0, 'b')
  plt.plot(t, apath0, 'b', linestyle='--')
  plt.plot(t, kpath1, 'r')
  plt.plot(t, apath1, 'r', linestyle='--')
  plt.savefig("k-approx-1.svg", transparent=True)
#+end_src

#+CAPTION: *Approximation of the transition*
#+LABEL: fig:k-approx-1
[[file:k-approx-1.svg]]

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.

* Transition of the speed of convergence

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}$ [fn:1: The distance is measured here in
logarithms, and not in levels by $(\hat k(t)-\hat k^{\star})/\hat k^{\star}$, so
that $\beta(t)$ is directly comparable with $\beta^{\star}$: the previous
section indeed defines $\beta^{\star}$ as the rate of decay of the
/logarithmic/ distance to the steady state. The two measures coincide in a
neighbourhood of the steady state, since $\log \zeta \sim \zeta-1$ as $\zeta$
tends to 1, but they diverge, and even lead to opposite conclusions, far from
it.].

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[fn:3: We know that \[ \zeta(t) = \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 can compute
the speed of convergence at every instant. Figure [[fig:k-speed-1]] 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).

#+begin_src python :session solow-simulation-en :exports none :results none
  plt.figure(6)
  # Initial conditions
  kinit0 = 0.5*kstar(alpha, s, n, x, delta)
  kinit1 = 1.5*kstar(alpha, s, n, x, delta)
  # Exact transitions
  kpath0 = vexact(t, kinit0, kstar(alpha, s, n, x, delta), alpha, (1-alpha)*(n+x+delta))
  kpath1 = vexact(t, kinit1, kstar(alpha, s, n, x, delta), alpha, (1-alpha)*(n+x+delta))
  zeta0 = kpath0/kstar(alpha, s, n, x, delta)
  zeta1 = kpath1/kstar(alpha, s, n, x, delta)
  # Transitions of the speed of convergence
  speed0 = -(n+x+delta)*(zeta0**(alpha-1)-1)/np.log(zeta0)
  speed1 = -(n+x+delta)*(zeta1**(alpha-1)-1)/np.log(zeta1)
  plt.plot(kpath0, speed0, 'b')
  plt.plot(kpath1, speed1, 'r')
  plt.savefig("k-speed-1.svg", transparent=True)
  plt.figure(7)
  plt.plot(t, speed0, 'b')
  plt.plot(t, speed1, 'r')
  plt.savefig("k-speed-2.svg", transparent=True)
#+end_src

#+CAPTION: *Transition of the speed of convergence*
#+LABEL: fig:k-speed-1
[[./k-speed-1.svg]]

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 [[fig:k-speed-2]], $\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).

#+CAPTION: *Transition of the speed of convergence*
#+LABEL: fig:k-speed-2
[[./k-speed-2.svg]]

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.