22 KiB
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é.
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) \]
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, à cause de la population et de l'efficacité du travail qui 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 :
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.
\[ \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 :
def dotk(k, t, alpha, s, n, x, delta):
return s*k**alpha-(n+x+delta)*k def kstar(alpha, s, n, x, delta):
return (s/(n+x+delta))**(1/(1-alpha))
Finalement le bloc de code suivant utilise la fonction odeint de la librairie
scipy pour résoudre l'équation différentielle :
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
# 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))
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α-1 - (n+x+δ) \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+δ)≤ft(\frac{s}{n+x+δ}\hat kα-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+δ)≤ft(≤ft(\frac{\hat k}{\hat k⋆}\right)α-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)$).
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 = (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)$.
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))
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) :
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édante. 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 \hat k^{\star}^{\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 mersure la distance à l'état stationnaire, peut 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=(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 t} \]
\[ \Leftrightarrow \log \hat k(t) \approx \left(1-e^{-\beta t}\right)\log \hat k^{\star} + e^{-\beta t}\log\hat k(0) \]
\[ \Leftrightarrow \hat k(t) \approx {\hat k^{\star}}^{1-e^{-\beta t}} {\hat k(0)}^{e^{-\beta 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 (traits pleins). 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.
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 = (1-\alpha)(n+x+\delta)$, obtenue en log-linéarisant la dynamique), sous estime la vitesse d'ajustement vers l'état stationnaire.