#+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: fr
#+TITLE: Simulation du modèle de Solow
#+DATE: Février 2022
#+AUTHOR: Stéphane Adjemian
#+EMAIL: stephane.adjemian@univ-lemans.fr

#+BEGIN_QUOTE
Cette note montre comment simuler le modèle de Solow, à l'aide de Python,
discute l'exactitude des simulations numériques et compare celles-ci avec ce que
nous obtiendrions en simulant le modèle linéarisé.
#+END_QUOTE

* Simulation numérique du modèle de Solow

On considère le modèle de Solow standard avec une fonction de production
Cobb-Douglas, une population et une efficacité du travail qui croissent à taux
constant et un taux d'épargne constant. La loi d'évolution du stock de capital
physique agrégé, $K$, est donc :

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

où $s\in]0,1[$ est le taux d'épargne, $\delta\in]0,1[$ le taux de dépréciation,
$L$ la population, $A$ l'efficacité du travail, et on notera
$Y=K^{\alpha}\bigl(AL\bigr)^{1-\alpha}$ la production (identique au revenu dans
ici). Cette équation nous dit que le stock de capital physique agrégé de
l'économie augmente si et seulement si l'investissement est supérieur à la
dépréciation du capital.

\\

On aimerait, pour une condition initiale arbitraire, déterminer la trajectoire
du capital ou de la production à partir de cette équation. Il s'agit d'une
équation différentielle non linéaire que nous pourrions résoudre, au moins
numériquement, mais en l'état l'analyse de la dynamique serait malaisée car
cette équation n'est pas autonome : le lien entre la variation de $K$ et son
niveau n'est pas stable dans le temps, car la population et
l'efficacité du travail augmentent à chaque instant. Au préalable on
transforme le modèle en éliminant ces deux sources exogènes de croissance. Si on note $n$
et $x$ les taux de croissance (supposés constants) de la population $L$ et de
l'efficacité du travail $A$, on peut facilement montrer que la dynamique du
stock de capital par tête efficace $\hat k(t) = \frac{K(t)}{A(t)L(t)}$ est
donnée par :

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

Le stock de capital par tête efficace augmente si et seulement si
l'investissement par tête efficace est supérieur à la dépréciation du capital
par tête efficace. Nous avons bien ici une relation entre la variation et
le niveau invariante dans le temps, il sera plus simple d'étudier la dynamique à
partir de cette équation. On pourra toujours revenir aux variables agrégées, en
notant que par définition on a $K(t)=\hat k(t)A(t)L(t)$.

\\

Cette équation différentielle admet un unique état stationnaire strictement
positif :

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

et puisque la production par tête est donnée par $\hat y = \hat k^{\alpha}$, on a aussi :

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

On peut montrer que cet état stationnaire est asymptotiquement globalement
stable, c'est-à-dire que le stock de capital par tête efficace converge vers
$\hat k^{\star}$ pour toute condition initiale $\hat k_0>0$. On peut aussi montrer
que cette convergence vers l'état stationnaire est monotone. Nous reviendrons
sur le comportement asymptotique de la dynamique plus loin, en résolvant
analytiquement l'équation différentielle, mais nous pourrions établir de façon
plus générale ce résultat (cela fera l'objet d'un autre article) sans identifier
explicitement la solution de l'équation différentielle.

\\

Nous utilisons Python, ainsi que les librairies =numpy= et =scipy= pour résoudre
numériquement l'équation différentielle. Nous devons d'abord définir le problème
à résoudre. La fonction suivante retourne la variation du stock de capital par
tête efficace associée à un niveau de capital par tête efficace :

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

#+RESULTS:

Notons que le second argument est inutile dans notre cas, il permet de traiter
les problèmes non autonome (notre problème est autonome puisque nous avons
éliminés les deux sources exogènes de croissance). La fonction suivante retourne l'état stationnaire :

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

#+RESULTS:

nous n'avons pas besoin de connaître l'état stationnaire pour résoudre
l'équation différentielle, mais cela nous permettra de choisir une condition
initiale en ayant les idées précises sur la distance initiale à l'état
stationnaire (qui en principe doit aussi être le niveau de long terme du stock
de capital par tête efficace).

\\

Finalement le bloc de code suivant utilise la fonction =odeint= de la librairie
=scipy= pour résoudre l'équation différentielle :

#+begin_src python :session :exports code
  import numpy as np
  from scipy.integrate import odeint
  # Discrétisation du temps
  t = np.linspace(0, 200, 10000)
  # Étalonnage du modèle
  alpha = .33
  s = .20;
  n = .02
  x = .02
  delta = .02
  # Conditions initiales (en déviation à l'état stationnaire)
  kinit0 = 0.05*kstar(alpha, s, n, x, delta)
  kinit1 = 1.95*kstar(alpha, s, n, x, delta)
  # Résolution de l'équation différentielle (transition vers l'état stationnaire)
  kpath0 = odeint(dotk, kinit0, t, args=(alpha, s, n, x, delta))
  kpath1 = odeint(dotk, kinit1, t, args=(alpha, s, n, x, delta))
#+end_src

#+RESULTS:

La figure suivante représente deux trajectoires pour le stock de capital par
tête efficace. La trajectoire rouge correspond au cas où la condition initiale
est au dessus de l'état stationnaire. Cette trajectoire est monotone
décroissante, comme attendu. La trajectoire bleue correspond au cas où la
condition initiale est inférieure à l'état stationnaire, elle est monotone
croissante. Les deux trajectoires se rejoignent (on ne distingue plus la
différence après 130), car elles convergent toutes les deux vers le même état
stationnaire.

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

#+RESULTS:
: None

#+CAPTION: *Trajectoires du capital physique par tête efficace.*
#+LABEL: fig:k-transition-1
[[file:img/k-transition-1.svg]]

Une autre façon de représenter la dynamique est d'illustrer la relation entre le
taux de croissance (sur l'axe des ordonnées) et le niveau (sur l'axe des
abscisses). Le taux de croissance est défini comme le rapport de la variation de
$\hat k$ et de son niveau. On a donc :

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

Clairement, à cause de l'hypothèse de rendement marginal décroissant du capital
($\alpha<1$), le taux de croissance est une fonction monotone décroissante du
niveau. On peut réécrire le taux de croissance de façon plus élégante en faisant
apparaître l'état stationnaire. En factorisant le taux de dépréciation du
capital par tête efficace :

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

puis en exprimant le rapport du taux d'épargne et du taux de dépréciation du capital
par tête efficace en fonction de l'état stationnaire, on obtient finalement
l'expression suivante du taux de croissance :

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

On vérifie facilement que le taux de croissance est positif si et seulement si
$\hat k<\hat k^{\star}$ et que le taux de croissance est une fonction monotone
décroissante du niveau de stock de capital par tête efficace (ces deux
propriétés sont des conséquences de l'hypothèse $\alpha<1$).

\\

La figure suivante représente la fonction $g(\hat k)$ (la courbe bleue), elle
passe bien par 0 à l'état stationnaire (intersection des droites rouges) et,
conformément à ce que nous pouvions voir plus haut, elle est monotone
décroissante (on devine l'asymptote verticale en 0 et l'asymptote horizonale en
$-(n+x+\delta)$).
#+begin_src python :session :exports 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("img/k-growth-1.svg", transparent=True)
#+end_src

#+RESULTS:
: None

#+CAPTION: *Taux de croissance du stock de capital physique par tête efficace.*
#+LABEL: fig:k-growth-1
[[file:img/k-growth-1.svg]]


Au total ces simulations numériques sont conformes à ce que nous attendions. Il
n'est pas forcément très intéressant de simuler ce modèle, mais l'approche
numérique peut se révéler intéressante, voire indispensable, si on complexifie
le modèle.

\\

Il convient néanmoins de s'interroger sur la pertinence de cette solution
numérique. D'un point de vue qualitatif nous venons de voir qu'il n'y avait a
priori pas de problème ici (nous retrouvons les propriétés attendues). Dans la
suite nous allons aborder ce problème d'un point de vue quantitatif. Sans
rentrer dans une description détaillée de l'algorithme derrière la fonction
=odeint,= il est important de garder à l'esprit qu'une approche numérique repose
toujours sur des approximations. Dès lors, il faut se demander si les erreurs
d'approximation sont importantes ou négligeables. Il est aussi utile de savoir
où ces erreurs sont éventuellement plus importantes.

Par chance, ce modèle dispose d'une solution analytique. Dans la prochaine
section nous présentons cette solution et comparons la solution exacte avec une
solution obtenue à l'aide de simulations numériques.

* Comparaison avec la solution exacte

Nous savons que la dynamique du stock de capital par tête obéit à l'équation différentielle suivante :

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

Cette équation non linéaire ne peut être résolue directement, mais à l'aide d'un
changement de variable simple on peut se ramener à une équation différentielle
linéaire (que nous savons résoudre). Plutôt que d'étudier la dynamique du stock
de capital par tête efficace, nous allons nous intéresser à la dynamique de
l'inverse de la productivité moyenne du capital, $z(t)=\hat k(t)^{1-\alpha}$. En
excluant le cas où $\hat k(t)=0$, on peut réécrire l'équation différentielle
sous la forme :

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

On reconnaît $z(t)$ sur le membre de droite. En remarquant que la dérivée de
l'inverse de la productivité moyenne du capital par rapport au temps est $\dot
z(t) = (1-\alpha) {\hat k}(t)^{-\alpha} \dot {\hat k}(t)$, on peut finalement
écrire l'équation différentielle sous la forme :

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

La dynamique de l'inverse de la productivité moyenne du capital est linéaire !
Nous savons résoudre cette équation différentielle (une EDL avec second membre à
coefficients constants). La solution générale de l'équation complète est la
somme d'une solution particulière (l'état stationnaire
$z^{\star}=\frac{s}{n+x+\delta}$) et de la solution générale de l'équation sans
second membre ($Ae^{-(1-\alpha)(n+x+\delta)t}$) :

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

Si on dispose d'une condition initiale $z(0)$, alors la solution qui satisfait cette condition est :

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

que nous pouvons réécrire avec $\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}}
\]

ou encore :

\[
\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}}
\]

On vérifie facilement que $\lim_{t\rightarrow\infty} \hat k(t) = \hat k^{\star}$
pour toute condition initiale positive. On vérifie tout aussi facilement que la
trajectoire de capital par tête efficace est monotone croissante si $\hat
k(0)<\hat k^{\star}$ et monotone décroissante si $\hat k(0)>\hat k^{\star}$.
Nous avons donc bien que, dans ce modèle, l'état stationnaire est
asymptotiquement globalement stable. En passant, on reconnaît sous
l'exponentielle la vitesse de convergence $\beta^{\star} = (1-\alpha)(n+x+\delta)$.

\\

Nous souhaitons maintenant comparer cette solution, par construction exacte,
avec la solution numérique. La fonction suivante permet de calculer $\hat k$ a
un instant $t$ quelconque étant donnée une condition initiale $\hat k(0)$.

#+begin_src python :session :exports code
  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

#+RESULTS:

Le graphique suivant montre l'erreur numérique, c'est-à-dire la différence entre
la solution exacte et la solution numérique, contre le niveau du stock de
capital par tête efficace. Clairement l'erreur numérique est très petite (de
l'ordre de $10^{-7}$ pour une variable prenant des valeurs entre 0 et 12). Nous
pourrions probablement la rendre encore plus petite en changeant les défauts de
la routine =odeint=.

#+begin_src python :session :exports none
  vexact = np.vectorize(exact)
  # Condition initiale (50% de l'état stationnaire)
  kinit0 = 0.0001*kstar(alpha, s, n, x, delta)
  kinit1 = 2.0000*kstar(alpha, s, n, x, delta)
  # Résolution de l'équation différentielle (transition vers l'état stationnaire)
  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("img/k-error-1.svg", transparent=True)
#+end_src

#+RESULTS:
: None

#+CAPTION: *Erreurs de la solution numérique*
#+LABEL: fig:accuracy-errors-1
[[file:img/k-error-1.svg]]


La conclusion est donc que la routine =odeint= fait ici un très bon travail. On
remarque, ce n'est peut-être pas un résultat général, que les erreurs sont plus
petites (en valeur absolue) quand on est proche de l'état stationnaire. En
cours, nous avons obtenu les prédictions quantitatives du modèle de Solow (par
exemple la vitesse de convergence) en recourant à une approximation (log)
linéaire du modèle de Solow. Dans la prochaine section nous utilisons la
solution exacte pour calculer les erreurs d'approximation dans ce cas.

* Comparaison avec l'approximation linéaire

L'approximation linéaire est construite en faisant une approximation de Taylor
d'ordre 1 de la fonction $g(\hat k)$ dans un voisinage de l'état stationnaire
$\hat k^{\star}$:

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

ou, puisque le taux de croissance est nul à l'état stationnaire, encore :

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

Graphiquement, cela revient à remplacer dans la figure [[fig:k-growth-1]] la
courbe bleue par sa tangente en $\hat k^{\star}$ (la droite verte dans la figure [[fig:k-growth-2]]) :

#+begin_src python :session :exports none
  def tangente(k, alpha, s, n, x, delta):
      kinf = kstar(alpha, s, n, x, delta)
      pente = -(1-alpha)*s*kinf**(alpha-2)
      return pente*(k-kinf)
#+end_src

#+RESULTS:

#+begin_src python :session :results file :exports none
  plt.figure(4)
  # Condition initiale (50% de l'état stationnaire)
  kinit0 = 0.05*kstar(alpha, s, n, x, delta)
  kinit1 = 1.95*kstar(alpha, s, n, x, delta)
  # Résolution de l'équation différentielle (transition vers l'état stationnaire)
  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, tangente(kpath0, alpha, s, n, x, delta), 'g')
  plt.plot(kpath1, tangente(kpath1, alpha, s, n, x, delta), 'g')
  plt.savefig("img/k-growth-2.svg", transparent=True)
#+end_src

#+RESULTS:
[[file:None]]

#+CAPTION: *Taux de croissance approximé du stock de capital physique par tête efficace.*
#+LABEL: fig:k-growth-2
[[file:img/k-growth-2.svg]]

On devine sur ce graphique que les erreurs d'approximation (la différence entre
la courbe bleue et la droite verte) sont ici beaucoup plus importantes que les
erreurs numériques calculées dans la section précédente. On remarque que les erreurs :

 - sont d'autant plus importantes que l'on s'éloigne de l'état stationnaire,

 - sont asymétriques, au sens où elles sont plus importantes à gauche qu'à droite ce qui
   s'explique par la courbure plus importante à gauche de la fonction $g(\hat
   k)$,

 - ont toujours le même signe, puisque $g(\hat k)$ est une fonction
   convexe, le taux de croissance approximé est toujours inférieur au taux de
   croissance théorique.

\\

Ces différences, que nous observons ici graphiquement, ont évidemment des
conséquences sur les trajectoires que nous pouvons calculer. Avec
l'approximation d'ordre 1, nous avons :

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

En substituant lé définition de l'état stationnaire, il vient :

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

En rappelant que le taux de croissance de $\hat k$ peut s'écrire comme la variation de $\log \hat k$, on a de façon équivalente :

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

Finalement, le terme $1-\frac{\hat k}{\hat k^{\star}}$, qui mesure la distance
à l'état stationnaire, peut être approximé par une fonction $\log$ dans un voisinage
de l'état stationnaire (on sait que $\log(1-x)\approx -x$) :


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

ou encore :


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

Cette équation nous dit que le taux de croissance de la distance à l'état
stationnaire, mesurée par $\log \frac{\hat k}{\hat k^{\star}}$, est constant et
égal à $-(1-\alpha)(n+x+\delta)$. On note $\beta^{\star}=(1-\alpha)(n+x+\delta)$ la
vitesse de convergence, c'est le taux de décroissance de la distance à l'état
stationnaire. Il s'agit d'une équation différentielle linéaire à coefficients
constants et sans second membre, que nous savons résoudre. Ainsi, la distance à
l'état stationnaire à un instant $t$ quelconque, étant donnée une distance
initiale à l'état stationnaire est donnée par :

\[
\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}}
\]

La dernière équation, issue de la log-linéarisation de la dynamique, doit être
comparée à l'équation exacte obtenue plus haut pour $\hat k(t)$. La figure
[[fig:k-approx-1]] compare les dynamiques approximées (tirets) et exactes (trait
plein). Comme nous pouvions l'anticiper la différence est plus importante quand
on considère une condition initiale inférieure à l'état stationnaire, car c'est
quand le stock de capital se rapproche de zéro que la courbure du modèle devient
plus importante.

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

#+RESULTS:

#+begin_src python :session :exports none
  plt.figure(5)
  # Conditions initiales
  kinit0 = 0.001*kstar(alpha, s, n, x, delta)
  kinit1 = 2.000*kstar(alpha, s, n, x, delta)
  # Transitions exactes
  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))
  # Transitions approximées
  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("img/k-approx-1.svg", transparent=True)
#+end_src

#+RESULTS:
: None

#+CAPTION: *Approximation de la transition*
#+LABEL: fig:k-approx-1
[[file:img/k-approx-1.svg]]

On remarque que la transition approximée, dans le cas où la condition initiale
est proche de zéro, n'a pas la même courbure que la transition exacte. Elle
semble aussi significativement plus lente à rejoindre l'état stationnaire. Ceci
suggère que la mesure de la vitesse d'ajustement que nous avons l'habitude
d'utiliser ($\beta^{\star} = (1-\alpha)(n+x+\delta)$, obtenue en log-linéarisant la
dynamique), sous estime la vitesse d'ajustement vers l'état stationnaire.

* Transition de la vitesse de convergence

La vitesse d'ajustement vers l'état stationnaire peut être définie, de façon plus
générale, comme le taux de décroissance de la distance à l'état stationnaire.
Sans recourir à des approximations, comme dans la section précédente, on posera :

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

En substituant la loi d'évolution du stock de capital physique, on obtient
facilement l'expression suivante pour la vitesse de convergence :

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

avec $\zeta(t) = \hat k(t) / \hat k^{\star}$. On sait que le long de la
transition $\zeta(t)$ va se rapprocher de 1 de façon monotone [fn:1: Si
l'économie est initialement au dessus de l'état stationnaire, $\zeta(t)\geq 1$
pour tout $t$ et décroît de façon monotone pour converger vers 1. Si l'économie
est initialement sous l'état stationnaire, $\zeta(t)\leq 1$ pour tout $t$ et
croît de façon monotone pour converger vers 1.], la vitesse d'ajustement va donc
varier pendant la transition. Lorsque $t$ tend vers l'infini (ou $\zeta$ tend vers 1), le numérateur et
le dénominateur du ratio dans la dernière équation tendent vers 0 ; nous avons
donc ici une forme indéterminée, mais à l'aide de la règle de l'Hôpital nous pouvons
montrer que ce ratio tend vers $\alpha-1$ lorsque $\zeta$ tend vers 1 et donc
que :

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

la vitesse de convergence (constante) que nous obtenons lorsque nous considérons
l'approximation log-linéaire dans un voisinage de l'état stationnaire (voir la
section précédente). En substituant la solution exacte dans l'expression de la
vitesse de convergence (exacte)[fn:2: On sait que \[ \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}} \]] on peut calculer la
vitesse de convergence à chaque instant. La figure [[fig:k-speed-1]] représente la
vitesse de convergence comme une fonction du niveau de capital par tête
efficace. Afin que le graphique soit plus lisible on a réduit les écarts à
l'état stationnaire (on considère des valeurs de $\hat k$ entre $50\%$ et $150\%$ de
l'état stationnaire).

#+begin_src python :session :exports none
  plt.figure(6)
  # Conditions initiales
  kinit0 = 0.5*kstar(alpha, s, n, x, delta)
  kinit1 = 1.5*kstar(alpha, s, n, x, delta)
  # Transitions exactes
  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 de la vitesse de convergence
  speed0 = -(n+x+delta)*(zeta0**(alpha-1)-1)/(zeta0-1)
  speed1 = -(n+x+delta)*(zeta1**(alpha-1)-1)/(zeta1-1)
  plt.plot(kpath0, speed0, 'b')
  plt.plot(kpath1, speed1, 'r')
  plt.savefig("img/k-speed-1.svg", transparent=True)
  plt.figure(7)
  plt.plot(t, speed0, 'b')
  plt.plot(t, speed1, 'r')
  plt.savefig("img/k-speed-2.svg", transparent=True)
#+end_src

#+RESULTS:
: None

#+CAPTION: *Transition de la vitesse de convergence*
#+LABEL: fig:k-speed-1
[[./img/k-speed-1.svg]]

On observe que la vitesse de convergence exacte s'écarte plus de la vitesse de
convergence constante $\beta^{\star}$, obtenue en approximant le modèle dans un
voisinage de l'état stationnaire ($4\%$ pour notre étalonnage du modèle), lorsque
l'économie est sous l'état stationnaire (la courbe bleue). Dans la figure
[[fig:k-speed-2]], on représente $\beta(t)$ contre le temps, avec $\hat k(0) =
.5\hat k^{\star}$ (courbe bleue) ou $\hat k(0) = 1.5\hat k^{\star}$ (courbe
rouge).

#+CAPTION: *Transition de la vitesse de convergence*
#+LABEL: fig:k-speed-2
[[./img/k-speed-2.svg]]

Quand la dotation initiale de l'économie est inférieure à l'état stationnaire,
on observe que la vitesse de convergence exacte est supérieure à $\beta^{\star}$ pour
tout $t$ et décroît de façon monotone. Les écarts à $\beta^{\star}$ sont plus importants
quand l'économie est sous l'état stationnaire. En effet, quand l'économie est au
dessus de l'état stationnaire la vitesse de convergence est bornée par 0 (elle
doit être positive) alors qu'elle peut devenir arbitrairement grande quand
$\zeta$ tend vers 0. Mais surtout, nous savons que la courbure du modèle est
plus importante quand le stock de capital se rapproche de 0, ce qui explique que
la différence entre la vitesse de convergence exacte et la vitesse de
convergence résultant d'une approximation locale soit plus importante. À la
limite, quand $\hat k$ (ou $\zeta$) tend vers 0, la condition d'Inada nous dit
que le rendement marginal du capital est infini. Dans ce cas, la vitesse de
convergence tend vers l'infini.
