diff --git a/.gitignore b/.gitignore
index 938b640..758030c 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1 +1,3 @@
-blog/representation-ma-du-processus-ar2/ltximg/*
\ No newline at end of file
+pages/posts/ltximg/*
+blog/representation-ma-du-processus-ar2/ltximg/*
+blog/simulation-du-modele-de-solow/ltximg/*
diff --git a/Makefile b/Makefile
index 842ff8d..ee7fd37 100644
--- a/Makefile
+++ b/Makefile
@@ -22,6 +22,7 @@ assets:
@rsync --recursive -avz assets/oldies output
@rsync assets/stepan.gpg-pub.asc output
@rsync assets/favicon.ico output
+ @scp -r blog/simulation-du-modele-de-solow/img output/posts/simulation-du-modele-de-solow
clean:
@rm -rf output/posts/representation-ma-du-processus-ar2/ltximg
diff --git a/blog/simulation-du-modele-de-solow/index.org b/blog/simulation-du-modele-de-solow/index.org
new file mode 100644
index 0000000..b9a4ea8
--- /dev/null
+++ b/blog/simulation-du-modele-de-solow/index.org
@@ -0,0 +1,520 @@
+#+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:
+#+HTML_HEAD:
+#+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)
+\]
+
+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 :
+
+#+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
+ # 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))
+#+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 = (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 du stock de capital physique par tête efficace (approximation).*
+#+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é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.
+
+#+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 = (1-\alpha)(n+x+\delta)$, obtenue en log-linéarisant la
+dynamique), sous estime la vitesse d'ajustement vers l'état stationnaire.
+
+# Local Variables:
+# org-confirm-babel-evaluate: nil
+# End:
diff --git a/pages/posts/index.org b/pages/posts/index.org
index c3675da..1098038 100644
--- a/pages/posts/index.org
+++ b/pages/posts/index.org
@@ -21,6 +21,8 @@
- [[./representation-ma-du-processus-ar2][Représentation MA($\infty$) d'un processus AR(2) à racines réelles]]
+- [[./simulation-du-modele-de-solow][Simulation du modèle de Solow]]
+
\\
\\
\\
diff --git a/publish.el b/publish.el
index b950822..f583139 100644
--- a/publish.el
+++ b/publish.el
@@ -7,8 +7,6 @@
(add-to-list 'package-archives '("melpa" . "http://melpa.org/packages/"))
(package-initialize)
-(use-package htmlize
- :ensure t)
(use-package ox-publish
:init
@@ -116,3 +114,29 @@
(add-to-list 'org-export-filter-link-functions 'sa-filter-local-links)
)
+
+ ;(use-package htmlize
+ ; :ensure t);
+ ;
+ ;(setq org-export-htmlize-output-type 'css);
+
+
+(org-babel-do-load-languages
+ 'org-babel-load-languages
+ '(
+ (emacs-lisp . t)
+ (org . t)
+ (shell . t)
+ (python . t)
+ (gnuplot . t)
+ (matlab . t)
+ (R . t)
+ (lilypond . t)
+ ))
+
+(setq org-confirm-babel-evaluate nil)
+;(setq org-src-fontify-natively t)
+(setq org-src-tab-acts-natively t)
+(setq org-babel-python-command "python3")
+(setq org-html-htmlize-output-type 'css)
+(setq org-html-head-include-default-style nil)