#+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 html5-fancy:t #+HTML_DOCTYPE: html5 #+HTML_HEAD: #+HTML_HEAD: #+HTML_HEAD: #+HTML_HEAD: #+HTML_HEAD: #+HTML_HEAD: #+LANGUAGE: fr #+STARTUP: latexpreview #+TITLE: Les vitesses de convergence en panel #+DATE: Septembre 2026 #+AUTHOR: Stéphane Adjemian #+EMAIL: stephane.adjemian@univ-lemans.fr #+PROPERTY: header-args:python :python /tmp/blog-panel/bin/python #+BEGIN_QUOTE La note sur [[https://stephane-adjemian.fr/posts/vitesses-de-convergence-propres-a-chaque-pays/][les vitesses de convergence propres à chaque pays]] s'achevait sur un doute\nbsp{}: le modèle de Solow y était rejeté, mais une simulation montrait que ce rejet pouvait venir de l'endogénéité du revenu initial, qu'une estimation en coupe ne permet pas de lever. Le panel le permet, en laissant chaque pays avoir son propre niveau d'efficacité du travail. Il en coûte pourtant plus cher qu'annoncé\nbsp{}: l'estimateur le plus usuel, l'estimateur intra, ne répond presque pas au paramètre qu'il est censé mesurer, et il faut le montrer par simulation avant de le croire. Cette note introduit le panel pour qui n'en a jamais fait, l'estime, mesure le biais, et le corrige par inférence indirecte, une méthode expliquée pas à pas. Elle répond ensuite aux questions de la note précédente. Les estimations de la vitesse — 0,10 en coupe, 0,48 par l'estimateur intra, 0,18 par le maximum de vraisemblance — sont trois lectures d'une même valeur, voisine de 0,36, soit une convergence de l'ordre de 3\nbsp{}% par an — un ordre de grandeur, pas un chiffre\nbsp{}: la correction est imprécise, et la note dit pourquoi. Le biais de la coupe venait bien de l'endogénéité du revenu initial, mais le rejet du modèle de Solow, lui, survit à sa correction. Une dernière section confronte aux données l'hypothèse sur laquelle tout cela repose\nbsp{}: les erreurs y apparaissent corrélées d'une période à l'autre, non pour une raison économique, mais parce que les données agrègent le temps. Les calculs sont faits en Python, et les corrigés des exercices sont donnés dans des blocs dépliables. #+END_QUOTE \\ \\ \\ #+BEGIN_SRC bash :results silent :exports none :async t python3 -m venv /tmp/blog-panel source /tmp/blog-panel/bin/activate pip install numpy scipy matplotlib pandas openpyxl #+END_SRC #+begin_src python :session panel :exports none :results none import os import urllib.request as request import urllib.error as error import numpy as np import pandas as pd import matplotlib.pyplot as plt from scipy.optimize import least_squares, minimize from scipy import stats #+end_src * De la coupe au panel :PROPERTIES: :CUSTOM_ID: panel :END: On reprend les notations de la note précédente\nbsp{}: $s$ est le taux d'épargne, $n$ le taux de croissance de la population, $x$ celui de l'efficacité du travail, $\delta$ le taux de dépréciation, $\alpha$ l'élasticité de la production au capital, et la vitesse de convergence du pays $i$ est $\beta_i = (1-\alpha)(n_i+x+\delta_i)$. En coupe transversale, on ne dispose que d'une observation par pays, et l'efficacité initiale du travail, $\log A_i(0) = a + \varepsilon_i$, doit être traitée comme une perturbation supposée indépendante du taux d'épargne et du revenu initial. C'est cette hypothèse qui était en cause. \\ Suivre les pays sur plusieurs périodes la rend inutile. Découpons 1963-2023 en périodes de longueur $\tau$, et notons $t_k = k\tau$ la date qui ouvre la période $k$. Rien n'oblige les taux à rester constants d'une période à l'autre\nbsp{}: notons $c_{i,k} = n_{i,k}+x+\delta_{i,k}$ leur somme sur la période. Le modèle de Solow en déduit la vitesse $\beta_{i,k} = (1-\alpha)c_{i,k}$\nbsp{}; comme dans la note précédente, nous n'imposons pas que le coefficient de $c_{i,k}$ soit $1-\alpha$, et écrivons $\beta_{i,k} = \theta\,c_{i,k}$, où $\theta$ est un paramètre libre. Le modèle de Solow complet correspond à $\theta = 1-\alpha$, ce qui fournira un test. Sur un intervalle pendant lequel le taux d'épargne vaut $s_{i,k}$, la solution log-linéarisée du modèle relie exactement la production par tête efficace en fin de période à celle du début\nbsp{}: $\log\hat y_{i,k+1} = \rho_{i,k}\log\hat y_{i,k} + (1-\rho_{i,k})\log\hat y^{\star}_{i,k}$, avec $\rho_{i,k} = e^{-\beta_{i,k}\tau} = e^{-\theta c_{i,k}\tau}$. En revenant au revenu par tête, avec $\log A_i(t) = \log A_i(0) + xt$, on obtient \begin{equation*} \log y_{i,k+1} = \rho_{i,k}\log y_{i,k} + (1-\rho_{i,k})\frac{\alpha}{1-\alpha}\bigl(\log s_{i,k} - \log c_{i,k}\bigr) + g_{i,k}\log A_i(0) + \gamma_{i,k} + v_{i,k+1} \end{equation*} où $g_{i,k} = 1-\rho_{i,k}$ et où \begin{equation*} \gamma_{i,k} = x\bigl[\tau + t_k\,g_{i,k}\bigr] \end{equation*} est un terme de tendance. L'efficacité initiale du travail, qui contenait $\varepsilon_i$ et faussait l'estimation en coupe, n'est plus une perturbation\nbsp{}: c'est une inconnue propre à chaque pays, et l'on ne suppose plus rien sur sa corrélation avec le taux d'épargne ou le revenu. \\ Cette équation peut surprendre qui vient de la note précédente\nbsp{}: la variable expliquée y était un *taux de croissance*, et le revenu initial une explicative. Ici le revenu de fin de période est expliqué par son *niveau* de début de période. Ce n'est pourtant que la même équation, réarrangée\nbsp{}: en retranchant $\log y_{i,k}$ des deux côtés, \begin{equation*} \underbrace{\log y_{i,k+1} - \log y_{i,k}}_{\text{croissance sur la période}} = -g_{i,k}\underbrace{\Bigl[\log y_{i,k} - \frac{\alpha}{1-\alpha}\bigl(\log s_{i,k} - \log c_{i,k}\bigr) - \log A_i(0)\Bigr]}_{\text{écart au sentier stationnaire}} + \gamma_{i,k} + v_{i,k+1} \end{equation*} on retrouve exactement la forme de la note en coupe\nbsp{}: la croissance est proportionnelle à la distance qui reste à parcourir, et le facteur de proportionnalité est $g_{i,k} = 1-e^{-\beta_{i,k}\tau}$. Déplacer $\log y_{i,k}$ d'un membre à l'autre change son coefficient de $\rho$ en $\rho-1$, et rien d'autre\nbsp{}: mêmes résidus, mêmes estimations, à la mêmes résidus et mêmes estimations. Si nous écrivons désormais le niveau, c'est pour une raison de lecture\nbsp{}: sous cette forme l'équation est un processus autorégressif d'ordre un affecté d'un effet pays, et tout ce qui suit — le biais, sa formule, les remèdes — s'énonce en fonction de $\rho$. La forme en taux de croissance le masquerait. \\ Trois remarques sur cette équation, qui commandent tout ce qui suit. D'abord, la tendance n'ajoute aucun paramètre\nbsp{}: $\gamma_{i,k}$ est entièrement connue dès que $\theta$ l'est, alors que la pratique courante consiste à la remplacer par des indicatrices temporelles, qui en coûtent une par période. Ensuite, l'hétérogénéité des vitesses reste gratuite, comme en coupe. Avec les 110 pays et les douze périodes de cinq ans que la section suivante construit, l'équation fait intervenir 1320 coefficients $\rho_{i,k}$, tous différents, mais aucun n'est estimé séparément\nbsp{}: chacun se calcule à partir du seul $\theta$ et du $c_{i,k}$ observé. Les inconnues sont donc $\theta$ et $\alpha$, qui sont l'objet de l'étude, et les 110 efficacités initiales $\log A_i(0)$, une par pays, qu'il faut estimer pour ne pas fausser les deux autres mais qui ne nous intéressent pas en elles-mêmes — on parle de *paramètres de nuisance*, et, lorsque leur nombre croît avec l'échantillon comme ici, de *paramètres incidents*. Enfin, l'équation est non linéaire en $\theta$ et $\alpha$ mais *linéaire* en $\log A_i(0)$, ce dont l'estimation va tirer parti — à ceci près que le coefficient qui la multiplie, $g_{i,k}$, change d'une période à l'autre, ce qui demandera une précaution que nous verrons le moment venu. #+BEGIN_remarque Le terme $v_{i,k+1}$ n'a pas été défini, et il le faut, car les deux estimateurs de cette note reposent sur ce qu'on en suppose. Puisque $\log A_i(0)$ absorbe désormais le *niveau* de l'efficacité du travail, $v_{i,k+1}$ est tout le reste de ce qui déplace le revenu entre $k$ et $k+1$\nbsp{}: les innovations de l'efficacité du travail autour de sa tendance, l'erreur de l'approximation log-linéaire lorsque le pays est loin de son état stationnaire, l'erreur de mesure du revenu, et la variation des taux à l'intérieur d'une période, que les moyennes écrasent. Nous le supposerons d'espérance nulle, de variance $\sigma^2$ constante, et d'abord non corrélé dans le temps. C'est une hypothèse forte, et c'est elle qui donne sa forme à la matrice de covariance de la section sur [[#condition][le maximum de vraisemblance]]\nbsp{}: si $v$ était autocorrélé, la vraisemblance que nous y écrirons serait fausse. Nous verrons qu'il l'est, et pour quelle raison, et nous l'autoriserons alors à suivre une moyenne mobile. #+END_remarque #+BEGIN_remarque L'inconnue de pays n'est pas seulement une nuisance à éliminer\nbsp{}: elle est interprétable. On récupère $\log\hat A_i(0)$ pour chaque pays, et on peut la confronter au taux d'épargne et à la croissance démographique, c'est-à-dire tester l'hypothèse que Mankiw, Romer et Weil posent et que la coupe transversale obligeait à admettre. C'est l'objet de la section [[#retour][Retour à la coupe transversale]]. #+END_remarque * Construire le panel :PROPERTIES: :CUSTOM_ID: donnees :END: Nous reprenons les Penn World Tables 11.0 et l'échantillon de la note précédente\nbsp{}: les 110 pays pour lesquels le PIB, la population, la part de l'investissement et le taux de dépréciation sont renseignés toutes les années de 1960 à 2023. #+begin_src python :session panel :exports none :results none PWT = 'pwt110.xlsx' class SansRedirection(request.HTTPRedirectHandler): def redirect_request(self, *args, **kwargs): return None def adresse(url): """Adresse de téléchargement directe d'un fichier DataverseNL.""" ouvreur = request.build_opener(SansRedirection) try: with ouvreur.open(request.Request(url, method='HEAD')) as reponse: return reponse.headers.get('Location') or url except error.HTTPError as echec: return echec.headers.get('Location') or url if not os.path.exists(PWT): request.urlretrieve(adresse('https://dataverse.nl/api/access/datafile/554105'), PWT) #+end_src La période retenue commence en 1963 et non en 1960\nbsp{}: soixante années se découpent exactement en douze périodes de cinq ans, et aussi en six périodes de dix ans, dont la section sur les limites du modèle statistique fera un usage ponctuel. #+begin_src python :session panel :exports code :results none DEBUT, FIN, X = 1963, 2023, 0.02 pwt = pd.read_excel('pwt110.xlsx', 'Data') brut = (pwt[(pwt['year'] >= 1960) & (pwt['year'] <= 2023)] .set_index(['countrycode', 'year'])[['rgdpo', 'pop', 'csh_i', 'delta']] .sort_index()) complet = brut.notna().groupby('countrycode').sum() == 64 codes = complet[complet.all(axis=1)].index positif = brut.loc[codes].groupby('countrycode')['csh_i'].mean() > 0 codes = positif[positif].index don = brut.loc[codes] annees = np.arange(DEBUT, FIN + 1) tableau = lambda v: don[v].unstack(0).loc[annees, codes].to_numpy() ly_an = np.log(tableau('rgdpo')/tableau('pop')) s_an, pop_an, delta_an = tableau('csh_i'), tableau('pop'), tableau('delta') N = len(codes) def panel(tau): """Revenus aux dates de découpage, et taux moyens sur chaque période.""" dates = np.arange(0, FIN - DEBUT + 1, tau) K = len(dates) - 1 LY = ly_an[dates, :] S = np.array([s_an[dates[k]:dates[k+1], :].mean(axis=0) for k in range(K)]) D = np.array([delta_an[dates[k]:dates[k+1], :].mean(axis=0) for k in range(K)]) NN = np.array([(pop_an[dates[k+1]]/pop_an[dates[k]])**(1/tau) - 1 for k in range(K)]) return LY, S, NN + X + D, K LY, S, C, K = panel(5) print(N, K, np.median(C), C.std(axis=0).mean(), C.mean(axis=0).std()) #+end_src L'échantillon compte 110 pays et douze périodes. Le $c_{i,k}$ médian vaut 8,1\nbsp{}% par an. Sa dispersion est pour l'essentiel entre les pays plutôt que dans le temps — écart-type de 1,53 point entre pays contre 0,81 point d'une période à l'autre à l'intérieur d'un pays — mais nous n'avons pas besoin de le supposer\nbsp{}: la section sur [[#intra][l'estimateur intra]] traite le cas général, celui de taux qui varient d'une période à l'autre. * Ce qu'un panel ajoute à une coupe :PROPERTIES: :CUSTOM_ID: panel-abc :END: Les notes précédentes travaillaient sur une *coupe transversale*\nbsp{}: une observation par pays, 110 lignes. Nous disposons maintenant d'un *panel*\nbsp{}: 110 pays observés à 13 dates, soit 1320 transitions d'une période à la suivante. Ce n'est pas seulement douze fois plus de données. C'est une seconde dimension, et elle change la nature du problème. \\ Toute grandeur observée dans un panel varie de deux façons. Elle varie *entre* les pays — le revenu français n'est pas le revenu indien — et elle varie *à l'intérieur* de chaque pays au fil du temps — la France de 2023 n'est pas la France de 1963. On appelle la première variation /inter/, la seconde /intra/, et la variance totale se partage exactement entre les deux. #+begin_src python :session panel :exports code :results none moyennes = LY.mean(axis=0) # une moyenne temporelle par pays v_tot, v_inter, v_intra = LY.var(), moyennes.var(), (LY - moyennes).var() print(round(v_tot, 3), round(v_inter, 3), round(v_intra, 3), round(v_inter/v_tot, 3), round(v_intra/v_tot, 3)) fra = LY[:, list(codes).index('FRA')] print(fra[[0, -1]].round(2), fra.mean().round(2), (fra - fra.mean())[[0, -1]].round(2)) croissance = (LY[1:, :] - LY[:-1, :])/5 # taux de croissance annuel moyen par période m_g = croissance.mean(axis=0) print(round(m_g.var()/croissance.var(), 3), round((croissance - m_g).var()/croissance.var(), 3), round(100*m_g.std(), 2), round(100*(croissance - m_g).std(), 2), round(croissance.mean(axis=1).var()/croissance.var(), 3)) def pente(a, b): """Pente des MCO de b sur a, toutes observations empilées.""" A = np.column_stack([np.ones(a.size), a.ravel()]) return np.linalg.lstsq(A, b.ravel(), rcond=None)[0][1] avant, apres = LY[:-1, :], LY[1:, :] p_confondu = pente(avant, apres) p_intra = pente(avant - avant.mean(axis=0), apres - apres.mean(axis=0)) print(round(p_confondu, 4), round(p_intra, 4), round(-np.log(p_intra)/5*100, 2)) #+end_src Dans nos données, 81\nbsp{}% de la variance du revenu par tête est /inter/ et 19\nbsp{}% seulement est /intra/. Autrement dit, les quatre cinquièmes de ce que l'on observe tiennent au fait que les pays sont durablement différents, et un cinquième au fait que chacun bouge dans le temps. **La coupe transversale n'exploitait que la première dimension** — et c'est exactement là que loge le problème, puisque l'efficacité du travail $\log A_i(0)$, inobservée, est une caractéristique durable du pays\nbsp{}: elle est tout entière dans la variation /inter/, et rien ne permettait de l'y distinguer du reste. \\ Le partage s'inverse lorsqu'on regarde le taux de croissance du revenu par tête sur chaque période plutôt que son niveau\nbsp{}: 14\nbsp{}% seulement de sa variance est /inter/, et 86\nbsp{}% est /intra/. L'écart-type des taux de croissance moyens des pays est de 1,65 point par an, quand celui des écarts de chaque période à la moyenne de son pays atteint 4,15 points. Les pays diffèrent durablement par leur *niveau* de revenu, beaucoup moins par leur *croissance*, qui fluctue surtout d'une période à l'autre — et les périodes y comptent elles-mêmes pour 7,5\nbsp{}%, signe de chocs communs à tous les pays. Or la variable expliquée de la note en coupe, le taux de croissance moyen de 1963 à 2023, n'est rien d'autre que la moyenne des taux de croissance de période d'un pays\nbsp{}: la coupe ne cherchait à expliquer que ces 14\nbsp{}%, et laissait de côté les 86\nbsp{}% restants, qui sont précisément ce que le panel exploite. C'est aussi à cette variation qu'il faudra rapporter l'ajustement du modèle. \\ La variation /intra/, elle, en est affranchie, et pour une raison presque triviale\nbsp{}: $\log A_i(0)$ est le *même nombre* à toutes les dates d'un pays donné. Comparer un pays à lui-même l'élimine donc, sans qu'on ait besoin de le connaître ni de supposer quoi que ce soit sur lui. C'est tout le principe. Concrètement, pour la France, dont le revenu par tête en logarithme passe de 9,54 en 1963 à 10,85 en 2023, de moyenne 10,33, on ne travaille plus sur les niveaux mais sur les écarts à cette moyenne, de −0,79 au début à +0,52 à la fin. Retrancher à chaque pays sa propre moyenne temporelle s'appelle la *transformation intra*, et l'estimateur qui en résulte l'*estimateur intra* — /within estimator/ en anglais. #+BEGIN_remarque Le vocabulaire du panel est un petit obstacle qu'il vaut mieux lever d'emblée. Une variable propre au pays et constante dans le temps, comme $\log A_i(0)$, s'appelle un *effet fixe* quand on la traite comme un paramètre à estimer, et un *effet aléatoire* quand on la traite comme une perturbation — c'est ce que faisait la coupe transversale. Un panel est dit *dynamique* lorsque la variable expliquée figure aussi, retardée, parmi les explicatives, ce qui est notre cas puisque $\log y_{i,k}$ explique $\log y_{i,k+1}$. Enfin, les $N$ effets fixes sont appelés *paramètres incidents*\nbsp{}: leur nombre croît avec l'échantillon, et c'est de là que viendront tous nos ennuis. #+END_remarque \\ Une expérience simple montre que les deux dimensions ne racontent pas la même histoire. Régressons le revenu d'une période sur celui de la précédente, d'abord en mélangeant toutes les observations, puis après avoir retranché à chaque pays sa moyenne. La pente confondue vaut 1,002\nbsp{}: aucune convergence inconditionnelle, les pays riches le restent. La pente intra vaut 0,906, ce qui, lu naïvement, ferait une vitesse de 2,0\nbsp{}% par an. #+begin_src python :session panel :exports none :results none fig, ax = plt.subplots(figsize=(6, 4.5)) ax.scatter(avant.ravel(), apres.ravel(), s=6, color='lightgrey', label='les 1320 transitions') couleurs = {'FRA': 'r', 'USA': 'b', 'JPN': 'g', 'IND': 'darkorange', 'NGA': 'purple'} for pays, couleur in couleurs.items(): j = list(codes).index(pays) ax.scatter(avant[:, j], apres[:, j], s=18, color=couleur, label=pays) bornes = np.array([avant.min(), avant.max()]) ax.plot(bornes, bornes, 'k:', linewidth=1, label='première bissectrice') ax.set_xlabel(r'$\log y_{i,k}$, début de période') ax.set_ylabel(r'$\log y_{i,k+1}$, fin de période') ax.legend(fontsize=7, ncol=2) fig.tight_layout() fig.savefig("panel-nuage.svg", transparent=True) #+end_src #+CAPTION: *Les 1320 transitions du panel, et cinq pays mis en évidence.* Le nuage d'ensemble s'aligne sur la bissectrice — entre pays, rien ne converge. À l'intérieur de chaque pays, les points montent régulièrement le long d'une trajectoire propre : c'est cette variation-là que le panel exploite. #+LABEL: fig:nuage [[file:panel-nuage.svg]] Cette lecture naïve serait pourtant une erreur, et il faut s'en convaincre tout de suite. Le revenu par tête n'est pas une variable stationnaire\nbsp{}: il croît, en espérance et en variance, et rien n'assure a priori qu'il revienne vers quoi que ce soit. Prenons donc le cas extrême d'un monde où *aucun* pays ne converge\nbsp{}: chaque revenu suit une marche aléatoire, avec la croissance moyenne observée et des chocs de la dispersion observée, à partir des revenus de 1963. Puis refaisons les deux régressions. #+begin_src python :session panel :exports code :results none dy = LY[1:, :] - LY[:-1, :] derive, choc = dy.mean(), (dy - dy.mean(axis=0)).std() alea = np.random.default_rng(5) marche = [] for _ in range(1000): Y = LY[0] + np.vstack([np.zeros(N), np.cumsum(derive + alea.normal(0, choc, (K, N)), axis=0)]) av, ap = Y[:-1, :], Y[1:, :] marche.append([pente(av, ap), pente(av - av.mean(axis=0), ap - ap.mean(axis=0))]) marche = np.array(marche) print(marche.mean(axis=0).round(3), np.percentile(marche[:, 1], [5, 95]).round(3), (marche[:, 1] <= p_intra).mean().round(2)) #+end_src La pente confondue vaut en moyenne 1,000, et la pente intra 0,906 — les deux chiffres observés, alors que ces données ne contiennent aucune convergence. La pente intra tombe entre 0,884 et 0,926 dans 90\nbsp{}% des cas, et la valeur observée est au milieu de cette distribution. Autrement dit, l'expérience ne permet pas de distinguer nos données d'un monde sans convergence. La remarque qui suit explique pourquoi retrancher la moyenne d'un pays fait mécaniquement baisser la pente, même quand la vraie vaut un, et [[ex-marche][un exercice]] en donne la formule exacte. \\ Ce que l'expérience montre, c'est seulement que les deux dimensions des données donnent des réponses différentes, et que la dimension intra ne peut pas être lue sans précaution. Pour établir une convergence, il faudra deux choses que cette régression n'a pas\nbsp{}: les déterminants de l'état stationnaire de chaque pays, que le modèle de Solow fournit, et une correction du biais, qui occupera la majeure partie de la note. #+BEGIN_remarque Un avertissement, pour éviter un contresens fréquent. La transformation intra est inoffensive dans un panel statique\nbsp{}: retrancher une moyenne ne corrèle pas les régresseurs à l'erreur si ceux-ci sont donnés d'avance. Elle cesse de l'être dès que la variable retardée figure au membre de droite, car la moyenne que l'on retranche contient, elle aussi, les perturbations de toutes les périodes — y compris celles qui ont déterminé le revenu retardé. Le régresseur transformé et l'erreur transformée deviennent alors corrélés. C'est le biais de Nickell, et il occupera la section sur [[#biais][le biais de l'estimateur intra]]. #+END_remarque * L'estimateur intra :PROPERTIES: :CUSTOM_ID: intra :END: Passons de l'illustration au modèle. La transformation intra de la section précédente retranchait à chaque pays sa moyenne temporelle\nbsp{}; il va falloir la généraliser, parce que notre équation n'est pas celle d'un panel ordinaire. Notons $r_{i,k}(\theta,\alpha)$ le résidu avant retrait de l'inconnue de pays, c'est-à-dire tout ce qui reste une fois la partie prédite par le modèle soustraite. L'équation s'écrit alors $r_{i,k} = g_{i,k}\log A_i(0) + v_{i,k+1}$\nbsp{}: l'inconnue de pays y est bien la même à toutes les dates, mais elle est multipliée par un coefficient $g_{i,k}$ qui, lui, change d'une période à l'autre. Retrancher la moyenne simple ne l'éliminerait donc pas. Que faire\nbsp{}? Regarder l'équation $r_{i,k} = g_{i,k}\log A_i(0) + v_{i,k+1}$ pour ce qu'elle est\nbsp{}: une régression, pays par pays, du résidu $r_{i,k}$ sur le coefficient $g_{i,k}$, et sans constante. L'inconnue $\log A_i(0)$ en est la pente. Pour $\theta$ et $\alpha$ donnés, la minimisation en $\log A_i(0)$ donne donc la pente des moindres carrés ordinaires de $r_{i,k}$ sur $g_{i,k}$\nbsp{}: \begin{equation*} \widehat{\log A_i(0)} = \frac{\sum_{k} g_{i,k}\,r_{i,k}}{\sum_{k} g_{i,k}^2} \end{equation*} et le critère concentré est la somme des carrés des $r_{i,k} - g_{i,k}\widehat{\log A_i(0)}$. On peut aussi la lire comme une moyenne pondérée\nbsp{}: en divisant l'équation par $g_{i,k}$, chaque période propose l'estimation $r_{i,k}/g_{i,k}$, et la formule les moyenne en donnant à chacune le poids $g_{i,k}^2$ — les périodes où le pays converge le plus vite, donc celles où l'efficacité initiale pèse le plus, comptent davantage. Lorsque les taux sont constants, $g_{i,k}$ ne dépend plus de $k$, la formule se réduit à la moyenne temporelle des résidus, et l'on retrouve la transformation intra habituelle. C'est donc bien une généralisation, et non un autre estimateur. \\ L'estimation elle-même se fait par les moindres carrés non linéaires (MCNL), comme dans la note en coupe\nbsp{}: on minimise sur $\theta$ et $\alpha$ la somme des carrés de ces résidus concentrés, par un algorithme de type Gauss-Newton — qui remplace à chaque itération le modèle par son approximation linéaire et applique les MCO — lancé depuis plusieurs points de départ. Comme l'équation est linéaire en $\log A_i(0)$, concentrer ainsi les 110 inconnues de pays donne exactement les mêmes estimations que les MCNL sur les 112 paramètres à la fois, pour un problème à deux dimensions seulement. Les écarts-types proviennent de la formule habituelle des MCNL, $\hat\sigma^2(J'J)^{-1}$, où $J$ est la matrice des dérivées des résidus concentrés par rapport aux deux paramètres — la /jacobienne/, qui joue ici le rôle de la matrice des régresseurs, et le test de la contrainte compare les sommes des carrés des résidus avec et sans elle. #+BEGIN_remarque Cette régression n'est pas un détail d'écriture. Si l'on retranchait la moyenne simple des résidus alors que le coefficient de $\log A_i(0)$ varie, on n'éliminerait pas l'inconnue de pays, et l'estimation serait fausse\nbsp{}: $\hat\theta$ tomberait à 0,07. Et elle doit être refaite *à l'intérieur* de la fonction qui calcule les résidus\nbsp{}: son régresseur $g_{i,k}$ dépend de $\theta$, donc il change à chaque itération de l'optimiseur. Centrer les données en amont, comme le ferait un réflexe de panel linéaire, donne un estimateur différent, et incorrect. #+END_remarque #+begin_src python :session panel :exports code :results none def bruts(p, LY, S, C, K, tau, contraint=False): """Résidus r_{i,k} avant retrait de l'inconnue de pays, et coefficients g_{i,k}.""" alpha = p[0] if contraint else p[1] theta = 1 - alpha if contraint else p[0] rho = np.exp(-theta*C*tau) tk = (np.arange(K)*tau)[:, None] pred = (rho*LY[:-1, :] + (1 - rho)*alpha/(1 - alpha)*(np.log(S) - np.log(C)) + X*(tau + tk*(1 - rho))) return LY[1:, :] - pred, 1 - rho def intra(LY, S, C, K, tau, contraint=False): """MCNL, log A_i(0) concentré par la régression de r sur g.""" def critere(p): r, g = bruts(p, LY, S, C, K, tau, contraint) return (r - g*(g*r).sum(axis=0)/(g*g).sum(axis=0)).ravel() if contraint: departs, bornes = [[0.4], [0.6], [0.8]], ([0.01], [0.99]) else: departs, bornes = [[0.4, 0.4], [0.1, 0.7], [0.8, 0.2], [0.03, 0.5]], \ ([1e-4, 0.01], [5, 0.99]) essais = [least_squares(critere, p0, bounds=bornes) for p0 in departs] return min(essais, key=lambda e: e.cost) def resultats(tau): """Estimation libre, estimation contrainte, écarts-types et test.""" A, B, Cx, Kx = panel(tau) f, fc = intra(A, B, Cx, Kx, tau), intra(A, B, Cx, Kx, tau, contraint=True) ddl = Kx*N - N - 2 V = np.linalg.inv(f.jac.T @ f.jac)*(f.fun @ f.fun)/ddl sigma = np.sqrt(f.fun @ f.fun/((Kx - 1)*N)) F = (2*fc.cost - 2*f.cost)/(2*f.cost/ddl) return f, np.sqrt(np.diag(V)), sigma, F, Kx, Cx f, se, sigma, F, Kx, Cx = resultats(5) print(Kx, f.x.round(4), se.round(4), np.median(f.x[0]*Cx.mean(axis=0)).round(4), sigma.round(4), round(F, 2), round(f.x[0] - (1 - f.x[1]), 3)) #+end_src Avec des périodes de cinq ans, l'estimation donne $\hat\theta = 0{,}478$ (écart-type 0,048) et $\hat\alpha = 0{,}370$ (0,037), soit une vitesse de convergence médiane de 3,9\nbsp{}% par an. C'est plus de quatre fois la vitesse estimée en coupe, 0,86\nbsp{}%, et l'ordre de grandeur que rapporte Islam (1995). La contrainte du modèle de Solow complet, $\theta = 1-\alpha$, est rejetée\nbsp{}: la statistique de Fisher vaut 12,6. #+begin_src python :session panel :exports code :results none y = LY[1:, :].ravel() D = np.tile(np.eye(N), (K, 1)) T = np.repeat(np.eye(K)[:, 1:], N, axis=0) Z = np.column_stack([LY[:-1, :].ravel(), np.log(S).ravel(), np.log(C).ravel(), D, T]) b = np.linalg.lstsq(Z, y, rcond=None)[0] print(b[0], -np.log(b[0])/5, b[1], b[1]/(1 - b[0] + b[1])) def ajustement(e, k): """SCR, R² sur la croissance intra, AIC et BIC d'une spécification à k paramètres.""" dy_intra = (LY[1:, :] - LY[:-1, :]) - (LY[1:, :] - LY[:-1, :]).mean(axis=0) scr, n = e @ e, e.size ll = -n/2*np.log(scr/n) return [round(scr, 2), k, round(1 - scr/(dy_intra**2).sum(), 3), round(-2*ll + 2*k, 1), round(-2*ll + k*np.log(n), 1)] r_s, g_s = bruts(resultats(5)[0].x, LY, S, C, K, 5) e_s = (r_s - g_s*(g_s*r_s).sum(axis=0)/(g_s*g_s).sum(axis=0)).ravel() Z0 = Z[:, :3 + N] b0 = np.linalg.lstsq(Z0, y, rcond=None)[0] print(ajustement(e_s, 2 + N)) print(ajustement(y - Z @ b, Z.shape[1])) print(ajustement(y - Z0 @ b0, Z0.shape[1]), round(b0[0], 3)) #+end_src À titre de repère, la régression linéaire habituelle, celle qui traite $\rho$ comme un coefficient commun et ajoute des indicatrices de pays et de période, donne $\hat\rho = 0{,}842$, soit 3,4\nbsp{}% par an, et une part du capital implicite de 0,44. \\ Comparer l'ajustement des deux approches est instructif, à condition de séparer ce qui les distingue. Le tableau rapporte la somme des carrés des résidus (SCR), la part de la variation intra du *taux de croissance* qui est expliquée, et deux critères d'information qui pénalisent le nombre de paramètres, ceux d'Akaike (AIC) et de Schwarz (BIC)\nbsp{}: plus ils sont bas, meilleure est la spécification. #+ATTR_HTML: :class numbers | spécification | paramètres | SCR | $R^2$ | AIC | BIC | |----------------------------------------------+------------+-------+-------+---------+---------| | structurelle | 112 | 50,99 | 0,102 | −4071,0 | −3490,2 | | linéaire, indicatrices de pays et de période | 124 | 45,51 | 0,199 | −4197,2 | −3554,2 | | linéaire, indicatrices de pays seulement | 113 | 50,62 | 0,109 | −4078,5 | −3492,6 | La régression habituelle ajuste nettement mieux\nbsp{}: sa SCR est plus basse de 11\nbsp{}%, et les deux critères la préfèrent. Mais la dernière ligne montre d'où vient l'écart. Sans indicatrices de période, la régression linéaire, qui laisse libres les coefficients de $\log s$ et de $\log c$, n'ajuste pas mieux que le modèle structurel, et son $\hat\rho$ monte à 0,884. Imposer la structure de Solow au lieu de coefficients libres ne coûte donc presque rien, et rapporte un test. Ce qui coûte, c'est de remplacer des chocs communs à chaque période par la tendance déterministe du modèle. La section sur [[#limites][les limites du modèle statistique]] y reviendra. #+begin_src python :session panel :exports code :results none f5b = resultats(5)[0] r5b, g5b = bruts(f5b.x, LY, S, C, K, 5) e5b = r5b - g5b*(g5b*r5b).sum(axis=0)/(g5b*g5b).sum(axis=0) for cible in (LY[1:, :], LY[1:, :] - LY[:-1, :]): variation = cible - cible.mean(axis=0) print(round(1 - (e5b*e5b).sum()/(variation*variation).sum(), 3)) #+end_src #+BEGIN_remarque Voici le seul endroit où l'écriture en niveau et l'écriture en taux de croissance ne se valent pas, et il vaut d'être signalé. Les résidus sont les mêmes, mais la part de variation *expliquée* dépend de ce qu'on prétend expliquer\nbsp{}: rapportée à la variation intra du niveau, elle vaut 0,84, alors que le tableau ci-dessus, rapporté au taux de croissance, donnait 0,10. Le premier chiffre ne dit rien du modèle — le revenu est si persistant que prédire celui de la période suivante est presque gratuit — et flatterait n'importe quelle spécification. Le second est le chiffre honnête, et il est modeste. C'est pourquoi un $R^2$ n'a de sens, dans un panel dynamique, que sur la grandeur que le modèle prétend réellement expliquer. #+END_remarque * Le biais de l'estimateur intra :PROPERTIES: :CUSTOM_ID: biais :END: Reste à savoir ce que vaut $\hat\theta = 0{,}478$. Retrancher à chaque pays sa moyenne pondérée, dans une équation où le revenu retardé figure au membre de droite, corrèle le régresseur transformé et l'erreur transformée\nbsp{}: c'est le biais que Nickell (1981) a calculé pour le modèle linéaire, d'ordre $-(1+\rho)/(K-1)$ sur $\rho$, et la note sur [[https://stephane-adjemian.fr/posts/biais-de-l-estimateur-mco-d-un-processus-ar1/][le biais de l'estimateur des MCO d'un AR(1)]] en expose le mécanisme. Un $\hat\rho$ trop petit est une vitesse trop grande. Mais l'ordre de grandeur de Nickell suppose un modèle linéaire, et il va se révéler très optimiste. \\ Nous pouvons mesurer ce biais par simulation. Le modèle fournit son propre générateur de données\nbsp{}: on se donne $\theta$, $\alpha$, l'écart-type $\sigma$ des perturbations, le paramètre de moyenne mobile $\psi$ introduit plus loin et les efficacités $\log A_i(0)$, on part des revenus observés en 1963, et on applique l'équation période après période. Les efficacités nous manquent encore\nbsp{}: nous les récupérons de l'estimation par la même régression de $r_{i,k}$ sur $g_{i,k}$, celle qui les a concentrées hors du critère. #+begin_src python :session panel :exports code :results none f5 = intra(LY, S, C, K, 5) theta5, alpha5 = f5.x r5, g5 = bruts(f5.x, LY, S, C, K, 5) lA0 = (g5*r5).sum(axis=0)/(g5*g5).sum(axis=0) cbar = C.mean(axis=0) print(lA0.std().round(3), np.median(cbar).round(4)) #+end_src Elles se dispersent largement d'un pays à l'autre — un écart-type de 1,11 en logarithme, soit un rapport de trois entre un pays typiquement efficace et un autre — et nous les traiterons dans les simulations comme des caractéristiques données, au même titre que les taux observés. #+begin_src python :session panel :exports code :results none def simule(theta, alpha, sigma, psi, lA0, LY, S, C, K, tau, rng): """Un panel simulé par le modèle, erreurs MA(1), efficacités et revenus initiaux donnés.""" rho = np.exp(-theta*C*tau) g = 1 - rho Z = g*(alpha/(1 - alpha)*(np.log(S) - np.log(C)) + lA0) eps = rng.normal(0, sigma, (K + 1, len(lA0))) v = eps[1:] + psi*eps[:-1] Y = np.empty_like(LY) Y[0] = LY[0] for k in range(K): Y[k+1] = rho[k]*Y[k] + Z[k] + X*(tau + k*tau*g[k]) + v[k] return Y def esperance(theta, alpha, sigma, psi, lA0, LY, S, C, K, tau, B, rng): """Moyenne de l'estimateur intra sur B panels simulés.""" return np.array([intra(simule(theta, alpha, sigma, psi, lA0, LY, S, C, K, tau, rng), S, C, K, tau).x for _ in range(B)]) #+end_src L'écart-type des perturbations estimé sur les données vaut 0,20, et c'est à ce niveau que l'expérience a un sens\nbsp{}: un estimateur peut fort bien être excellent dans un monde presque déterministe et muet dans le nôtre. Précisons au passage ce que nous cherchons ici. Ce n'est pas la convergence de l'estimateur, qui se joue à écart-type donné quand la taille de l'échantillon grandit, et qui fait l'objet du dernier des [[#exercices][exercices]]\nbsp{}; c'est son biais à la taille de notre échantillon, 110 pays et douze périodes. \\ Le tableau suivant donne l'expérience au niveau de bruit des données. Deux valeurs y sont empruntées à la section suivante, faute de mieux\nbsp{}: la part du capital, $\alpha = 0{,}508$, et le paramètre de moyenne mobile, $\psi = 0{,}294$, tous deux estimés plus loin. Simuler demande un jeu complet de paramètres, et ce sont les seuls dont nous disposions. #+begin_src python :session panel :exports code :results none GRILLE = (0.10, 0.176, 0.30, 0.60) des_biais = np.random.default_rng(2014) courbe = {th: esperance(th, 0.508, 0.2045, 0.294, lA0, LY, S, C, K, 5, 200, des_biais) for th in GRILLE} for th in GRILLE: print(th, courbe[th][:, 0].mean().round(3), courbe[th][:, 0].std().round(3)) #+end_src #+ATTR_HTML: :class numbers | vrai $\theta$ | 0,10 | 0,176 | 0,30 | 0,60 | | $\mathbb{E}[\hat\theta]$ | 0,451 | 0,456 | 0,492 | 0,686 | La courbe est *plate*\nbsp{}: entre 0,10 et 0,176 — l'intervalle où se situe, nous le verrons, l'estimation par maximum de vraisemblance — l'espérance de $\hat\theta$ ne bouge presque pas, et il faut tripler la vraie valeur pour la voir monter de 0,451 à 0,492. Sur tout cet intervalle, l'estimateur intra ne distingue rien. Ce n'est plus un biais que l'on corrigerait en retranchant une constante, c'est une statistique qui ne répond presque pas au paramètre. L'écart-type de $\hat\theta$, 0,048, n'a dans ces conditions aucune valeur informative. #+begin_src python :session panel :exports none :results none fig, ax = plt.subplots(figsize=(6, 4.5)) vrais = np.array(GRILLE) moyennes = np.array([courbe[th][:, 0].mean() for th in GRILLE]) ax.plot(vrais, moyennes, 'o-', color='r', label=r'espérance de $\hat\theta$') ax.plot(vrais, vrais, 'k--', linewidth=1, label='estimateur sans biais') ax.axhline(theta5, color='grey', linewidth=1, label='estimation sur les données') ax.set_xlabel(r'vraie valeur de $\theta$') ax.set_ylabel(r'espérance de $\hat\theta$') ax.legend(fontsize=8) fig.tight_layout() fig.savefig("panel-biais.svg", transparent=True) #+end_src #+CAPTION: *Espérance de l'estimateur intra en fonction de la vraie valeur du paramètre, au niveau de bruit des données.* La courbe devrait suivre la diagonale ; elle est plate. #+LABEL: fig:biais [[file:panel-biais.svg]] * Un estimateur par le maximum de vraisemblance :PROPERTIES: :CUSTOM_ID: condition :END: L'estimateur intra ne répondant presque pas au paramètre, il en faut un autre. La voie classique pour se débarrasser d'un effet fixe sans retrancher de moyenne consiste à *différencier*\nbsp{}: comparer chaque période à la précédente plutôt qu'à la moyenne de toutes. Mais ici la différence simple $r_{i,k} - r_{i,k-1}$ ne supprime pas l'inconnue de pays, puisque son coefficient change d'une période à l'autre. Procédons donc en deux temps, lentement, car c'est le cœur technique de la note. \\ *Premier temps, normaliser.* Divisons chaque résidu par son propre coefficient\nbsp{}: \begin{equation*} u_{i,k} \equiv \frac{r_{i,k}}{g_{i,k}} = \log A_i(0) + \frac{v_{i,k+1}}{g_{i,k}} \end{equation*} L'inconnue de pays apparaît maintenant avec le même coefficient, un, à toutes les périodes\nbsp{}; en échange, l'erreur a été déformée, et sa variance vaut désormais $\sigma^2 h_{i,k}$ avec $h_{i,k} = 1/g_{i,k}^2$. *Second temps, différencier.* Les écarts successifs $z_{i,k} = u_{i,k} - u_{i,k-1}$ ne contiennent plus $\log A_i(0)$ du tout, pour $k = 2,\ldots,K$. C'est, à un facteur $1/g_{i,k}$ près, la quasi-différence $r_{i,k} - (g_{i,k}/g_{i,k-1})\,r_{i,k-1}$ que Holtz-Eakin, Newey et Rosen (1988) ont proposée pour les effets fixes à coefficient variable. Restent $K-1$ quantités débarrassées de l'inconnue, et une première observation, $z_{i,1} = r_{i,1}$, qui la contient encore. Plutôt que de la jeter, on suit Hsiao, Pesaran et Tahmiscioglu (2002)\nbsp{}: on lui donne une moyenne libre — projetée sur le taux d'épargne initial et sur $\log c_{i,1}$ — et une variance libre $\omega\sigma^2$, qui absorbe ce que l'on ignore de la condition initiale. \\ *Troisième temps, dire ce que l'on suppose sur les erreurs.* L'indépendance des $v_{i,k}$ serait commode, mais la façon dont les données sont construites la rend peu plausible. Les taux d'épargne et de dépréciation sont des moyennes de période, et le produit intérieur brut est lui-même un flux mesuré sur l'année\nbsp{}: nos observations sont donc des agrégats temporels, et moyenner un flux crée une corrélation positive entre périodes voisines, même lorsque le processus sous-jacent n'en a aucune. C'est un résultat ancien, dû à Working (1960), qui donne +0,25 au premier retard pour une marche aléatoire agrégée. Plutôt que d'imposer l'indépendance, laissons donc au modèle la liberté correspondante, celle d'une moyenne mobile d'ordre un\nbsp{}: \begin{equation*} v_{i,k} = \varepsilon_{i,k} + \psi\,\varepsilon_{i,k-1},\qquad \varepsilon_{i,k} \sim \mathcal{N}(0,\sigma^2)\ \text{indépendantes} \end{equation*} de sorte que $\mathbb{V}(v) = \sigma^2(1+\psi^2)$ et $\mathrm{Cov}(v_{i,k},v_{i,k-1}) = \sigma^2\psi$. L'indépendance reste le cas particulier $\psi = 0$, et les données diront ce qu'il en est. \\ Écrivons maintenant la loi du vecteur $z_i = (z_{i,1},\ldots,z_{i,K})$. Les deux opérations, normaliser puis différencier, sont linéaires\nbsp{}: en notant $T_i$ la matrice qui laisse la première ligne intacte et forme les écarts $r_{i,k}/g_{i,k} - r_{i,k-1}/g_{i,k-1}$ pour $k\geq2$, on a $z_i = T_i\,r_i$, et donc \begin{equation*} \mathbb{V}(z_i) = \sigma^2\Sigma_i,\qquad \Sigma_i = T_i\,V\,T_i'\ \text{ où }\ V = \begin{pmatrix} 1+\psi^2 & \psi & & \\ \psi & 1+\psi^2 & \psi & \\ & \psi & 1+\psi^2 & \ddots \\ & & \ddots & \ddots \end{pmatrix} \end{equation*} est la covariance des $v$, à $\sigma^2$ près. Une seule entrée échappe à cette formule\nbsp{}: celle de la première observation, dont la variance est le paramètre libre $\omega\sigma^2$ et non ce que le modèle prédirait. Pour lire cette matrice, le plus simple est de regarder le cas $\psi = 0$, où elle est tridiagonale\nbsp{}: \begin{equation*} \Sigma_i = \begin{pmatrix} \omega & -1/g_{i,1} & & \\ -1/g_{i,1} & h_{i,2}+h_{i,1} & -h_{i,2} & \\ & -h_{i,2} & h_{i,3}+h_{i,2} & \ddots \\ & & \ddots & \ddots \end{pmatrix} \end{equation*} Chaque écart partage une erreur avec son voisin, d'où les termes hors diagonale, et la matrice est *propre à chaque pays* puisque les $h$ le sont. Avec $\psi \neq 0$, les erreurs se recouvrent d'une période de plus et $\Sigma_i$ cesse d'être tridiagonale, sans que rien change à la construction. La log-vraisemblance s'en déduit. En y remplaçant $\sigma^2$ par l'estimateur qui la maximise, $\hat\sigma^2 = \frac{1}{NK}\sum_i z_i'\Sigma_i^{-1}z_i$, on obtient le critère concentré \begin{equation*} \ell(\theta,\alpha,\omega,\psi) = -\frac{NK}{2}\log\hat\sigma^2(\theta,\alpha,\omega,\psi) - \frac{1}{2}\sum_{i}\log|\Sigma_i| + \sum_{i}\sum_{k\geq2}\log g_{i,k} \end{equation*} que l'on maximise numériquement. #+BEGIN_remarque Le dernier terme est le plus facile à oublier, et l'oublier ruine tout. Nous avons remplacé les données $r$ par des quantités $z$ qui en dérivent, or la densité d'une variable transformée fait intervenir le jacobien de la transformation. Ici la transformation divise les lignes $2$ à $K$ par $g_{i,2},\ldots,g_{i,K}$, si bien que son déterminant vaut $\prod_{k\geq2} 1/g_{i,k}$ et que la log-vraisemblance des données gagne $+\sum_{k\geq2}\log g_{i,k}$. Ce terme dépend de $\theta$, puisque $g_{i,k} = 1-e^{-\theta c_{i,k}\tau}$. Sans lui, l'optimiseur découvre qu'il peut améliorer le critère en déplaçant $\theta$ pour gonfler les $g$, sans rien expliquer de plus, et l'estimation part alors à la borne supérieure. Lorsque les taux sont constants, ce terme et la division par $g$ se compensent exactement, et l'on retrouve la vraisemblance en différences premières usuelle. #+END_remarque \\ #+begin_src python :session panel :exports code :results none def enveloppe(q): """(theta, alpha, omega, psi) à partir de paramètres non contraints.""" return (np.exp(q[0]), 0.01 + 0.98/(1 + np.exp(-q[1])), np.exp(q[2]), 0.95*np.tanh(q[3])) def covariance_ma(psi, K): """Covariance des erreurs, à sigma^2 près : gamma_h = somme_j psi_j psi_{j+h}.""" p = np.concatenate([[1.0], np.atleast_1d(psi)]) V = np.zeros((K, K)) for h in range(len(p)): gh = float(np.dot(p[:len(p)-h], p[h:])) V += np.eye(K, k=h)*gh + (np.eye(K, k=-h)*gh if h else 0) return V def morceaux(theta, alpha, omega, psi, LY, S, C, K, tau, projetee=True): """Vecteur transformé z et sa covariance Sigma, empilée sur les pays, à sigma^2 près.""" r, g = bruts([theta, alpha], LY, S, C, K, tau) if np.any(g <= 1e-8): return None, None, None n = r.shape[1] V = covariance_ma(psi, K) T = np.zeros((n, K, K)) T[:, 0, 0] = 1 for k in range(1, K): T[:, k, k], T[:, k, k-1] = 1/g[k], -1/g[k-1] Sigma = T @ V @ np.transpose(T, (0, 2, 1)) Sigma[:, 0, 0] = omega z = np.empty((n, K)) z[:, 0] = r[0] z[:, 1:] = (r[1:]/g[1:] - r[:-1]/g[:-1]).T W = (np.column_stack([np.ones(n), np.log(S[0]), np.log(C[0])]) if projetee else np.ones((n, 1))) z[:, 0] = z[:, 0] - W @ np.linalg.lstsq(W, z[:, 0], rcond=None)[0] return z, Sigma, g def critere(theta, alpha, omega, psi, LY, S, C, K, tau, projetee=True): """Critère concentré, jacobien compris (au signe près, à minimiser).""" if theta > 5 or omega > 1e6 or omega < 1e-6: return 1e10 z, Sigma, g = morceaux(theta, alpha, omega, psi, LY, S, C, K, tau, projetee) if z is None: return 1e10 try: L = np.linalg.cholesky(Sigma) except np.linalg.LinAlgError: return 1e10 y = np.linalg.solve(L, z[:, :, None])[:, :, 0] quad = np.sum(y*y) if quad <= 0: return 1e10 ldet = 2*np.sum(np.log(np.diagonal(L, axis1=1, axis2=2))) n = z.shape[0] return K*n/2*np.log(quad/(K*n)) + 0.5*ldet + np.sum(np.log(g[1:])) def logv(q, LY, S, C, K, tau, projetee=True): """Le même critère, à partir des paramètres non contraints.""" return critere(*enveloppe(q), LY, S, C, K, tau, projetee) def depart(theta, alpha, omega=1.2, psi=0.25): a = (alpha - 0.01)/0.98 return [np.log(theta), np.log(a/(1 - a)), np.log(omega), np.arctanh(psi/0.95)] def maxv(LY, S, C, K, tau, projetee=True, departs=((0.2, 0.55, 1.2, 0.25), (0.4, 0.4, 1.2, 0.0), (0.1, 0.7, 1.2, 0.4))): """Maximum de vraisemblance, plusieurs points de départ.""" essais = [minimize(logv, depart(*p), args=(LY, S, C, K, tau, projetee), method='Nelder-Mead', options=dict(maxiter=4000, xatol=1e-7, fatol=1e-9)) for p in departs] return enveloppe(min(essais, key=lambda e: e.fun).x) def ecart_type(theta, alpha, omega, psi, LY, S, C, K, tau): """Écart-type des innovations, concentré hors de la vraisemblance.""" z, Sigma, _ = morceaux(theta, alpha, omega, psi, LY, S, C, K, tau) L = np.linalg.cholesky(Sigma) y = np.linalg.solve(L, z[:, :, None])[:, :, 0] return np.sqrt(np.sum(y*y)/(K*z.shape[0])) for projetee in (False, True): print(projetee, np.round(maxv(LY, S, C, K, 5, projetee), 4)) libre = logv(depart(*maxv(LY, S, C, K, 5)), LY, S, C, K, 5) nul = minimize(lambda q: logv([q[0], q[1], q[2], np.arctanh(0.0)], LY, S, C, K, 5), depart(0.2, 0.55)[:3], method='Nelder-Mead').fun print('psi = 0 impose :', round(nul - libre, 1), round(2*(nul - libre), 1)) #+end_src Le choix de la moyenne de cette première observation est décisif, et c'est tout l'intérêt de l'exercice. Si on lui donne une moyenne commune à tous les pays, on obtient $\hat\theta = 0{,}027$ et $\hat\alpha = 0{,}801$, c'est-à-dire moins encore que l'estimation en coupe\nbsp{}: en supposant la condition initiale indépendante des régresseurs, on a réintroduit exactement l'endogénéité que le panel devait éliminer. Si on la projette sur le taux d'épargne initial et sur $\log c_{i,1}$, on obtient $\hat\theta = 0{,}176$, $\hat\alpha = 0{,}508$ et $\hat\omega = 1{,}27$. Toute l'information sur le niveau est dans la condition initiale, et c'est là que loge l'endogénéité. \\ Quant au paramètre de moyenne mobile, il confirme ce que l'agrégation laissait attendre\nbsp{}: $\hat\psi = 0{,}29$, un ordre de grandeur voisin du +0,25 de Working, et l'imposer nul dégrade le critère de 51 points, soit un rapport de vraisemblance de 102 pour un seul degré de liberté — deux fois l'écart de critère, la quantité dont la loi est tabulée par le $\chi^2$ lorsque la contrainte est vraie. La liberté que nous avons donnée aux erreurs n'était donc pas une précaution de style, et la [[#limites][section sur les limites]] établira que l'agrégation en est bien la cause. #+BEGIN_remarque La démonstration de Hsiao, Pesaran et Tahmiscioglu suppose un $\rho$ commun à tous les pays, alors qu'ici $\rho_{i,k}$ varie dans les deux dimensions. Rien ne garantit donc que leur résultat de convergence s'applique tel quel. Nous ne l'invoquerons pas\nbsp{}: la qualité de cet estimateur sera jugée sur le Monte Carlo, comme celle du précédent. On notera au passage qu'un test de recouvrement à faible bruit, décisif pour l'estimateur intra, ne dit rien de celui-ci\nbsp{}: quand $\sigma$ tend vers zéro, les contrastes s'annulent identiquement, quelle que soit la valeur de $\theta$. Cet estimateur tire son information de la structure de covariance des contrastes, non de leur espérance. #+END_remarque * Corriger le biais par simulation :PROPERTIES: :CUSTOM_ID: correction :END: Soumettons donc le nouvel estimateur à la même épreuve. #+begin_src python :session panel :exports code :results none theta_mv, alpha_mv, omega_mv, psi_mv = maxv(LY, S, C, K, 5) sigma_mv = ecart_type(theta_mv, alpha_mv, omega_mv, psi_mv, LY, S, C, K, 5) def efficacites(theta, alpha, LY, S, C, K, tau): """Efficacités initiales : régression sans constante de r sur g, à θ et α donnés.""" r, g = bruts([theta, alpha], LY, S, C, K, tau) return (g*r).sum(axis=0)/(g*g).sum(axis=0) lA0_mv = efficacites(theta_mv, alpha_mv, LY, S, C, K, 5) def esperance_mv(theta, alpha, sigma, psi, LY, S, C, K, tau, B, rng): """Moyenne du MV sur B panels simulés, efficacités déduites de θ et α.""" lA0 = efficacites(theta, alpha, LY, S, C, K, tau) return np.array([maxv(simule(theta, alpha, sigma, psi, lA0, LY, S, C, K, tau, rng), S, C, K, tau, departs=((theta, alpha, omega_mv, psi),))[:2] for _ in range(B)]) des_lien = np.random.default_rng(2015) courbe_mv = {th: esperance_mv(th, 0.508, sigma_mv, psi_mv, LY, S, C, K, 5, 40, des_lien) for th in GRILLE} for th in GRILLE: print(th, courbe_mv[th][:, 0].mean().round(3)) #+end_src #+ATTR_HTML: :class numbers | vrai $\theta$ | 0,10 | 0,176 | 0,30 | 0,60 | | $\mathbb{E}[\hat\theta_{\text{MV}}]$ | 0,093 | 0,129 | 0,182 | 0,484 | Il est biaisé lui aussi, et vers le bas, d'un facteur qui n'a rien de constant\nbsp{}: 7\nbsp{}% en $\theta = 0{,}10$, 27\nbsp{}% à hauteur de l'estimation observée, 39\nbsp{}% en $\theta = 0{,}30$, puis 19\nbsp{}% en $\theta = 0{,}60$. Mais, à la différence du précédent, il est *monotone*, et nettement croissant\nbsp{}: à chaque valeur de son espérance correspond une seule valeur du paramètre. Cela suffit pour le corriger, par une méthode appelée *inférence indirecte*. \\ *L'idée.* Imaginons une balance qui ne pèse pas juste, sans que l'on sache de combien elle se trompe. On peut l'étalonner\nbsp{}: poser dessus des poids dont on connaît la masse, noter ce qu'elle affiche, et tracer la courbe de l'affichage en fonction de la vraie masse. Pour peser ensuite un objet inconnu, on lit l'affichage, on cherche sur la courbe la masse qui le produit, et c'est cette masse que l'on retient. La balance n'a pas été réparée\nbsp{}; on a seulement appris à la lire. Notre estimateur est la balance, et les poids connus sont des panels simulés\nbsp{}: pour une valeur de $\theta$ choisie, on simule le modèle un grand nombre de fois, on applique l'estimateur à chaque panel, et l'on fait la moyenne. On obtient ainsi, point par point, la fonction \begin{equation*} b(\theta) = \mathbb{E}_{\theta}\bigl[\hat\theta\bigr] \end{equation*} que Gouriéroux, Monfort et Renault (1993) appellent la *fonction de lien* — c'est exactement la ligne du tableau ci-dessus. Un estimateur sans biais aurait $b(\theta) = \theta$. Le nôtre non, mais peu importe\nbsp{}: on retient la valeur $\theta^{\star}$ qui aurait produit, en moyenne, l'estimation obtenue sur les vraies données, \begin{equation*} b(\theta^{\star}) = \hat\theta_{\text{obs}} \end{equation*} Graphiquement, on part de $\hat\theta_{\text{obs}} = 0{,}176$ sur l'axe vertical, on rejoint horizontalement la courbe, et on lit l'abscisse. #+begin_src python :session panel :exports none :results none fig, ax = plt.subplots(figsize=(6, 4.5)) vrais = np.array(GRILLE) lien_intra = np.array([courbe[th][:, 0].mean() for th in GRILLE]) lien_mv = np.array([courbe_mv[th][:, 0].mean() for th in GRILLE]) ax.plot(vrais, lien_intra, 'o-', color='grey', label='estimateur intra') ax.plot(vrais, lien_mv, 'o-', color='r', label='maximum de vraisemblance') ax.plot(vrais, vrais, 'k--', linewidth=1, label='estimateur sans biais') lu = np.interp(theta_mv, lien_mv, vrais) ax.plot([0.08, lu, lu], [theta_mv, theta_mv, 0], 'r:', linewidth=1.2) ax.annotate(r'$\hat\theta_{\mathrm{obs}}$', (0.08, theta_mv), xytext=(2, 4), textcoords='offset points', fontsize=9, color='r') ax.annotate(r'$\theta^{\star}$', (lu, 0), xytext=(4, 4), textcoords='offset points', fontsize=9, color='r') ax.set_xlim(0.08, 0.62) ax.set_ylim(0, 0.72) ax.set_xlabel(r'vraie valeur de $\theta$') ax.set_ylabel(r"espérance de l'estimateur, $b(\theta)$") ax.legend(fontsize=8, loc='upper left') fig.tight_layout() fig.savefig("panel-inversion.svg", transparent=True) print(round(lu, 3)) #+end_src #+CAPTION: *Les deux fonctions de lien, et la lecture inverse.* On part de l'estimation obtenue sur les données, on rejoint la courbe rouge, on lit la valeur du paramètre qui la produit en moyenne. La même lecture sur la courbe grise est impossible : presque plate, elle ne passe que de 0,452 à 0,487 quand $\theta$ triple. #+LABEL: fig:inversion [[file:panel-inversion.svg]] La lecture graphique, à $\alpha$ fixé, donne déjà $\theta^{\star}\approx 0{,}29$. \\ *Ce qu'il faut, et ce qu'il ne faut pas.* La méthode demande trois choses. Il faut d'abord pouvoir simuler le modèle pour n'importe quelle valeur des paramètres\nbsp{}; c'est le cas, l'équation du modèle fournit son propre générateur. Il faut ensuite appliquer aux données simulées *exactement* le même estimateur qu'aux données réelles, avec les mêmes revenus de 1963, les mêmes taux et le même découpage du temps\nbsp{}: c'est ce qui garantit que le biais, quelle qu'en soit la source, frappe les deux de la même façon, et se neutralise lorsqu'on les compare. Il faut enfin que la fonction de lien soit strictement monotone, et suffisamment pentue\nbsp{}: là où elle est plate, une petite variation de l'estimation observée correspond à une grande variation du paramètre, et l'inversion n'a plus de sens. En revanche — et c'est ce qui surprend d'abord — il n'est *pas* nécessaire que l'estimateur que l'on inverse, qu'on appelle l'estimateur *auxiliaire*, soit sans biais, ni même convergent. Il n'a pas à être bon, seulement à être informatif. Voilà pourquoi l'estimateur intra, dont la courbe est plate, ne peut pas servir, alors que le maximum de vraisemblance, pourtant biaisé lui aussi, le peut. Gouriéroux, Phillips et Yu (2010) ont appliqué la méthode au panel dynamique en inversant l'estimateur intra\nbsp{}; dans leur modèle linéaire il est pentu, dans le nôtre il ne l'est pas. \\ *Avec deux paramètres.* La correction ne peut pas se faire à $\alpha$ fixé, car l'estimateur de $\alpha$ est biaisé lui aussi, et son biais dépend de $\theta$. La fonction de lien associe donc un couple d'espérances $(\mathbb{E}[\hat\theta],\mathbb{E}[\hat\alpha])$ à chaque couple $(\theta,\alpha)$, et il faut résoudre un système de deux équations à deux inconnues. On ne peut plus le lire sur un graphique, mais on le résout par itérations. Partant de $p_0 = (\hat\theta_{\text{obs}},\hat\alpha_{\text{obs}})$, on calcule \begin{equation*} p_{j+1} = p_j + \lambda\bigl(\hat p_{\text{obs}} - \hat b(p_j)\bigr) \end{equation*} où $\hat b(p_j)$ est la moyenne des estimations sur $B$ panels simulés avec les paramètres $p_j$. Une précision s'impose sur ce que contient $p_j$\nbsp{}: l'itération ne porte que sur $\theta$ et $\alpha$, mais les 110 efficacités n'en sont pas moins remises à jour à chaque pas. Elles ne sont pas des inconnues supplémentaires\nbsp{}: à $\theta$ et $\alpha$ donnés, elles se lisent dans les données par la régression sans constante de $r_{i,k}$ sur $g_{i,k}$, celle-là même qui sert à concentrer la vraisemblance. Les figer au point non corrigé reviendrait à simuler un monde dont les efficacités contrediraient sa propre vitesse de convergence, et nous verrons plus bas que ce n'est pas une subtilité d'écriture. Restent fixés l'écart-type des perturbations, la variance de la première observation et le paramètre de moyenne mobile\nbsp{}; pour $\sigma$, qui a lui aussi une forme explicite une fois $\omega$ et $\psi$ donnés, le remettre à jour ne déplace pas $\theta^{\star}$ au-delà du millième. Une précision encore\nbsp{}: les courbes de biais tracées plus haut, calculées avant que le maximum de vraisemblance ne soit construit, utilisent les efficacités de l'estimateur intra. La première itération est la correction de biais que chacun ferait spontanément\nbsp{}: mesurer le biais au point estimé et le retrancher. Si le biais était le même pour toutes les valeurs du paramètre, elle suffirait. Il ne l'est pas — de 7 à 39\nbsp{}% selon la valeur visée — et les itérations suivantes corrigent la correction, jusqu'à ce que les paramètres simulés reproduisent l'estimation observée. Le coefficient $\lambda = 0{,}9$ amortit les oscillations dues au bruit de simulation. #+begin_src python :session panel :exports code :results none def correction(vise, sigma, psi, LY, S, C, K, tau, rng, iterations=20, B=40, pas=0.9, garde=8): """θ et α tels que l'estimateur simulé reproduise l'estimation observée.""" p, trace = np.array(vise, dtype=float), [] for _ in range(iterations): out = esperance_mv(p[0], p[1], sigma, psi, LY, S, C, K, tau, B, rng).mean(axis=0) p = p + pas*(np.array(vise) - out) p = np.array([min(max(p[0], 1e-3), 3), min(max(p[1], 0.05), 0.95)]) trace.append(p.copy()) return np.mean(trace[-garde:], axis=0) # moyenne de fin de course, contre le bruit etoile = correction([theta_mv, alpha_mv], sigma_mv, psi_mv, LY, S, C, K, 5, np.random.default_rng(2016)) lA0_etoile = efficacites(*etoile, LY, S, C, K, 5) print(etoile.round(3), np.median(etoile[0]*cbar).round(4), (etoile[0] - (1 - etoile[1])).round(3)) #+end_src La correction donne $\theta^{\star} = 0{,}36$ et $\alpha^{\star} = 0{,}32$, soit une vitesse de convergence médiane de 2,9\nbsp{}% par an\nbsp{}: le paramètre corrigé vaut le double de l'estimation brute. Remarquons que la correction déplace aussi la part du capital, de 0,51 à 0,32\nbsp{}: c'est pourquoi il ne suffisait pas de lire la courbe de biais à $\alpha$ fixé, et pourquoi le point fixe porte sur les deux paramètres à la fois. \\ *Quelle précision\nbsp{}?* Puisque $p^{\star} = b^{-1}(\hat p_{\text{obs}})$, l'incertitude sur $\hat p_{\text{obs}}$ se transmet à $p^{\star}$ à travers la réciproque de la fonction de lien. Au premier ordre, si $J$ désigne la matrice des dérivées de $b$ et $V$ la matrice de variance de l'estimateur auxiliaire, toutes deux évaluées en $p^{\star}$, \begin{equation*} \mathbb{V}\bigl(p^{\star}\bigr) \simeq \Bigl(1+\frac{1}{B}\Bigr)J^{-1}\,V\,{J^{-1}}' \end{equation*} C'est la méthode delta, c'est-à-dire la propagation d'une incertitude à travers une fonction que l'on remplace par sa tangente, ici celle de la réciproque de $b$. Le facteur $1+1/B$ tient compte de ce que $b$ n'est connue qu'à travers $B$ simulations, et c'est le $B$ de la correction, quarante, qui entre ici. La formule dit aussi, en une ligne, pourquoi la pente compte\nbsp{}: diviser par une dérivée proche de zéro, celle de l'estimateur intra, rendrait la variance infinie. Les deux ingrédients s'obtiennent encore par simulation, $V$ en répétant l'estimation en $p^{\star}$, $J$ par différences finies, en réutilisant les mêmes tirages de part et d'autre du point pour que le bruit de simulation ne se confonde pas avec la pente. #+begin_src python :session panel :exports code :results none graine, h, Bp = 1, 0.03, 60 base = esperance_mv(*etoile, sigma_mv, psi_mv, LY, S, C, K, 5, Bp, np.random.default_rng(graine)) J = np.empty((2, 2)) for j in range(2): pas_j = h*np.eye(2)[j] haut, bas = (esperance_mv(*(etoile + signe*pas_j), sigma_mv, psi_mv, LY, S, C, K, 5, Bp, np.random.default_rng(graine)).mean(axis=0) for signe in (1, -1)) J[:, j] = (haut - bas)/(2*h) Ji = np.linalg.inv(J) et_etoile = np.sqrt(np.diag((1 + 1/40)*Ji @ np.cov(base.T) @ Ji.T)) print(J.round(2), et_etoile.round(3), (np.median(cbar)*(etoile[0] + np.array([-1.645, 1.645])*et_etoile[0])).round(4)) #+end_src Et c'est ici que la correction montre son prix. La pente de $\mathbb{E}[\hat\theta]$ par rapport à $\theta$ ne vaut que 0,53\nbsp{}: doubler la vraie vitesse ne déplace l'estimateur auxiliaire que de moitié autant, si bien que la lecture inverse *multiplie par deux* toute incertitude. L'écart-type de $\theta^{\star}$ est de 0,167, celui de $\alpha^{\star}$ de 0,134, et l'intervalle de confiance à 90\nbsp{}% de la vitesse médiane va de 0,7 à 5,2\nbsp{}% par an. C'est large, et il faut le dire ainsi\nbsp{}: les données situent la vitesse au-dessus de ce que l'estimation brute affiche, mais elles ne la situent pas précisément. \\ Un mot sur ce que coûterait la facilité. Si l'on figeait les efficacités au point non corrigé, au lieu de les recalculer, la pente afficherait 1,17 et l'écart-type 0,038 — une précision quatre fois meilleure, et entièrement illusoire. Elle vient de ce que le gel empêche les efficacités de répondre à $\theta$\nbsp{}: la fonction de lien paraît alors trois fois plus raide qu'elle ne l'est, et son inversion bien mieux déterminée. \\ La contrainte du modèle de Solow, $\theta = 1-\alpha$, se teste alors sans invoquer aucune loi asymptotique, en simulant sous la contrainte. #+begin_src python :session panel :exports code :results none sous_h0 = esperance_mv(1 - etoile[1], etoile[1], sigma_mv, psi_mv, LY, S, C, K, 5, 150, np.random.default_rng(2017)) ecart0 = sous_h0[:, 0] - (1 - sous_h0[:, 1]) ecart = theta_mv - (1 - alpha_mv) print(ecart0.mean().round(3), np.percentile(ecart0, [5, 95]).round(3), round(ecart, 3), (ecart0 <= ecart).mean()) #+end_src Sous l'hypothèse nulle, l'écart $\hat\theta - (1-\hat\alpha)$ n'est pas centré — il vaut −0,109 en moyenne, puisque les deux estimations qui le composent sont elles-mêmes biaisées — et tient dans l'intervalle [−0,224\nbsp{}; +0,016] dans 90\nbsp{}% des échantillons. C'est précisément pourquoi le test se simule au lieu de se lire dans une table. L'écart observé vaut −0,316, et aucun échantillon simulé ne l'atteint\nbsp{}: le modèle de Solow complet est rejeté. Cette fois, le verdict ne repose pas sur une approximation du $\chi^2$ dont la section sur le biais a montré qu'elle ne vaudrait rien ici. \\ Une question plus élémentaire précède pourtant celle-là, et la section [[#panel-abc][Ce qu'un panel ajoute à une coupe]] l'a rendue pressante\nbsp{}: y a-t-il seulement convergence\nbsp{}? Le revenu par tête n'est pas stationnaire, et nous avons vu qu'une marche aléatoire, qui ne converge vers rien, reproduit exactement la pente intra des données. Soumettons donc le maximum de vraisemblance à la même épreuve. Simulons des mondes sans convergence — chaque revenu suit une marche aléatoire, avec une dérive commune à tous les pays ou propre à chacun, et, dans le second cas, les erreurs corrélées de notre modèle — et estimons le modèle sur chacun d'eux. #+begin_src python :session panel :exports code :results none sans_convergence = np.random.default_rng(11) choc_rw = (dy - dy.mean(axis=0)).std() def monde_sans_convergence(propre, psi): """Marches aléatoires calibrées sur les données : dérive commune ou propre à chaque pays.""" e = sans_convergence.normal(0, choc_rw/np.sqrt(1 + psi*psi), (K + 1, N)) pas_rw = (dy.mean(axis=0) if propre else dy.mean()) + e[1:] + psi*e[:-1] return LY[0] + np.vstack([np.zeros(N), np.cumsum(pas_rw, axis=0)]) for propre, psi in ((False, 0.0), (True, psi_mv)): R = [] for _ in range(300): Y = monde_sans_convergence(propre, psi) R.append([maxv(Y, S, C, K, 5)[0], intra(Y, S, C, K, 5).x[0]]) R = np.array(R) print(propre, R[:, 0].mean().round(3), np.percentile(R[:, 0], 99).round(3), (R[:, 0] >= theta_mv).mean(), R[:, 1].mean().round(3), np.median(R[:, 1].mean()*cbar).round(4)) #+end_src Le verdict est sans ambiguïté. Dans un monde sans convergence, le maximum de vraisemblance donne en moyenne $\hat\theta = 0{,}025$ avec une dérive commune, 0,058 avec des dérives propres à chaque pays, et ne dépasse respectivement 0,10 et 0,13 que dans 1\nbsp{}% des cas. Sur 600 panels simulés, aucun n'atteint la valeur obtenue sur les données, 0,176. Contrairement à l'estimateur intra, le maximum de vraisemblance distingue donc nos données d'un monde sans convergence\nbsp{}; et l'estimateur intra, appliqué aux mêmes panels, y trouve en moyenne 0,63 et 0,23, c'est-à-dire des vitesses médianes de 5,2 et 1,9\nbsp{}% par an là où il n'y en a aucune — des chiffres que personne n'hésiterait à publier, ce qui montre assez qu'une valeur plausible n'est pas une preuve. La convergence que mesure cette note n'est pas un artefact de la non-stationnarité du revenu. Ce test a toutefois une portée limitée\nbsp{}: il écarte la marche aléatoire, pas toutes les formes de non-convergence, et en particulier pas des tendances déterministes propres à chaque pays, que nous n'avons pas simulées. #+begin_src python :session panel :exports code :results none reconcilie = esperance(theta_mv, alpha_mv, sigma_mv, psi_mv, lA0_mv, LY, S, C, K, 5, 400, np.random.default_rng(2019))[:, 0] print(reconcilie.mean().round(3), np.percentile(reconcilie, [5, 95]).round(3), theta5.round(3), (reconcilie >= theta5).mean().round(3)) #+end_src #+BEGIN_remarque Voici l'épreuve la plus exigeante que l'on puisse faire subir à l'ensemble. Simulons le modèle avec les paramètres estimés par le maximum de vraisemblance, et appliquons aux données ainsi produites l'estimateur intra\nbsp{}: sur 400 réplications, il donne en moyenne 0,387, et tombe entre 0,314 et 0,471 dans 90\nbsp{}% des cas. Or l'estimation intra obtenue sur les vraies données vaut 0,478 — au-dessus, de peu, mais au-dessus\nbsp{}: le modèle produit une valeur au moins aussi grande dans 3,5\nbsp{}% des échantillons. Il faut lire ce résultat dans les deux sens. D'un côté, l'essentiel de l'écart entre 0,18 et 0,48 est bien expliqué\nbsp{}: une même valeur de $\theta$, vue par deux estimateurs, produit deux nombres que tout semblait opposer — et la section [[#retour][Retour à la coupe transversale]] fera la même démonstration pour l'estimation en coupe. De l'autre, la réconciliation n'est pas parfaite. Le modèle peine à engendrer un estimateur intra aussi élevé que celui des données\nbsp{}; il reste donc, dans le comportement de l'estimateur biaisé, quelque chose que notre spécification ne reproduit pas tout à fait. #+END_remarque * Mesure-t-on vraiment l'efficacité du travail ? :PROPERTIES: :CUSTOM_ID: efficacite :END: Le panel livre une estimation de $\log A_i(0)$ pour chacun des 110 pays\nbsp{}: c'est en principe un acquis précieux, puisqu'il s'agit du niveau d'efficacité du travail en 1963, que la coupe transversale ne pouvait qu'absorber dans une perturbation. Encore faut-il vérifier que ce nombre mesure ce que son nom promet. \\ La comparaison demande une précaution. Notre $A$ augmente le travail, au sens de Harrod\nbsp{}: $Y = K^{\alpha}(AL)^{1-\alpha}$. Les Penn World Tables publient une productivité globale des facteurs au sens de Hicks, `ctfp`, définie par $Y = A_H K^{\alpha}(hL)^{1-\alpha}$ où $h$ est le capital humain par tête. Ce ne sont pas le même objet. En égalant les deux écritures, \begin{equation*} A^{1-\alpha} = A_H\,h^{1-\alpha} \qquad\Longleftrightarrow\qquad \log A = \frac{\log A_H}{1-\alpha} + \log h \end{equation*} si bien que la contrepartie empirique de notre estimation n'est pas $\log$ `ctfp`, mais $\log$ `ctfp` divisé par $1-\alpha$, plus le logarithme du capital humain. Avec $\hat\alpha = 0{,}51$ — la valeur non corrigée, celle à laquelle les $\log\hat A_i(0)$ ont été estimés — le diviseur vaut environ 2\nbsp{}: les écarts de productivité globale se trouvent doublés lorsqu'on les exprime en efficacité du travail. #+begin_src python :session panel :exports code :results none mesures = (pwt[pwt['year'] == DEBUT].set_index('countrycode') .reindex(codes)[['ctfp', 'hc']]) bon = ((mesures['ctfp'] > 0) & mesures['hc'].notna()).to_numpy() cible = np.log(mesures['ctfp'].to_numpy(float)[bon])/(1 - alpha_mv) \ + np.log(mesures['hc'].to_numpy(float)[bon]) notre = lA0_mv[bon] rangs = lambda v: np.argsort(np.argsort(v)) hc = np.log(mesures['hc'].to_numpy(float)[bon]) print(np.corrcoef(notre, hc)[0, 1].round(3)) # lien avec le capital humain print(bon.sum(), np.corrcoef(notre, cible)[0, 1].round(3), np.corrcoef(rangs(notre), rangs(cible))[0, 1].round(3), np.polyfit(cible, notre, 1)[0].round(3), cible.std().round(3), notre.std().round(3)) u, iu = list(np.array(codes)[bon]), list(np.array(codes)[bon]).index('USA') for pays in ('JPN', 'IND'): i = u.index(pays) print(pays, np.exp(cible[i] - cible[iu]).round(2), np.exp(notre[i] - notre[iu]).round(2)) def norme1(x, y): """Droite qui minimise la somme des écarts absolus, moins sensible aux points isolés.""" return minimize(lambda p: np.abs(y - p[0] - p[1]*x).sum(), np.polyfit(x, y, 1)[::-1], method='Nelder-Mead', options=dict(xatol=1e-8, fatol=1e-10)).x mco, l1 = np.polyfit(cible, notre, 1)[0], norme1(cible, notre)[1] jeté = np.array([np.polyfit(np.delete(cible, i), np.delete(notre, i), 1)[0] for i in range(len(cible))]) # pente en retirant un pays dedans = lambda v: (v > np.percentile(v, 5)) & (v < np.percentile(v, 95)) centre = dedans(cible) & dedans(notre) # sans les 5 % extrêmes de chaque axe print(round(mco, 3), round(l1, 3), round(jeté.min(), 3), round(jeté.max(), 3), u[int(np.argmax(np.abs(jeté - mco)))]) haut = np.argsort(notre)[::-1][:4] # les quatre plus hautes estimations garde = np.ones(len(cible), bool) garde[haut] = False print([u[i] for i in haut], int(centre[haut].sum()), round(np.polyfit(cible[garde], notre[garde], 1)[0], 3), round(norme1(cible[garde], notre[garde])[1], 3)) print(centre.sum(), round(np.polyfit(cible[centre], notre[centre], 1)[0], 3), round(norme1(cible[centre], notre[centre])[1], 3), round(stats.spearmanr(cible[centre], notre[centre]).statistic, 3), round(stats.spearmanr(cible[centre], notre[centre]).pvalue, 3), round(stats.spearmanr(cible, notre).pvalue, 3)) #+end_src Le résultat est un avertissement. Sur les 62 pays pour lesquels les deux mesures existent, la corrélation entre notre $\log\hat A_i(0)$ et sa contrepartie est de *zéro*, et la pente de la régression vaut 0,00 là où la théorie en exige une de 1. Les dispersions, en revanche, sont du même ordre — 1,27 pour notre estimation sur ces 62 pays, 0,99 pour la contrepartie — de sorte que le problème n'est pas d'échelle mais de classement. Les cas particuliers le disent mieux qu'un coefficient\nbsp{}: les PWT situent le Japon de 1963 à 27\nbsp{}% de l'efficacité américaine, notre estimation à 119\nbsp{}%\nbsp{}; l'Inde à 5\nbsp{}% contre 20\nbsp{}%. \\ Une pente nulle sur 62 points doit toutefois être auscultée avant d'être crue, car quelques pays suffisent à la produire. Trois vérifications. D'abord, la pente des moindres carrés est *fragile*\nbsp{}: retirer un seul pays la déplace entre −0,08 et +0,15, et c'est Malte, que les PWT placent très bas et notre estimation très haut, qui pèse le plus. Ensuite, on peut ajuster la droite en minimisant la somme des écarts *absolus* plutôt que celle de leurs carrés — la norme L1, qui ne laisse pas un point éloigné peser par son carré\nbsp{}: la pente monte alors à 0,17. Enfin, en écartant les 5\nbsp{}% extrêmes de chaque axe, il reste 48 pays, et la pente des moindres carrés vaut 0,34, celle de la norme L1 0,43, avec une corrélation des rangs de 0,32 qui, cette fois, n'est pas fortuite (probabilité critique de 0,03). Sur l'ensemble des 62 pays, la même corrélation des rangs ne vaut que 0,13 et pourrait très bien venir du hasard (0,31). \\ Ces points écartés ne sont pas quelconques, et il vaut la peine de les nommer. Les quatre plus hautes estimations d'efficacité initiale sont celles de Malte, de la Corée du Sud, de l'Irlande et de Taïwan, toutes quatre très basses selon les PWT\nbsp{}; ce sont elles, en haut à gauche du graphique, qui aplatissent la droite, et les retirer suffit à porter la pente des moindres carrés de 0,00 à 0,29, celle de la norme L1 de 0,17 à 0,41. Ce sont aussi, et ce n'est pas un hasard, les quatre grandes réussites de rattrapage de la période. Le modèle n'a qu'un seul levier pour rendre compte d'une croissance durablement supérieure à $x$\nbsp{}: relever le niveau du sentier vers lequel le pays se dirige, c'est-à-dire $\log A_i(0)$. Notre estimation mesure donc une *cible* là où les PWT mesurent un *niveau atteint* en 1963, et c'est précisément pour les économies de rattrapage que les deux notions s'écartent le plus. #+begin_src python :session panel :exports none :results none fig, ax = plt.subplots(figsize=(6, 4.5)) u = list(np.array(codes)[bon]) iu = u.index('USA') px, py = cible - cible[iu], notre - notre[iu] ax.scatter(px[centre], py[centre], s=12, color='grey') ax.scatter(px[~centre], py[~centre], s=18, marker='x', color='grey', label='pays écartés du petit échantillon') ecarts = {'JPN': (5, 4), 'FRA': (6, -9), 'GBR': (-22, 3), 'IND': (5, 3), 'MLT': (7, -2), 'KOR': (7, -2), 'IRL': (7, -2), 'TWN': (7, -2), 'VEN': (6, 2)} for pays, ecart in ecarts.items(): if pays in u: i = u.index(pays) ax.annotate(pays, (px[i], py[i]), xytext=ecart, textcoords='offset points', fontsize=8, color='k') lim = np.array([px.min() - 0.2, px.max() + 0.2]) ax.plot(lim, lim, 'k:', linewidth=1, label='ce que la théorie prédit') for masque, trait, taille in ((slice(None), '-', '62 pays'), (centre, '--', '48 pays')): x, y = px[masque], py[masque] ax.plot(lim, np.polyval(np.polyfit(x, y, 1), lim), trait, color='tab:blue', linewidth=1, label=f'moindres carrés, {taille}') c = norme1(x, y) ax.plot(lim, c[0] + c[1]*lim, trait, color='tab:red', linewidth=1, label=f'norme L1, {taille}') ax.set_xlabel(r'efficacité impliquée par les PWT (USA = 0)') ax.set_ylabel(r'notre estimation, $\log \hat A_i(0)$ (USA = 0)') ax.legend(fontsize=7) fig.tight_layout() fig.savefig("panel-efficacite.svg", transparent=True) #+end_src #+CAPTION: *L'efficacité initiale estimée par le panel, confrontée à celle qu'impliquent les Penn World Tables.* Les points devraient s'aligner sur la diagonale. En trait plein, les ajustements sur les 62 pays ; en tirets, ceux qui écartent les 5 % extrêmes de chaque axe, marqués d'une croix. Le bleu minimise les carrés des écarts, le rouge leur valeur absolue. #+LABEL: fig:efficacite [[file:panel-efficacite.svg]] #+BEGIN_remarque Que retenir de ces trois lectures\nbsp{}? Non pas qu'il n'y a aucune relation — l'ajustement robuste en trouve une —, mais qu'elle est faible, très imprécise, et trois fois trop plate pour la théorie, qui exige une pente de un. Il serait donc malhonnête de présenter cette estimation du niveau comme un acquis du panel. Ce que le panel identifie bien, c'est la *dynamique*\nbsp{}: la vitesse à laquelle les pays rejoignent leur sentier. Le *niveau* $\log A_i(0)$, lui, est un paramètre de nuisance qui absorbe tout ce que le modèle ne sait pas expliquer — l'absence de capital humain, les erreurs de mesure du revenu initial, l'écart à l'état stationnaire en 1963 — et rien ne garantit que la somme de ces résidus ressemble à une efficacité. Le signe encourageant est qu'il corrèle à +0,44 avec le capital humain, ce qui est cohérent puisque notre modèle n'en contient pas et que $A$ doit l'absorber. Cela tempère aussi le verdict que la section suivante tirera de ces mêmes inconnues\nbsp{}: leur corrélation avec le taux d'épargne réfute bien l'hypothèse d'orthogonalité, mais elle ne dit rien de la valeur de $\log A_i(0)$ comme mesure. #+END_remarque * Retour à la coupe transversale :PROPERTIES: :CUSTOM_ID: retour :END: Nous pouvons maintenant revenir aux questions que [[https://stephane-adjemian.fr/posts/vitesses-de-convergence-propres-a-chaque-pays/][la note en coupe]] avait dû laisser ouvertes. Elles étaient au nombre de trois, et chacune tenait à ce qu'une seule observation par pays ne permet pas de séparer l'efficacité du travail du reste. 1. *Le revenu initial contient-il l'efficacité du travail\nbsp{}?* Si oui, l'estimateur en coupe est biaisé. Les simulations de la note précédente encadraient l'estimation observée, $\hat\theta = 0{,}102$, entre deux mondes extrêmes — revenu initial sans efficacité (toutes les estimations au-dessus de 0,17), revenu initial la contenant tout entière (toutes au-dessous de 0,09) — et rien ne permettait de dire où se situaient les données. 2. *L'efficacité du travail est-elle indépendante du taux d'épargne\nbsp{}?* Mankiw, Romer et Weil le supposent, et la coupe ne pouvait que l'admettre. 3. *Le rejet du modèle de Solow est-il réel\nbsp{}?* Il pouvait n'être que l'effet de l'endogénéité, et il reposait de surcroît sur une loi du $\chi^2$ dont la taille réelle était de 8\nbsp{}% et non de 5\nbsp{}%. *Commençons par la deuxième*, la plus facile, car le panel y répond sans rien simuler et sans rien estimer de plus. L'efficacité initiale de chaque pays est désormais un nombre\nbsp{}: il suffit de regarder s'il va de pair avec le taux d'épargne. L'hypothèse de Mankiw, Romer et Weil affirme qu'il n'y a aucune relation, de sorte qu'une corrélation non nulle suffit à la réfuter. #+begin_src python :session panel :exports code :results none ls, lc = np.log(S.mean(axis=0)), np.log(cbar) tirage = np.random.default_rng(8) def liaison(v, x): """Corrélations de Pearson et de Spearman, et intervalle par rééchantillonnage.""" r, p = stats.pearsonr(v, x) bs = [np.corrcoef(v[i], x[i])[0, 1] for i in (tirage.integers(0, N, N) for _ in range(2000))] return round(r, 3), round(p, 5), round(stats.spearmanr(v, x).statistic, 3), \ np.percentile(bs, [2.5, 97.5]).round(3) def partielle(v, x, z): """Corrélation de v et x une fois retirée de chacun la part expliquée par z.""" Z = np.column_stack([np.ones(N), z]) reste = lambda u: u - Z @ np.linalg.lstsq(Z, u, rcond=None)[0] return round(np.corrcoef(reste(v), reste(x))[0, 1], 3) print(liaison(lA0_mv, ls), liaison(lA0_mv, lc)) print(liaison(lA0, ls), liaison(lA0, lc)) # avec les inconnues de l'estimateur intra print(round(np.corrcoef(ls, lc)[0, 1], 3), # les deux variables sont elles-mêmes liées partielle(lA0_mv, ls, lc), partielle(lA0_mv, lc, ls)) #+end_src Il y en a une, et elle est solide. #+ATTR_HTML: :class numbers | | corrélation | des rangs | partielle | |-+-+-+-| | $\log\bar s_i$ | +0,35 | +0,42 | +0,29 | | $\log\bar c_i$ | −0,21 | −0,29 | −0,08 | La première colonne est la corrélation usuelle, celle de Pearson. La deuxième porte sur les rangs plutôt que sur les valeurs, et son intérêt est d'être insensible à quelques pays extrêmes\nbsp{}: elle confirme la première. La troisième est la *corrélation partielle*, celle qui subsiste entre deux variables une fois retirée de chacune la part que la troisième explique\nbsp{}: elle est utile ici parce que le taux d'épargne et $\bar c_i$ sont eux-mêmes corrélés, à −0,40, de sorte que chacun pourrait n'apparaître que par l'entremise de l'autre. \\ La lecture est nette. Le taux d'épargne va bien avec l'efficacité, et cette liaison ne passe pas par $\bar c_i$, puisqu'elle survit à sa neutralisation. La corrélation apparente de $\bar c_i$, en revanche, s'évanouit dès que l'on tient l'épargne fixe. La liaison avec l'épargne, elle, est franche\nbsp{}: sa probabilité critique vaut 0,0002 et son intervalle de confiance à 95\nbsp{}%, obtenu par rééchantillonnage des pays, va de 0,16 à 0,52. L'hypothèse d'orthogonalité est donc fausse dans ces données, et l'estimation en coupe était bien biaisée. #+begin_src python :session panel :exports none :results none fig, ax = plt.subplots(figsize=(6, 4.5)) ax.scatter(100*S.mean(axis=0), lA0_mv, s=12, color='grey') for pays in ('FRA', 'USA', 'JPN', 'IND', 'NGA', 'CHN'): if pays in list(codes): i = list(codes).index(pays) ax.scatter(100*S.mean(axis=0)[i], lA0_mv[i], s=16, color='r') ax.annotate(pays, (100*S.mean(axis=0)[i], lA0_mv[i]), xytext=(5, 3), textcoords='offset points', fontsize=8, color='r') ajust = np.polyfit(ls, lA0_mv, 1) grille = np.linspace(ls.min(), ls.max(), 20) ax.plot(100*np.exp(grille), np.polyval(ajust, grille), 'k--', linewidth=1, label='ajustement linéaire') ax.set_xscale('log') ax.set_xlabel("taux d'épargne moyen du pays (%, échelle logarithmique)") ax.set_ylabel(r'efficacité initiale estimée, $\log \hat A_i(0)$') ax.legend(fontsize=8) fig.tight_layout() fig.savefig("panel-epargne.svg", transparent=True) #+end_src #+CAPTION: *L'efficacité initiale estimée et le taux d'épargne moyen, pays par pays.* Sous l'hypothèse de Mankiw, Romer et Weil, ce nuage devrait être sans pente. #+LABEL: fig:epargne [[file:panel-epargne.svg]] On aurait pu régresser l'efficacité sur les deux taux, mais cela n'aurait rien appris de plus\nbsp{}: la théorie ne prédit pas une pente particulière, elle prédit l'absence de relation, et la valeur d'un coefficient de régression n'aurait ici aucune interprétation structurelle. La corrélation partielle donne le même renseignement que la régression sans en suggérer une. #+BEGIN_remarque Deux réserves sur ces chiffres. D'abord, l'efficacité initiale n'est pas une donnée mais une estimation\nbsp{}: les probabilités critiques et l'intervalle ci-dessus ignorent cette incertitude-là. Ensuite, ils dépendent un peu de l'estimateur retenu\nbsp{}: les inconnues de l'estimateur intra donnent des corrélations plus fortes, +0,47 avec le taux d'épargne et −0,40 avec $\bar c_i$. La conclusion, elle, ne change pas. #+END_remarque \\ *Les deux autres questions* demandent une dernière expérience. #+begin_src python :session panel :exports code :results none def coupe(Y): """Estimateur M2 de la note en coupe, sur les revenus de 1963 et de 2023.""" t = FIN - DEBUT gy, ly0, sb = (Y[-1] - Y[0])/t, Y[0], S.mean(axis=0) def moins_logv(p): theta, a, alpha = p w = (1 - np.exp(-theta*cbar*t))/t u = (gy - X - w*(a + alpha/(1 - alpha)*(np.log(sb) - np.log(cbar)) - ly0))/w return N/2*np.log(u @ u/N) + np.sum(np.log(w)) essais = [minimize(moins_logv, p0, method='L-BFGS-B', bounds=[(1e-4, 3), (-50, 50), (0.01, 0.99)]) for p0 in ([0.1, 7, 0.7], [0.3, 6, 0.5], [0.05, 8, 0.8], [1, 8, 0.3])] return min(essais, key=lambda e: e.fun).x[[0, 2]] des = np.random.default_rng(3) en_coupe = coupe(LY) monde = np.array([coupe(simule(*etoile, sigma_mv, psi_mv, lA0_etoile, LY, S, C, K, 5, des)) for _ in range(300)]) commune = np.full(N, lA0_etoile.mean()) # aucune différence permanente entre pays sans_lien = np.array([coupe(simule(*etoile, sigma_mv, psi_mv, commune, LY, S, C, K, 5, des)) for _ in range(300)]) print(en_coupe.round(3), np.corrcoef(lA0_mv, LY[0])[0, 1].round(2)) for nom, tirages in (('monde', monde), ('sans lien', sans_lien)): print(nom, tirages.mean(axis=0).round(3), np.percentile(tirages, [5, 95], axis=0).round(3), np.median(tirages[:, 0]).round(3), tirages[:, 0].std().round(3), (tirages[:, 0] <= en_coupe[0]).mean().round(3)) print(np.median(np.outer([en_coupe[0], theta5, theta_mv, etoile[0]], cbar), axis=1).round(4)) #+end_src Appliqué aux années 1963 et 2023 du panel, l'estimateur de la note en coupe donne $\hat\theta = 0{,}105$ et $\hat\alpha = 0{,}708$, à peu près ce qu'il donnait sur 1960-2023. Simulons maintenant le panel avec les paramètres corrigés, $\theta^{\star} = 0{,}36$ et $\alpha^{\star} = 0{,}32$, les revenus observés en 1963 et les efficacités qu'ils impliquent, puis appliquons à chaque panel simulé ce même estimateur en coupe, sans rien y changer. Il donne en moyenne $\hat\theta = 0{,}089$ et $\hat\alpha = 0{,}74$, avec un écart-type de 0,022 pour le premier, dont l'intervalle central à 90\nbsp{}% va de 0,057 à 0,126. L'estimation obtenue sur les vraies données, 0,105, tombe à l'intérieur de cet intervalle, au-dessus de la médiane simulée, 0,086. **L'estimation en coupe est donc exactement ce qu'un monde où $\theta$ vaut 0,36 produit lorsqu'on le regarde avec une seule observation par pays.** \\ Reste à savoir d'où vient l'écart entre 0,36 et 0,09. Refaisons la même expérience dans un monde qui ne diffère du précédent que sur un point\nbsp{}: tous les pays y ont la même efficacité initiale, si bien que le revenu de 1963 ne peut plus en contenir la moindre part. L'estimateur en coupe retrouve alors $\hat\theta = 0{,}43$ et $\hat\alpha = 0{,}31$ en moyenne, cette fois du bon ordre de grandeur, et même un peu au-dessus de la vraie valeur. Tout le biais de la coupe vient donc des différences permanentes d'efficacité entre pays et de leur présence dans le revenu initial, qui est bien réelle dans les données\nbsp{}: l'efficacité estimée et le revenu de 1963 ont une corrélation de 0,36. C'est la réponse à la première question\nbsp{}: le revenu initial contient une partie de l'efficacité du travail, et cette partie suffit à diviser par quatre la vitesse estimée en coupe. La note précédente avait diagnostiqué le bon mal\nbsp{}; le panel en mesure la dose. #+BEGIN_remarque La section sur [[#efficacite][la mesure de l'efficacité du travail]] a montré que $\log\hat A_i(0)$ est une mauvaise mesure de l'efficacité au sens des comptables de la croissance. Cela n'affaiblit pas l'argument. Ce qui biaise la coupe, c'est la composante *permanente* propre à chaque pays, quel que soit le nom qu'on lui donne, et c'est précisément elle que l'inconnue de pays absorbe. #+END_remarque La troisième question reçoit une réponse plus nuancée, et c'est peut-être le résultat le plus important des deux notes. Le tableau suivant met côte à côte les quatre estimations. #+ATTR_HTML: :class numbers | | Coupe | Panel, intra | Panel, MV | Panel, corrigé | |----------------------------+-------+--------------+-----------+----------------| | $\theta$ | 0,105 | 0,478 | 0,176 | 0,357 | | $\alpha$ | 0,708 | 0,370 | 0,508 | 0,322 | | $1-\alpha$ | 0,292 | 0,630 | 0,492 | 0,678 | | vitesse médiane (% par an) | 0,86 | 3,91 | 1,44 | 2,92 | Corriger l'endogénéité a déplacé *les deux* paramètres, et dans le sens attendu\nbsp{}: la vitesse a plus que triplé, et la part du capital est descendue de 0,71 à 0,32, à hauteur du tiers que mesurent les comptes nationaux. L'endogénéité expliquait donc les *niveaux* aberrants de la coupe. Elle n'explique pas la *tension* entre eux\nbsp{}: $\theta^{\star}$ reste à 0,36 quand le modèle de Solow exige $1-\alpha^{\star} = 0{,}68$, et le test par simulation, qui ne doit rien au $\chi^2$, rejette la contrainte. Le soupçon de la note précédente — que le rejet ne soit qu'un artefact — est levé\nbsp{}: la part du capital qui rend compte des écarts de revenu implique toujours une convergence plus rapide que celle que l'on observe. \\ Ces chiffres éclairent aussi la dispersion des estimations publiées que la note précédente recensait. Les quelque 4\nbsp{}% par an qu'Islam (1995) tire de l'estimateur intra sont ce que nous obtenons avec lui, et la fonction de lien montre qu'ils ne disent presque rien de la vraie vitesse. La valeur corrigée, voisine de 3\nbsp{}%, tombe au milieu de la fourchette de 2 à 4\nbsp{}% que Bond, Hoeffler et Temple (2001) obtiennent par une tout autre voie, les instruments. * Les limites du modèle statistique :PROPERTIES: :CUSTOM_ID: limites :END: Tout ce qui précède repose sur une hypothèse posée au début et jamais vérifiée. Elle porte sur les innovations $\varepsilon_{i,k}$, et sur elles seules\nbsp{}: normales, indépendantes entre elles comme d'un pays à l'autre, et de même variance $\sigma^2$ partout. Tout le reste en découle\nbsp{}: la covariance $\Sigma_i$ des contrastes, donc l'estimateur retenu, et la loi sous l'hypothèse nulle du test de Solow, qui est simulée à partir du modèle. Si cette hypothèse est fausse, ce n'est pas un détail de second ordre. \\ *Comment regarder $\varepsilon$, qui ne s'observe pas\nbsp{}?* En deux étapes. La première donne la perturbation de l'équation\nbsp{}: avec les estimations de la section précédente, $\hat\theta = 0{,}176$ et $\hat\alpha = 0{,}508$, on forme le résidu $r_{i,k}$, qui vaut $g_{i,k}\log A_i(0) + v_{i,k+1}$, et l'on en retire l'inconnue de pays — l'efficacité initiale $\log A_i(0)$, que l'on ne connaît pas — par la régression de $r_{i,k}$ sur $g_{i,k}$\nbsp{}: \begin{equation*} \hat v_{i,k+1} = r_{i,k} - g_{i,k}\widehat{\log A_i(0)} \end{equation*} La seconde étape *blanchit* ce résidu, c'est-à-dire défait la moyenne mobile. La relation $v_k = \varepsilon_k + \psi\varepsilon_{k-1}$ s'inverse de proche en proche, \begin{equation*} \hat\varepsilon_{i,k} = \hat v_{i,k} - \hat\psi\,\hat\varepsilon_{i,k-1} \end{equation*} en partant de $\hat\varepsilon_{i,1} = \hat v_{i,1}$. Nous obtenons ainsi une estimation des innovations elles-mêmes, et ce sont elles que nous soumettrons aux trois hypothèses\nbsp{}: pas de corrélation d'une période à l'autre, variance constante dans le temps, variance constante entre les pays. \\ Une précaution avant de lire les chiffres. Ces innovations estimées ne sont pas les vraies\nbsp{}: l'inconnue de pays y a été remplacée par une estimation tirée des mêmes résidus, ce qui les corrèle mécaniquement, et négativement, même lorsque le modèle est exact\nbsp{}; le blanchiment, qui part d'une condition initiale approchée, ajoute sa propre déformation. Comparer leur autocorrélation à zéro n'aurait donc aucun sens. La bonne référence se simule plutôt qu'elle ne se lit. Nous fabriquons donc 200 jeux de données artificiels — 200 panels de 110 pays et treize dates, produits par le modèle lui-même — et nous leur appliquons *toute* la procédure\nbsp{}: estimation par maximum de vraisemblance, résidus, blanchiment. Les paramètres qui servent à les fabriquer sont ceux que nous tenons pour les plus proches de la vérité, à savoir les valeurs *corrigées* $\theta^{\star} = 0{,}36$ et $\alpha^{\star} = 0{,}32$, avec $\hat\sigma = 0{,}20$, $\hat\psi = 0{,}29$ et les efficacités estimées\nbsp{}: si le modèle est vrai, c'est de là que nos données viennent, et les 200 panels nous disent alors ce que chaque statistique devrait valoir. #+begin_src python :session panel :exports code :results none def blanchi(v, psi): """Innovations retrouvées en inversant la moyenne mobile.""" e = np.empty_like(v) e[0] = v[0] for k in range(1, v.shape[0]): e[k] = v[k] - psi*e[k-1] return e def residus_et_contrastes(Y): """Résidus, innovations blanchies et contrastes, au maximum de vraisemblance de Y.""" th, al, _, ps = maxv(Y, S, C, K, 5, departs=((theta_mv, alpha_mv, omega_mv, psi_mv),)) r, gg = bruts([th, al], Y, S, C, K, 5) A = (gg*r).sum(axis=0)/(gg*gg).sum(axis=0) v, u = r - gg*A, r/gg return v, blanchi(v, ps), u[1:] - u[:-1] def autocorr(x, retard): """Autocorrélation empirique, mise en commun sur les pays.""" return np.corrcoef(x[:-retard].ravel(), x[retard:].ravel())[0, 1] RETARDS = (1, 2, 3, 4) v_obs, e_obs, z_obs = residus_et_contrastes(LY) av_obs = np.array([autocorr(v_obs, l) for l in RETARDS]) # la perturbation v ac_obs = np.array([autocorr(e_obs, l) for l in RETARDS]) # les innovations zc_obs = np.array([autocorr(z_obs, l) for l in (1, 2)]) sp = e_obs.std(axis=0) des = np.random.default_rng(7) AV, AC, ZC, SDP, RAPPORT = [], [], [], [], [] for _ in range(200): Y = simule(*etoile, sigma_mv, psi_mv, lA0_etoile, LY, S, C, K, 5, des) v, e, z = residus_et_contrastes(Y) AV.append([autocorr(v, l) for l in RETARDS]) AC.append([autocorr(e, l) for l in RETARDS]) ZC.append([autocorr(z, l) for l in (1, 2)]) SDP.append(e.std(axis=1)) RAPPORT.append(e.std(axis=0).max()/np.median(e.std(axis=0))) AV, AC, ZC = np.array(AV), np.array(AC), np.array(ZC) SDP, RAPPORT = np.array(SDP), np.array(RAPPORT) print(ac_obs.round(3), AC.mean(axis=0).round(3), np.percentile(AC, [5, 95], axis=0).round(3)) print(av_obs.round(3), AV.mean(axis=0).round(3), np.percentile(AV, [5, 95], axis=0).round(3)) print(zc_obs.round(3), np.percentile(ZC, [5, 95], axis=0).round(3)) print(e_obs.std(axis=1).round(3), np.percentile(SDP, [5, 95]).round(3), sp.min().round(3), codes[sp.argmin()], np.median(sp).round(3), sp.max().round(3), codes[sp.argmax()]) print((sp.max()/np.median(sp)).round(2), np.percentile(RAPPORT, 95).round(2), (e_obs.mean(axis=1).var()/e_obs.var()).round(4)) #+end_src **Les innovations restent corrélées, au troisième retard.** Aux retards 1 à 4, l'autocorrélation des innovations estimées vaut −0,070, −0,066, −0,230 et −0,053, quand le modèle, simulé, en attend −0,098, −0,106, −0,095 et −0,086. Le premier retard est presque rentré dans l'ordre\nbsp{}: −0,070 contre une bande de [−0,115\nbsp{}; −0,083], dont il ne dépasse la borne que de treize millièmes, là où la perturbation $v$ affichait +0,194. C'est le travail de la moyenne mobile, et c'est ce qu'on lui demandait. Le deuxième et le quatrième retards sont dans leur bande. Le troisième, en revanche, vaut −0,230 contre [−0,149\nbsp{}; −0,036]\nbsp{}: il reste là une dépendance que la moyenne mobile ne capte pas, et donc des innovations qui ne sont pas indépendantes. \\ Reste à vérifier l'explication avancée plus haut, celle de l'agrégation temporelle, car rien n'obligeait la corrélation à venir de là plutôt que de l'économie. La prédiction est précise, et donc testable. Échantillonner tous les cinq ans un processus qui évolue continûment n'induit aucune corrélation — les innovations accumulées sur des intervalles disjoints restent indépendantes\nbsp{}; *moyenner* un flux en induit une, positive. Simulons donc le modèle en fréquence annuelle, puis reconstruisons le panel quinquennal des deux façons, par prélèvement ponctuel et par moyenne de période, et comparons. #+begin_src python :session panel :exports code :results none n_an = pop_an[1:]/pop_an[:-1] - 1 c_an = np.clip(n_an + X + delta_an[:-1], 1e-3, None) s_pos = np.clip(s_an[:-1], 1e-3, None) rho_an, a2 = np.exp(-theta_mv*c_an), alpha_mv/(1 - alpha_mv) def annuel(u): """Trajectoire annuelle du même modèle, sans agrégation, pour un choc u donné.""" Y = np.empty((FIN - DEBUT + 1, N)) Y[0] = ly_an[0] for t in range(FIN - DEBUT): Y[t+1] = (rho_an[t]*Y[t] + (1 - rho_an[t])*(a2*(np.log(s_pos[t]) - np.log(c_an[t])) + lA0_mv) + X*(1 + t*(1 - rho_an[t])) + u[t]) return Y def sur_grille(Yan, moyenne, tau=5): """Grille de pas tau : prélèvement ponctuel, ou moyenne autour de la date.""" d, demi = np.arange(0, FIN - DEBUT + 1, tau), tau//2 if not moyenne: return Yan[d, :] return np.array([Yan[max(k-demi, 0):k+demi+1, :].mean(axis=0) for k in d]) for moyenne in (False, True): des, AG = np.random.default_rng(11), [] for _ in range(60): # 0,10 : écart-type annuel qui, agrégé, reproduit le bruit observé bruit = des.normal(0, 0.10, (FIN - DEBUT, N)) e, _, _ = residus_et_contrastes(sur_grille(annuel(bruit), moyenne)) AG.append([autocorr(e, l) for l in RETARDS]) print(moyenne, np.mean(AG, axis=0).round(3)) #+end_src Le verdict est net. Reconstruit par prélèvement ponctuel, le modèle ne produit aucune corrélation positive\nbsp{}; reconstruit par moyennes de période, il en produit une au premier retard et creuse les suivants, c'est-à-dire exactement la signature observée. L'agrégation temporelle n'explique donc pas seulement le signe, elle explique la forme. \\ La prédiction se vérifie une seconde fois, de façon indépendante. Si $\psi$ vient bien de la moyenne de période, il doit s'atténuer lorsque les périodes s'allongent, puisqu'une année moyennée y pèse moins lourd. Refaisons donc l'estimation sur le même échantillon découpé en six périodes de dix ans, et traçons le profil du critère autour de l'optimum pour savoir si le $\psi$ trouvé est identifié ou simplement indéterminé. #+begin_src python :session panel :exports code :results none LY10, S10, C10, K10 = panel(10) # le même échantillon, en six périodes psi10 = maxv(LY10, S10, C10, K10, 10)[3] profil = {psi: minimize(lambda q: logv([q[0], q[1], q[2], np.arctanh(psi/0.95)], LY10, S10, C10, K10, 10), depart(0.2, 0.55)[:3], method='Nelder-Mead').fun for psi in (-0.2, -0.1, 0.02, 0.1, 0.2)} creux = min(profil.values()) print(round(psi10, 3), {k: round(v - creux, 1) for k, v in profil.items()}) #+end_src $\hat\psi$ tombe à 0,02, et le profil montre que ce zéro est bien déterminé\nbsp{}: imposer $\psi = -0{,}1$ coûte 1,6 point de log-vraisemblance, $\psi = +0{,}1$ en coûte 0,8, et $\psi = \pm0{,}2$ plus de quatre. Le paramètre $\psi$ mesure donc une convention de mesure, pas un mécanisme économique. #+begin_src python :session panel :exports none :results none fig, (h1, h2) = plt.subplots(1, 2, figsize=(9, 3.6)) x = np.arange(len(RETARDS)) h1.bar(x - 0.18, ac_obs, width=0.36, color='r', label='observé') h1.bar(x + 0.18, AC.mean(axis=0), width=0.36, color='grey', label='attendu sous le modèle') h1.vlines(x + 0.18, np.percentile(AC, 5, axis=0), np.percentile(AC, 95, axis=0), color='k', linewidth=1) h1.axhline(0, color='k', linewidth=0.8) h1.set_xticks(x, [str(l) for l in RETARDS]) h1.set_xlabel('retard') h1.set_ylabel('autocorrélation des innovations') h1.legend(fontsize=8) debuts = DEBUT + 5*np.arange(K) h2.plot(debuts, e_obs.std(axis=1), 'o-', color='r', label='observé') h2.fill_between(debuts, np.percentile(SDP, 5, axis=0), np.percentile(SDP, 95, axis=0), color='grey', alpha=0.35, label='attendu sous le modèle') h2.set_xlabel('début de période') h2.set_ylabel('écart-type des innovations') h2.legend(fontsize=8) fig.tight_layout() fig.savefig("panel-residus.svg", transparent=True) #+end_src #+CAPTION: *Autocorrélation des innovations estimées et leur écart-type par période, comparés à ce que le modèle estimé produit.* La bande grise couvre 90 % des échantillons simulés. #+LABEL: fig:residus [[file:panel-residus.svg]] \\ *Faut-il alors enrichir le processus d'erreur\nbsp{}?* La question se pose, et elle s'estime\nbsp{}: il suffit d'autoriser des coefficients de moyenne mobile à d'autres retards, ce que la covariance $V$ accepte sans changer une ligne du reste. #+begin_src python :session panel :exports code :results none def mv_retards(retards): """Maximum de vraisemblance avec des coefficients de moyenne mobile aux retards donnés.""" def moins_l(x): psis = np.zeros(max(retards)) psis[np.array(retards) - 1] = 0.95*np.tanh(x[3:]) return critere(np.exp(x[0]), 0.01 + 0.98/(1 + np.exp(-x[1])), np.exp(x[2]), psis, LY, S, C, K, 5) a = (0.508 - 0.01)/0.98 x0 = np.array([np.log(0.176), np.log(a/(1 - a)), np.log(1.272)] + [np.arctanh(0.29/0.95)]*len(retards)) essais = [minimize(moins_l, x0 + d, method='Nelder-Mead', options=dict(maxiter=8000, xatol=1e-7, fatol=1e-9)) for d in (0.0, 0.1)] x = min(essais, key=lambda e: e.fun).x return (-min(essais, key=lambda e: e.fun).fun, np.exp(x[0]), 0.01 + 0.98/(1 + np.exp(-x[1])), (0.95*np.tanh(x[3:])).round(3)) for retards in ([1], [1, 2, 3], [1, 3]): l, th, al, ps = mv_retards(retards) print(retards, round(l, 2), round(th, 4), round(al, 4), ps) #+end_src #+ATTR_HTML: :class numbers | erreurs | log-vraisemblance | $\hat\theta$ | $\hat\alpha$ | coefficients | |-+-+-+-+-| | MA(1) | 2042,27 | 0,176 | 0,508 | $\psi_1 = 0{,}294$ | | MA(3) | 2050,19 | 0,179 | 0,516 | 0,299 ; 0,004 ; −0,123 | | MA aux retards 1 et 3 | 2050,18 | 0,179 | 0,516 | 0,297 ; −0,125 | Trois enseignements. D'abord, le troisième retard est bien réel\nbsp{}: le gagner vaut 15,8 points de rapport de vraisemblance pour un seul paramètre, ce qu'aucun hasard ne produit. Ensuite, le deuxième retard ne sert à rien\nbsp{}: le MA(3) complet et sa version contrainte à $\psi_2 = 0$ atteignent la même vraisemblance, à un centième près, la seconde avec un paramètre de moins. La spécification que les données préfèrent est donc $v_{i,k} = \varepsilon_{i,k} + \psi_1\varepsilon_{i,k-1} + \psi_3\varepsilon_{i,k-3}$, avec $\hat\psi_1 = 0{,}297$ et $\hat\psi_3 = -0{,}125$. Enfin, et c'est ce qui importe pour cette note, le paramètre d'intérêt ne bouge pas\nbsp{}: 0,176 avec le MA(1), 0,179 avec le retard supplémentaire, soit le douzième de l'écart-type de $\hat\theta$. Nous gardons donc le MA(1) comme référence — l'agrégation des flux le justifie, il est plus parcimonieux, et les conclusions de la note n'en dépendent pas — mais il faut savoir qu'il est rejeté par les données. \\ *D'où peut venir une dépendance à trois périodes\nbsp{}?* Nous n'en savons rien. Trois périodes de cinq ans font quinze ans, et nous ne voyons aucun mécanisme, économique ou de mesure, qui produirait une corrélation à cet horizon-là et à lui seul. Nous laissons donc le fait tel quel\nbsp{}: c'est une réserve sur le modèle statistique, pas un résultat de la note — et, comme le tableau ci-dessus le montre, la vitesse de convergence n'en dépend pas. **La variance n'est pas constante dans le temps.** L'écart-type des innovations estimées passe de 0,143 sur la période ouverte en 1973 à 0,316 sur celle ouverte en 2013, quand la bande simulée va de 0,178 à 0,223, et les moyennes de période expliquent 8,0\nbsp{}% de leur variance. Le modèle ne prévoit rien de tel, puisque sa tendance $\gamma_{i,k}$ est commune mais déterministe. C'est un argument sérieux en faveur des indicatrices temporelles que la section [[#panel][De la coupe au panel]] écartait un peu vite\nbsp{}: elles coûtent $K-1$ paramètres, mais elles absorbent ces chocs communs, et la comparaison d'ajustement de la section sur [[#intra][l'estimateur intra]] montre qu'elles les paient largement, puisqu'elles doublent la part expliquée de la croissance. **Elle ne l'est pas davantage entre les pays.** L'écart-type des innovations va de 5,0\nbsp{}% en Australie à 95,9\nbsp{}% au Venezuela, pour une médiane de 13,5\nbsp{}%. Le rapport du maximum à la médiane vaut 7,1, quand le modèle ne le fait dépasser 1,8 que dans 5\nbsp{}% des échantillons. Un seul $\sigma$ pour 110 pays est une fiction commode. \\ Le graphique laisse deviner une régularité, et elle se vérifie\nbsp{}: les économies riches sont les plus calmes. La corrélation entre $\log\hat\sigma_i$ et le revenu par tête moyen du pays, lui aussi en logarithme, vaut −0,27, soit un rapport de Student de −2,9, et elle résiste au retrait des cinq pays les plus turbulents. La pente vaut −0,12\nbsp{}: doubler le revenu par tête réduit l'écart-type des innovations de 8\nbsp{}%. Ce n'est toutefois pas un dégradé régulier. Les médianes de $\hat\sigma_i$ par quartile de revenu valent 14,7\nbsp{}%, 15,0\nbsp{}%, 15,2\nbsp{}% et 8,6\nbsp{}%\nbsp{}: les trois premiers quarts sont indiscernables et seul le plus riche se détache, si bien que le niveau de développement ne rend compte que de 7\nbsp{}% de la dispersion. Nous n'avons pas davantage d'explication à proposer ici qu'au paragraphe précédent — diversification des productions, stabilisation macroéconomique, qualité des comptes nationaux sont toutes plausibles, et aucune n'est testée dans cette note. #+begin_src python :session panel :exports code :results none niveau = LY.mean(axis=0) # revenu par tête moyen du pays, en logarithme pente_revenu, ordonnee = np.polyfit(niveau, np.log(sp), 1) quartile = np.digitize(niveau, np.quantile(niveau, [0.25, 0.5, 0.75])) residu = np.log(sp) - (ordonnee + pente_revenu*niveau) print(round(np.corrcoef(niveau, np.log(sp))[0, 1], 3), round(pente_revenu, 3), round(1 - residu.var()/np.log(sp).var(), 3), [round(100*float(np.median(sp[quartile == j])), 1) for j in range(4)]) #+end_src #+begin_src python :session panel :exports none :results none fig, ax = plt.subplots(figsize=(6, 4.5)) ax.plot(np.exp(niveau), 100*sp, '.', color='grey') for pays, decal in zip(('VEN', 'FRA', 'AUS', 'USA'), (6, -12, -12, 6)): i = list(codes).index(pays) ax.plot(np.exp(niveau[i]), 100*sp[i], 'o', color='r') ax.annotate(pays, (np.exp(niveau[i]), 100*sp[i]), xytext=(6, decal), textcoords='offset points', fontsize=8, color='r') bornes = np.array([niveau.min(), niveau.max()]) ax.plot(np.exp(bornes), 100*np.exp(ordonnee + pente_revenu*bornes), color='b', linewidth=1, label='ajustement, pente −0,12') ax.axhline(100*sigma_mv, color='k', linestyle='--', linewidth=1, label=r"écart-type commun supposé, $\sigma$") ax.set_xscale('log') ax.set_yscale('log') ax.set_xlabel('revenu par tête moyen du pays (dollars, échelle logarithmique)') ax.set_ylabel('écart-type des innovations (%, échelle logarithmique)') ax.legend(fontsize=8) fig.tight_layout() fig.savefig("panel-pays.svg", transparent=True) #+end_src #+CAPTION: *Écart-type des innovations estimées contre le revenu par tête moyen du pays.* Les deux échelles sont logarithmiques : un facteur de plus de vingt sépare l'Australie du Venezuela, alors que le modèle suppose un écart-type unique. #+LABEL: fig:pays [[file:panel-pays.svg]] \\ Reste une objection, qu'il faut examiner parce qu'elle semble tout absoudre. La vraisemblance ne porte pas sur les innovations, mais sur les contrastes \begin{equation*} z_{i,k} = u_{i,k} - u_{i,k-1} = \frac{v_{i,k+1}}{g_{i,k}} - \frac{v_{i,k}}{g_{i,k-1}} \end{equation*} dont elle suppose la covariance $\sigma^2\Sigma_i = \sigma^2 T_i V T_i'$. Contrôlons donc cet objet-là plutôt que les innovations\nbsp{}: sa corrélation vaut −0,49 au premier retard et −0,01 au second, et les deux tombent dans les bandes simulées, [−0,50\nbsp{}; −0,32] et [−0,19\nbsp{}; 0,00]. Ce que l'estimateur regarde est donc compatible avec les données. \\ Il serait tentant d'y voir un démenti de tout ce qui précède. Ce serait mal lire, pour deux raisons. La première est que ce contrôle a peu de puissance\nbsp{}: différencier impose à soi seul une corrélation voisine de −0,5 au premier retard, et les bandes sont larges. La seconde est plus fondamentale. Tout, dans ces contrastes, vient des innovations\nbsp{}: \begin{equation*} z_i = T_i r_i,\qquad r_i = g_i\log A_i(0) + M\varepsilon_i \end{equation*} où $M$ est la matrice de la moyenne mobile, celle qui vérifie $V = MM'$. Comme $T_i$ élimine $\log A_i(0)$ de toutes les composantes sauf la première, les contrastes sont, à celle-ci près, une fonction linéaire des seules $\varepsilon_{i,k}$. Les hypothèses portent donc en dernier ressort sur $\varepsilon$, et sur lui seul\nbsp{}: si $\varepsilon$ les viole, $z$ les viole aussi. Un contrôle sur $z$ qui ne décèle rien n'innocente pas $\varepsilon$\nbsp{}; il dit seulement que la transformation — différencier, diviser par $g_{i,k}$, mettre en commun des pays d'écarts-types très inégaux — a dilué la trace au point de la rendre indétectable. \\ La bonne question est donc ailleurs\nbsp{}: que devient l'estimateur lorsque la loi supposée pour $\varepsilon$ est fausse, alors que l'équation qui en donne la moyenne, elle, reste correcte\nbsp{}? La vraisemblance cesse d'être la vraie vraisemblance et n'est plus qu'un critère parmi d'autres — on parle de *pseudo*-vraisemblance. Elle garde un maximum, mais la courbe de lien qui sert à corriger le biais et la variance qui sert à tester sont l'une et l'autre calculées sous une loi qui n'est pas celle des données. #+BEGIN_remarque Que faut-il en conclure pour les résultats de la note\nbsp{}? La correction par inférence indirecte garde son sens\nbsp{}: elle ne suppose rien d'autre que la capacité du modèle à reproduire une statistique auxiliaire, et elle serait seulement moins efficace si la loi des erreurs est mal spécifiée. Le test de Solow est plus exposé, puisque sa loi sous l'hypothèse nulle est simulée avec ces mêmes erreurs mal spécifiées. Son verdict reste vraisemblablement bon — l'écart observé, −0,32, vaut plus de quatre fois la borne de la bande simulée, ±0,07, et il faudrait une erreur de spécification énorme pour combler un tel fossé — mais « vraisemblablement » est ici le mot juste, et une version robuste du test resterait à écrire. Les trois remèdes naturels sont un $\sigma_i$ par pays, des effets temporels communs, et une variance robuste à la mauvaise spécification, dite en sandwich, pour l'inférence\nbsp{}; aucun n'est gratuit, et tous mériteraient leur propre note. #+END_remarque * Des vitesses entièrement libres :PROPERTIES: :CUSTOM_ID: libres :END: Avant de conclure, un dernier point de comparaison, qui situe la restriction structurelle par rapport à l'autre extrême. Tout ce qui précède impose au pays $i$ la vitesse $\theta\,c_{i,k}$\nbsp{}; on peut au contraire n'imposer rien du tout et estimer une vitesse par pays, sur les données annuelles. #+begin_src python :session panel :exports code :results none def vitesses_libres(tendance): """Une régression par pays sur données annuelles, et la vitesse implicite.""" rho = np.empty(N) for i in range(N): cols = [np.ones(FIN - DEBUT), ly_an[:-1, i]] if tendance: cols.append(np.arange(FIN - DEBUT)) rho[i] = np.linalg.lstsq(np.column_stack(cols), ly_an[1:, i], rcond=None)[0][1] return rho libres = {} for tendance in (False, True): rho = vitesses_libres(tendance) fini = rho < 1 beta = -np.log(np.clip(rho, 1e-6, None)) libres[tendance] = (beta, fini) A = np.column_stack([np.ones(fini.sum()), cbar[fini]]) co = np.linalg.lstsq(A, beta[fini], rcond=None)[0] res = beta[fini] - A @ co V = np.linalg.inv(A.T @ A)*(res @ res)/(fini.sum() - 2) print(tendance, (~fini).sum(), np.median(beta[fini]).round(4), beta[fini].mean().round(4), beta[fini].std().round(4), co[1].round(3), np.sqrt(V[1, 1]).round(3)) #+end_src Le résultat est un avertissement. Sans tendance, 26 pays sur 110 ont un $\hat\rho_i$ supérieur ou égal à un, c'est-à-dire ne convergent pas du tout\nbsp{}; pour les autres, la vitesse médiane est de 2,7\nbsp{}% par an, la moyenne de 5,0\nbsp{}% et l'écart-type de 5,8 points. La régression de ces vitesses sur $\bar c_i$, que le modèle voudrait de pente $\theta$, donne 0,45 avec un écart-type de 0,41\nbsp{}: la prédiction structurelle n'est ni confirmée ni infirmée, l'estimation libre est trop imprécise pour trancher. Et tout bascule si l'on ajoute une tendance à chaque régression\nbsp{}: il ne reste que 3 pays non convergents, la vitesse médiane monte à 8,4\nbsp{}% et la pente tombe à 0,14 avec un écart-type de 0,55. Lee, Pesaran et Smith (1997) ont fait la même expérience en allant plus loin\nbsp{}: en laissant libres, pays par pays, le niveau et la tendance, ils obtiennent une vitesse de l'ordre de 30\nbsp{}% par an, qu'ils jugent eux-mêmes très imprécise et difficile à interpréter. Le sens est le même que le nôtre\nbsp{}: plus on libère de paramètres par pays, plus la vitesse estimée monte, et moins elle est précise. #+begin_src python :session panel :exports none :results none beta_l, fini_l = libres[False] fig, ax = plt.subplots(figsize=(6, 4.5)) ax.scatter(100*cbar[fini_l], 100*beta_l[fini_l], s=10, color='grey', label='vitesses estimées pays par pays') grille_c = np.linspace(cbar.min(), cbar.max(), 20) ax.plot(100*grille_c, 100*etoile[0]*grille_c, 'r-', linewidth=1, label=r'prédiction du modèle, $\theta^{\star}\bar c_i$') ax.set_xlabel(r'$\bar c_i$ (% par an)') ax.set_ylabel(r'vitesse de convergence estimée (% par an)') ax.legend(fontsize=8) fig.tight_layout() fig.savefig("panel-libres.svg", transparent=True) #+end_src #+CAPTION: *Vitesses estimées pays par pays, sans contrainte, et prédiction du modèle.* Les 26 pays dont l'estimation n'implique aucune convergence ne figurent pas sur le graphique. #+LABEL: fig:libres [[file:panel-libres.svg]] \\ On mesure ici ce que la restriction structurelle achète. Arnold, Bassanini et Scarpetta (2011) partent de la même observation que cette note — le modèle implique des vitesses différentes d'un pays à l'autre — mais en tirent la conclusion opposée\nbsp{}: ils estiment une vitesse libre par pays, et se servent de $n_i+x+\delta_i$ comme d'une simple variable de contrôle. La relation structurelle est alors une motivation, jamais une contrainte. En l'imposant, on paie une hypothèse et on gagne deux choses\nbsp{}: des vitesses estimables, et un test. * Conclusion Le panel tient l'essentiel de ce qu'il promettait. La note en coupe s'achevait sur un doute qu'elle ne pouvait lever\nbsp{}: son rejet du modèle de Solow venait-il du modèle, ou de l'endogénéité du revenu initial\nbsp{}? La réponse est maintenant double. L'endogénéité est réelle et mesurable\nbsp{}: l'efficacité du travail récupérée des inconnues de pays croît avec le taux d'épargne, ce qui réfute l'hypothèse de Mankiw, Romer et Weil, et elle est corrélée au revenu initial, ce qui suffit à diviser par quatre la vitesse estimée en coupe — des panels simulés avec la vitesse corrigée redonnent, vus en coupe, l'estimation de la note précédente. Mais le rejet du modèle de Solow n'en était pas un artefact\nbsp{}: une fois l'endogénéité éliminée, la part du capital tombe à 0,32 et la vitesse monte à 2,9\nbsp{}% par an, sans que l'une rejoigne l'autre. \\ Le chemin a demandé plus qu'un changement de données. Le panel introduit un problème que la coupe n'avait pas, et d'une sévérité qu'on sous-estime en récitant l'ordre de grandeur de Nickell\nbsp{}: au niveau de bruit de ces données, l'estimateur intra est presque muet sur la vitesse de convergence. L'inférence indirecte le contourne, à condition d'inverser un estimateur qui réponde au paramètre, même biaisé. Le verdict qui en sort ne repose plus sur l'hypothèse d'exogénéité ni sur une approximation asymptotique. Il repose en revanche sur la loi que l'on prête aux perturbations. Ces perturbations sont corrélées d'une période à l'autre, la [[#limites][section sur les limites]] montre que la cause en est la façon dont les données sont datées plutôt qu'un mécanisme économique, et le modèle en tient compte. Ce qu'il ne corrige pas, c'est l'hétéroscédasticité\nbsp{}: l'écart-type des perturbations varie dans un rapport de plus de vingt d'un pays à l'autre, et le rejet de Solow, quoique large, n'a pas été éprouvé contre cette mauvaise spécification-là. Il reste enfin à soumettre le modèle aux instruments, ce qui fera l'objet d'une prochaine note. * Exercices :PROPERTIES: :CUSTOM_ID: exercices :END: #+NAME: ex-tendance #+BEGIN_exercice Établissez que $\gamma_{i,k} = x\bigl[\tau + t_k\,g_{i,k}\bigr]$. Que suppose-t-on en la remplaçant par une indicatrice par période, comme il est d'usage\nbsp{}? #+END_exercice #+ATTR_HTML: :class corrige #+BEGIN_details #+BEGIN_summary Corrigé #+END_summary En écrivant $\log y_{i,k} = \log A_i(0) + xt_k + \log\hat y_{i,k}$ et en reportant dans la dynamique de $\log\hat y$, le terme en $x$ qui subsiste au membre de droite est $xt_{k+1} - \rho_{i,k} xt_k$. Comme $t_{k+1} = t_k+\tau$, il vaut $x\tau + xt_k(1-\rho_{i,k})$. Il dépend donc du pays et de la période, par $\rho_{i,k}$, et croît avec la date. Une indicatrice par période impose au contraire le même terme à tous les pays\nbsp{}: elle absorbe les chocs communs, ce qui est souvent souhaitable, mais aussi une partie de la dynamique structurelle, et elle consomme $K-1$ paramètres que le modèle fournissait gratuitement. #+END_details #+NAME: ex-ponderation #+BEGIN_exercice Montrez que la régression de $r_{i,k}$ sur $g_{i,k}$ se réduit à la moyenne temporelle des résidus lorsque les taux sont constants. Vérifiez numériquement que les deux écritures donnent alors la même estimation, puis mesurez l'écart lorsque les taux varient. #+END_exercice #+ATTR_HTML: :class corrige #+BEGIN_details #+BEGIN_summary Corrigé #+END_summary Si $g_{i,k} = g_i$ ne dépend pas de $k$, alors $\sum_k g_i r_{i,k}/\sum_k g_i^2 = \bar r_i/g_i$, et le résidu concentré $r_{i,k} - g_i\bar r_i/g_i = r_{i,k} - \bar r_i$ est exactement celui de la transformation intra. #+begin_src python :session panel :exports code :results none def intra_naif(LY, S, C, K, tau): """Variante fautive : moyenne simple des résidus, sans le coefficient g.""" def critere(p): r, _ = bruts(p, LY, S, C, K, tau) return (r - r.mean(axis=0)).ravel() return min([least_squares(critere, p0, bounds=([1e-4, 0.01], [5, 0.99])) for p0 in ([0.4, 0.4], [0.1, 0.7], [0.03, 0.5])], key=lambda e: e.cost) Cfixe = np.tile(C.mean(axis=0), (K, 1)) print(intra(LY, S, Cfixe, K, 5).x.round(4), intra_naif(LY, S, Cfixe, K, 5).x.round(4)) print(intra(LY, S, C, K, 5).x.round(4), intra_naif(LY, S, C, K, 5).x.round(4)) #+end_src Avec des taux figés à leur moyenne, les deux estimateurs coïncident aux quatre décimales. Avec des taux par période, la version qui ignore $g_{i,k}$ s'effondre à 0,07 contre 0,48\nbsp{}: elle n'élimine pas l'inconnue de pays, et l'écart n'est pas un détail numérique mais une erreur de spécification. #+END_details #+NAME: ex-jacobien #+BEGIN_exercice La vraisemblance de la section sur [[#condition][le maximum de vraisemblance]] porte le terme $\sum_{k\geq2}\log g_{i,k}$. D'où vient-il, et que se passe-t-il si on l'oublie\nbsp{}? Vérifiez qu'il se réduit à une constante lorsque les taux sont figés. #+END_exercice #+ATTR_HTML: :class corrige #+BEGIN_details #+BEGIN_summary Corrigé #+END_summary C'est le jacobien du changement de variable. La vraisemblance est écrite pour le vecteur transformé $z$, alors que les données sont les résidus $r$\nbsp{}; la densité du premier est celle du second divisée par la valeur absolue du déterminant de la transformation, ici $\prod_{k\geq2} 1/g_{i,k}$. Comme les $g$ dépendent de $\theta$, ce terme ne se réduit pas à une constante et doit être porté. L'omettre offre à l'optimiseur une façon gratuite d'améliorer le critère\nbsp{}: gonfler les $g$ augmente le déterminant sans rien expliquer, et $\theta$ part à sa borne supérieure. Lorsque les taux sont constants, $g_{i,k} = g_i$ et le terme devient $(K-1)\sum_i\log g_i$, exactement compensé par la variation du déterminant de $\Sigma_i$\nbsp{}; les deux écritures donnent alors la même estimation. #+END_details #+NAME: ex-nickell #+BEGIN_exercice Calculez l'ordre de grandeur du biais de Nickell, $-(1+\rho)/(K-1)$, pour $\rho = 0{,}84$ et $K = 12$, puis traduisez-le en vitesse de convergence. Comparez au Monte Carlo de la section sur [[#biais][le biais de l'estimateur intra]]. #+END_exercice #+ATTR_HTML: :class corrige #+BEGIN_details #+BEGIN_summary Corrigé #+END_summary Le biais vaut $-1{,}84/11 \simeq -0{,}17$, soit $\hat\rho \simeq 0{,}67$ au lieu de 0,84, et une vitesse de $-\log(0{,}67)/5 = 8{,}0\nbsp{}\%$ par an au lieu de 3,4\nbsp{}%. Le sens est le bon, mais le Monte Carlo montre que l'analogie est trompeuse\nbsp{}: dans le modèle linéaire, le biais est une translation que l'on peut retrancher, alors qu'ici la fonction qui relie le vrai paramètre à l'espérance de son estimateur est plate sur tout un intervalle. Ce n'est pas la même maladie, et elle ne se soigne pas de la même façon. #+END_details #+NAME: ex-convergence #+BEGIN_exercice Un panel peut grandir dans deux directions\nbsp{}: en ajoutant des pays ou en ajoutant des périodes. Simulez des panels de taille croissante dans chacune d'elles, en répliquant les taux et les efficacités observés, et appliquez-leur les deux estimateurs. Dans quelle direction le biais s'efface-t-il, et pourquoi\nbsp{}? Qu'en conclure pour nos données\nbsp{}? #+END_exercice #+ATTR_HTML: :class corrige #+BEGIN_details #+BEGIN_summary Corrigé #+END_summary #+begin_src python :session panel :exports code :results none grandir, VRAI = np.random.default_rng(19), 0.176 # mêmes paramètres que la fonction de lien def agrandi(m, q): """Panel simulé avec m copies des pays et q copies des périodes.""" Sg, Cg, lA0g = np.tile(S, (q, m)), np.tile(C, (q, m)), np.tile(lA0_mv, m) Kg, gabarit = q*K, np.empty((q*K + 1, m*N)) gabarit[0] = np.tile(LY[0], m) Y = simule(VRAI, 0.508, sigma_mv, psi_mv, lA0g, gabarit, Sg, Cg, Kg, 5, grandir) return Y, Sg, Cg, Kg for m, q in ((1, 1), (2, 1), (8, 1), (1, 2), (1, 4), (1, 8)): W, M = [], [] for _ in range(15): Y, Sg, Cg, Kg = agrandi(m, q) W.append(intra(Y, Sg, Cg, Kg, 5).x[0]) M.append(maxv(Y, Sg, Cg, Kg, 5, departs=((VRAI, 0.508, omega_mv, psi_mv),))[0]) print(110*m, Kg, round(np.mean(W), 3), round(np.std(W), 3), round(np.mean(M), 3), round(np.std(M), 3)) #+end_src Le tableau donne la moyenne des estimations et, entre parenthèses, leur écart-type, pour un vrai $\theta$ de 0,176 et les paramètres de nuisance de la correction, efficacités du maximum de vraisemblance comprises. La colonne MV est donc comparable à la fonction de lien du corps de la note\nbsp{}; celle de l'estimateur intra ne l'est pas, le tableau du biais ayant été simulé, faute de mieux, avec les efficacités de l'intra. #+ATTR_HTML: :class numbers | pays | périodes | intra | MV | |-+-+-+-| | 110 | 12 | 0,365 (0,038) | 0,112 (0,028) | | 220 | 12 | 0,388 (0,040) | 0,107 (0,024) | | 880 | 12 | 0,384 (0,012) | 0,106 (0,010) | | 110 | 24 | 0,258 (0,018) | 0,160 (0,033) | | 110 | 48 | 0,190 (0,013) | 0,175 (0,019) | | 110 | 96 | 0,155 (0,005) | 0,173 (0,007) | *Multiplier les pays ne corrige rien.* De 110 à 880 pays, les estimations se resserrent — l'écart-type est divisé par trois — mais autour d'une valeur fausse\nbsp{}: 0,384 pour l'intra et 0,106 pour le maximum de vraisemblance, quand la vérité est 0,176. C'est le propre des paramètres incidents\nbsp{}: chaque pays supplémentaire apporte douze observations et une inconnue de plus, si bien que $\log A_i(0)$ n'est jamais mieux estimé. Or c'est la direction dans laquelle nos données ressemblent le plus à un grand échantillon. *Allonger la durée fonctionne, pour le maximum de vraisemblance.* À nombre de pays donné, son estimation passe de 0,112 à 0,160, puis 0,175 et 0,173, et son écart-type tombe de 0,028 à 0,007. L'estimateur intra, lui, descend de 0,365 à 0,155 en traversant la vraie valeur au passage\nbsp{}: il ignore la corrélation des erreurs, qui est pour lui une mauvaise spécification, et une mauvaise spécification ne s'efface pas quand l'échantillon grandit. Deux précautions de lecture. Voir une distribution se resserrer sur une grille de tailles n'est pas une démonstration de convergence en probabilité, et une moyenne qui se rapproche n'y suffirait pas — encore faut-il que la dispersion s'annule, ce que l'écart-type permet de surveiller. Et les deux expériences sont artificielles, puisqu'elles fabriquent des données que nous n'avons pas. La conclusion pratique n'en est pas moins claire\nbsp{}: avec douze périodes, il faut corriger le biais plutôt que l'espérer négligeable. #+END_details #+NAME: ex-marche #+BEGIN_exercice Supposons que chaque pays suive une marche aléatoire avec dérive, $y_{i,t} = y_{i,t-1} + \mu + \varepsilon_{i,t}$, les $\varepsilon$ étant indépendants, d'espérance nulle et de variance $\sigma^2$. Il n'y a donc aucune convergence. Montrez que la pente intra tend, lorsque le nombre de pays grandit, vers \begin{equation*} 1 - \frac{3\sigma^2}{(T+1)\left(\sigma^2 + \mu^2 T/2\right)} \end{equation*} où $T$ est le nombre de transitions. Que vaut cette limite sans dérive\nbsp{}? Que devient-elle quand $T$ grandit\nbsp{}? Vérifiez-la par simulation. #+END_exercice #+ATTR_HTML: :class corrige #+BEGIN_details #+BEGIN_summary Corrigé #+END_summary Écrivons $x_{i,t} = y_{i,t-1}$ le régresseur et $\tilde x$ son écart à la moyenne du pays. Comme $y_t = x_t + \mu + \varepsilon_t$ et que la constante disparaît au centrage, \begin{equation*} \hat\rho - 1 = \frac{\sum_i\sum_t \tilde x_{i,t}\,\tilde\varepsilon_{i,t}}{\sum_i\sum_t \tilde x_{i,t}^2} \end{equation*} Divisons haut et bas par le nombre de pays. Chaque somme devient une moyenne de termes indépendants et de même loi, un par pays. Ces termes sont des formes quadratiques en des chocs de variance finie, donc d'espérance finie\nbsp{}: la loi des grands nombres s'applique séparément au numérateur et au dénominateur, et chacun converge en probabilité vers l'espérance d'un terme. Celle du dénominateur étant strictement positive, le lemme de Slutsky autorise le passage au rapport. Deux précautions sur cette étape. D'abord, pour une suite indépendante et de même loi, l'existence de l'espérance suffit à la loi des grands nombres. La démonstration élémentaire, par l'inégalité de Bienaymé-Tchebychev — la variance d'une moyenne de $N$ termes indépendants vaut celle d'un terme divisée par $N$, et tend donc vers zéro — est parfaitement valable, mais elle exige en plus des chocs de moment d'ordre quatre fini, ce que l'énoncé ne suppose pas. Ensuite, c'est bien une limite en probabilité quand le nombre de pays grandit, à $T$ fixé, et non l'espérance de $\hat\rho$, qui n'est pas le rapport des espérances. *Le numérateur.* On a $\sum_t \tilde x_t\tilde\varepsilon_t = \sum_t x_t\varepsilon_t - T\bar x\bar\varepsilon$. Le premier terme est d'espérance nulle, puisque $x_t$ ne dépend que des chocs antérieurs à $t$. Le second ne l'est pas, et c'est tout le problème\nbsp{}: la moyenne $\bar x$ contient les chocs, avec $\mathbb{E}[\bar x\bar\varepsilon] = \sigma^2(T-1)/(2T)$, de sorte que \begin{equation*} \mathbb{E}\Bigl[\sum_t \tilde x_t\tilde\varepsilon_t\Bigr] = -\frac{\sigma^2(T-1)}{2} \end{equation*} Retrancher la moyenne du pays, c'est retrancher une quantité qui contient les chocs futurs\nbsp{}: le régresseur transformé et l'erreur transformée deviennent corrélés, négativement. *Le dénominateur.* La dérive et les chocs y contribuent séparément. Pour la partie aléatoire, $\mathbb{E}[x_tx_s] = \sigma^2\min(t-1,s-1)$ et $\sum_{a,b=0}^{T-1}\min(a,b) = (T-1)T(2T-1)/6$, d'où une contribution de $\sigma^2(T^2-1)/6$. Pour la dérive, $x_t$ contient $\mu(t-1)$, dont l'écart à la moyenne a pour somme des carrés $\mu^2T(T^2-1)/12$. Le rapport des deux limites donne la formule annoncée. *Les deux cas.* Sans dérive, elle se réduit à $1 - 3/(T+1)$, soit 0,769 pour nos douze transitions\nbsp{}: un biais considérable. Avec la dérive, le dénominateur grossit et le biais fond. Et quand $T$ grandit, la limite tend vers un dans les deux cas, en $1/T$\nbsp{}: c'est encore la leçon de l'[[ex-convergence][exercice précédent]], allonger la durée aide, multiplier les pays non. #+begin_src python :session panel :exports code :results none mu_rw, sig_rw = dy.mean(), (dy - dy.mean(axis=0)).std() limite = lambda mu, sig: 1 - 3*sig**2/((K + 1)*(sig**2 + mu**2*K/2)) print(round(mu_rw, 4), round(sig_rw, 4), round(limite(0, sig_rw), 4), round(limite(mu_rw, sig_rw), 4)) essai = np.random.default_rng(5) for derive in (0.0, mu_rw): P = [] for _ in range(400): Y = np.vstack([np.zeros(N), np.cumsum(derive + essai.normal(0, sig_rw, (K, N)), axis=0)]) av, ap = Y[:-1, :], Y[1:, :] P.append(pente(av - av.mean(axis=0), ap - ap.mean(axis=0))) print(round(derive, 4), round(float(np.mean(P)), 4)) #+end_src Sur nos données, la dérive vaut 0,1025 par période de cinq ans et l'écart-type des chocs 0,2074. La formule donne alors 0,9064, et quatre cents panels simulés en donnent 0,9063 — la pente observée sur les vraies données valant 0,906. Sans dérive, la formule donne 0,769 et la simulation 0,7691. Le biais de l'estimateur intra sur un monde sans convergence n'est donc pas un mystère numérique\nbsp{}: il s'écrit. #+END_details * Références L'hypothèse d'orthogonalité discutée ici est celle de Mankiw, N. G., Romer, D. et Weil, D. (1992), « A Contribution to the Empirics of Economic Growth », /Quarterly Journal of Economics/, 107(2), 407-437. Le cadre de cette note est celui de Islam, N. (1995), « Growth Empirics: A Panel Data Approach », /Quarterly Journal of Economics/, 110(4), 1127-1170, qui le premier a traité la régression de convergence comme un panel dynamique. Sur le biais de l'estimateur intra, voir Nickell, S. (1981), « Biases in Dynamic Models with Fixed Effects », /Econometrica/, 49(6), 1417-1426, et, pour sa correction analytique, Kiviet, J. (1995), « On Bias, Inconsistency and Efficiency of Various Estimators in Dynamic Panel Data Models », /Journal of Econometrics/, 68(1), 53-78. La vraisemblance en différences est celle de Hsiao, C., Pesaran, M. H. et Tahmiscioglu, A. K. (2002), « Maximum Likelihood Estimation of Fixed Effects Dynamic Panel Data Models Covering Short Time Periods », /Journal of Econometrics/, 109(1), 107-150. L'inférence indirecte est due à Gouriéroux, C., Monfort, A. et Renault, E. (1993), « Indirect Inference », /Journal of Applied Econometrics/, 8(S1), S85-S118, et son application au panel dynamique à Gouriéroux, C., Phillips, P. C. B. et Yu, J. (2010), « Indirect Inference for Dynamic Panel Models », /Journal of Econometrics/, 157(1), 68-77. \\ Sur l'hétérogénéité des vitesses, voir Lee, K., Pesaran, M. H. et Smith, R. (1997), « Growth and Convergence in a Multi-Country Empirical Stochastic Solow Model », /Journal of Applied Econometrics/, 12(4), 357-392, et Pesaran, M. H. et Smith, R. (1995), « Estimating Long-Run Relationships from Dynamic Heterogeneous Panels », /Journal of Econometrics/, 68(1), 79-113. La comparaison de la fin de la note porte sur Arnold, J., Bassanini, A. et Scarpetta, S. (2011), « Solow or Lucas? Testing Speed of Convergence on a Panel of OECD Countries », /Research in Economics/, 65(2), 110-123. Une voie entièrement différente, qui contourne les effets fixes en passant par la dynamique du ratio capital/production, est celle de McQuinn, K. et Whelan, K. (2007), « Conditional Convergence and the Dynamics of the Capital-Output Ratio », /Journal of Economic Growth/, 12(2), 159-184. La quasi-différence des panneaux à coefficients variables est celle de Holtz-Eakin, D., Newey, W. et Rosen, H. (1988), « Estimating Vector Autoregressions with Panel Data », /Econometrica/, 56(6), 1371-1395. Les estimations par instruments citées dans le retour à la coupe sont celles de Bond, S., Hoeffler, A. et Temple, J. (2001), « GMM Estimation of Empirical Growth Models », /CEPR Discussion Paper/ 3048. La corrélation induite par l'agrégation temporelle est établie par Working, H. (1960), « Note on the Correlation of First Differences of Averages in a Random Chain », /Econometrica/, 28(4), 916-918. \\ Les données proviennent de la version 11.0 des Penn World Tables (Groningen Growth and Development Centre, 2025, [[https://doi.org/10.34894/FABVLR][doi:10.34894/FABVLR]]), décrites dans Feenstra, R., Inklaar, R. et Timmer, M. (2015), « The Next Generation of the Penn World Table », /American Economic Review/, 105(10), 3150-3182. \\ \\ \\