Simulation 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 :

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

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 :

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

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 :

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

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.

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 :

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 :

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 :

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.

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

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

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.

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.

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 2 la courbe bleue par sa tangente en \(\hat k^{\star}\) (la droite verte dans la figure 4) :

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 5 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.

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.

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 ¹, 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)² on peut calculer la vitesse de convergence à chaque instant. La figure 6 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).

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

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.