Les vitesses de convergence en panel

On reprend les notations de la note précédente : \(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 : 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}\) ; 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 : \(\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

où \(g_{i,k} = 1-\rho_{i,k}\) et où

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 : 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 : 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 : en retranchant \(\log y_{i,k}\) des deux côtés,

on retrouve exactement la forme de la note en coupe : 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 : 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 : 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 : \(\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 : 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.

Nous reprenons les Penn World Tables 11.0 et l'échantillon de la note précédente : 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.

La période retenue commence en 1963 et non en 1960 : 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.

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

L'échantillon compte 110 pays et douze périodes. Le \(c_{i,k}\) médian vaut 8,1 % 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 : la section sur l'estimateur intra traite le cas général, celui de taux qui varient d'une période à l'autre.

Les notes précédentes travaillaient sur une coupe transversale : une observation par pays, 110 lignes. Nous disposons maintenant d'un panel : 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.

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

Dans nos données, 81 % de la variance du revenu par tête est inter et 19 % 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 : 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 : 14 % seulement de sa variance est inter, et 86 % 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 %, 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 : la coupe ne cherchait à expliquer que ces 14 %, et laissait de côté les 86 % 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 : \(\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.

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 : 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 % par an.

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

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

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 % 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 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 : 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.

Passons de l'illustration au modèle. La transformation intra de la section précédente retranchait à chaque pays sa moyenne temporelle ; 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}\) : 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 ? Regarder l'équation \(r_{i,k} = g_{i,k}\log A_i(0) + v_{i,k+1}\) pour ce qu'elle est : 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}\) :

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

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

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 % par an. C'est plus de quatre fois la vitesse estimée en coupe, 0,86 %, et l'ordre de grandeur que rapporte Islam (1995). La contrainte du modèle de Solow complet, \(\theta = 1-\alpha\), est rejetée : la statistique de Fisher vaut 12,6.

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

À 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 % 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) : plus ils sont bas, meilleure est la spécification.

La régression habituelle ajuste nettement mieux : sa SCR est plus basse de 11 %, 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 les limites du modèle statistique y reviendra.

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

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

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

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.

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

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 : 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 ; 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 : 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.

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

La courbe est plate : 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.

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

L'inconnue de pays apparaît maintenant avec le même coefficient, un, à toutes les périodes ; 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) : 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 : 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 :

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

est la covariance des \(v\), à \(\sigma^2\) près. Une seule entrée échappe à cette formule : 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 :

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é

que l'on maximise numériquement.

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

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 : 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 : \(\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 section sur les limites établira que l'agrégation en est bien la cause.

Soumettons donc le nouvel estimateur à la même épreuve.

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

Il est biaisé lui aussi, et vers le bas, d'un facteur qui n'a rien de constant : 7 % en \(\theta = 0{,}10\), 27 % à hauteur de l'estimation observée, 39 % en \(\theta = 0{,}30\), puis 19 % en \(\theta = 0{,}60\). Mais, à la différence du précédent, il est monotone, et nettement croissant : à 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 : 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 ; on a seulement appris à la lire. Notre estimateur est la balance, et les poids connus sont des panels simulés : 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

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 : on retient la valeur \(\theta^{\star}\) qui aurait produit, en moyenne, l'estimation obtenue sur les vraies données,

Graphiquement, on part de \(\hat\theta_{\text{obs}} = 0{,}176\) sur l'axe vertical, on rejoint horizontalement la courbe, et on lit l'abscisse.

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

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\) : 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 : à \(\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 ; 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 : 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 : 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 % 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.

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

La correction donne \(\theta^{\star} = 0{,}36\) et \(\alpha^{\star} = 0{,}32\), soit une vitesse de convergence médiane de 2,9 % par an : 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 : 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 ? 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}\),

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

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

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 : 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 % de la vitesse médiane va de 0,7 à 5,2 % par an. C'est large, et il faut le dire ainsi : 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\) : 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.

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

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 ; +0,016] dans 90 % 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 : 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 Ce qu'un panel ajoute à une coupe l'a rendue pressante : y a-t-il seulement convergence ? 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.

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

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 % 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 ; 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 % 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 : 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.

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

Le panel livre une estimation de \(\log A_i(0)\) pour chacun des 110 pays : 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 : \(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,

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 : les écarts de productivité globale se trouvent doublés lorsqu'on les exprime en efficacité du travail.

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

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 : les PWT situent le Japon de 1963 à 27 % de l'efficacité américaine, notre estimation à 119 % ; l'Inde à 5 % contre 20 %.

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 : 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é : la pente monte alors à 0,17. Enfin, en écartant les 5 % 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 ; 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\) : 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.

Nous pouvons maintenant revenir aux questions que 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 ? 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 ? 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 ? 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 % et non de 5 %.

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

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

Il y en a une, et elle est solide.

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 : 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 : 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 : sa probabilité critique vaut 0,0002 et son intervalle de confiance à 95 %, 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.

On aurait pu régresser l'efficacité sur les deux taux, mais cela n'aurait rien appris de plus : 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.

Les deux autres questions demandent une dernière expérience.

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

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 % 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 : 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 : 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 : 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 ; le panel en mesure la dose.

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.

Corriger l'endogénéité a déplacé les deux paramètres, et dans le sens attendu : 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 : \(\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é : 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 % 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 %, tombe au milieu de la fourchette de 2 à 4 % que Bond, Hoeffler et Temple (2001) obtiennent par une tout autre voie, les instruments.

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 : 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 : 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 ? En deux étapes. La première donne la perturbation de l'équation : 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}\) :

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,

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 : 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 : 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 ; 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 : 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 : 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.

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

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 : −0,070 contre une bande de [−0,115 ; −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 ; −0,036] : 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 ; 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.

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

Le verdict est net. Reconstruit par prélèvement ponctuel, le modèle ne produit aucune corrélation positive ; 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é.

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

\(\hat\psi\) tombe à 0,02, et le profil montre que ce zéro est bien déterminé : 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.

Faut-il alors enrichir le processus d'erreur ? La question se pose, et elle s'estime : il suffit d'autoriser des coefficients de moyenne mobile à d'autres retards, ce que la covariance \(V\) accepte sans changer une ligne du reste.

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)

Trois enseignements. D'abord, le troisième retard est bien réel : 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 : 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 : 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 ? 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 : 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 % 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 De la coupe au panel écartait un peu vite : elles coûtent \(K-1\) paramètres, mais elles absorbent ces chocs communs, et la comparaison d'ajustement de la section sur 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 % en Australie à 95,9 % au Venezuela, pour une médiane de 13,5 %. Le rapport du maximum à la médiane vaut 7,1, quand le modèle ne le fait dépasser 1,8 que dans 5 % des échantillons. Un seul \(\sigma\) pour 110 pays est une fiction commode.

Le graphique laisse deviner une régularité, et elle se vérifie : 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 : doubler le revenu par tête réduit l'écart-type des innovations de 8 %. Ce n'est toutefois pas un dégradé régulier. Les médianes de \(\hat\sigma_i\) par quartile de revenu valent 14,7 %, 15,0 %, 15,2 % et 8,6 % : 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 % 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.

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

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

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 : 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 ; −0,32] et [−0,19 ; 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 : 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 :

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 : si \(\varepsilon\) les viole, \(z\) les viole aussi. Un contrôle sur \(z\) qui ne décèle rien n'innocente pas \(\varepsilon\) ; 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 : 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 ? 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.

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}\) ; on peut au contraire n'imposer rien du tout et estimer une vitesse par pays, sur les données annuelles.

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

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 ; pour les autres, la vitesse médiane est de 2,7 % par an, la moyenne de 5,0 % 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 : 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 : il ne reste que 3 pays non convergents, la vitesse médiane monte à 8,4 % 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 : en laissant libres, pays par pays, le niveau et la tendance, ils obtiennent une vitesse de l'ordre de 30 % 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 : plus on libère de paramètres par pays, plus la vitesse estimée monte, et moins elle est précise.

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 : 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 : des vitesses estimables, et un test.

Le panel tient l'essentiel de ce qu'il promettait. La note en coupe s'achevait sur un doute qu'elle ne pouvait lever : son rejet du modèle de Solow venait-il du modèle, ou de l'endogénéité du revenu initial ? La réponse est maintenant double. L'endogénéité est réelle et mesurable : 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 : une fois l'endogénéité éliminée, la part du capital tombe à 0,32 et la vitesse monte à 2,9 % 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 : 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 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é : 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.

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