Comment réaliser un test de causalité de Granger dans R

La question de la causalité constitue depuis des siècles l’un des défis intellectuels et méthodologiques les plus ardus de la philosophie des sciences et de l’analyse empirique. De la conception classique formalisée par David Hume, qui réduisait la causalité perçue à une conjonction constante d’événements observables, jusqu’aux modélisations structurelles contemporaines développées par Judea Pearl, la recherche scientifique s’est constamment heurtée à l’impossibilité de déduire formellement un lien de causalité ontologique à partir de simples corrélations observationnelles. Dans l’analyse des séries temporelles, cette problématique est décuplée par la présence d’inerties dynamiques, de tendances partagées et de rétroactions réciproques qui trompent fréquemment l’observateur non averti.

Pour dépasser l’impasse épistémologique consistant à confondre simultanéité et causalité active, l’économiste et statisticien britannique Clive Granger a proposé en 1969 une définition purement opérationnelle et probabiliste de la dépendance temporelle. Fondée sur le postulat d’asymétrie temporelle introduit initialement par Norbert Wiener, cette approche redéfinit la causalité non pas comme un mécanisme physique ou métaphysique sous-jacent, mais comme une capacité prédictive incrémentale : une variable temporelle cause une autre variable si la connaissance du passé de la première permet d’améliorer significativement la prévision de la seconde, par rapport à une prévision fondée exclusivement sur le passé de cette dernière. Ce concept, couronné par le prix de la Banque de Suède en sciences économiques en mémoire d’Alfred Nobel en 2003, est devenu le test de causalité de Granger, un outil d’investigation statistique incontournable.

Dans l’écosystème du logiciel statistique R, la mise en œuvre de ce test s’inscrit au sein d’une chaîne de traitement rigoureuse qui dépasse la simple exécution d’une ligne de commande. Entre l’évaluation formelle de la stationnarité, le calibrage minutieux des retards temporels via des critères d’information stochastiques, la spécification de modèles autorégressifs bivariés et la conduite de diagnostics post-estimation portant sur les résidus, l’analyste doit faire preuve d’une grande rigueur méthodologique. Ce guide exhaustif propose une démarche pas à pas pour maîtriser l’ensemble des fondements théoriques, des formulations mathématiques et des applications pratiques du test de causalité au sens de Granger sous R, en articulant économétrie fondamentale et modélisation empirique des données dynamiques.

1. Fondements théoriques et épistémologiques de la causalité au sens de Granger

1.1 Définition opérationnelle de la prédictibilité temporelle

La conceptualisation de la causalité introduite par Clive W. J. Granger repose de manière fondamentale sur la flèche du temps et le principe d’antériorité temporelle. Dans l’univers physique et empirique tel que nous l’observons, une cause doit impérativement précéder son effet dans le déroulement chronologique des événements. Granger a transposé cette intuition philosophique fondamentale dans le domaine des processus stochastiques en affirmant qu’un événement passé peut contenir une information utile pour anticiper un événement futur, tandis que le futur ne saurait modifier le passé. Ce principe interdit de considérer comme causal tout mécanisme opérant de manière rétrograde et impose une séparation séquentielle stricte entre la variable prédictive et la variable cible.

Il importe d’insister sur la distinction fondamentale et irréductible qui sépare la causalité ontologique de la prédictibilité statistique au sens de Granger. La causalité ontologique suppose un transfert d’énergie, une force mécanique ou un mécanisme substantiel direct reliant deux phénomènes, assurant qu’une manipulation directe de la cause produira un changement de l’effet dans le cadre d’un modèle contrefactuel. À l’inverse, la causalité au sens de Granger est une propriété strictement informationnelle et statistique. Elle indique simplement si l’inclusion de l’historique d’une série chronologique réduit la variance de l’erreur de prévision d’une autre série au-delà de ce que permet l’historique de cette dernière. Ainsi, deux variables peuvent présenter une causalité de Granger significative sans qu’aucun lien physique direct n’existe entre elles, notamment lorsqu’elles sont toutes deux gouvernées par une dynamique sous-jacente commune non observée.

Sur le plan mathématique, Granger a formalisé cette intuition à travers le concept d’asymétrie de l’information disponible dans l’univers stochastique. Soit un ensemble d’informations universel disponible à l’instant t, noté Omega indice t, et soit un sous-ensemble d’informations contenant tout l’historique de l’univers à l’exclusion des valeurs passées et présentes d’une série temporelle spécifique notée x. Si la distribution conditionnelle de la variable y à l’horizon t+1 sachant Omega indice t est différente de la distribution conditionnelle de y à l’instant t+1 sachant Omega privé de l’historique de x, alors x contient une information unique qui permet de mieux cerner le comportement futur de y. En pratique empirique, l’ensemble universel d’informations étant inaccessible, les chercheurs réduisent cet univers informationnel à un ensemble bivarié restreint aux trajectoires temporelles observables de x et de y, introduisant l’hypothèse opérationnelle selon laquelle l’asymétrie temporelle observée dans cet espace clos reflète adéquatement la dynamique étudiée.

1.2 Pertinence méthodologique en sciences du comportement et psychologie

Bien que le test de causalité au sens de Granger trouve son origine historique dans l’économétrie quantitative et l’analyse macroéconomique de séries comme le produit intérieur brut et la masse monétaire, son cadre méthodologique a profondément transformé les sciences comportementales contemporaines. L’avènement des méthodologies d’évaluation écologique momentanée (Ecological Momentary Assessment, ou EMA) et du recueil continu de données physiologiques par capteurs connectés permet aujourd’hui d’enregistrer des séries temporelles intensives au niveau individuel. Les chercheurs en psychologie clinique, cognitive et sociale ne se contentent plus de mesurer des corrélations transversales figées ; ils étudient désormais la dynamique intra-individuelle des processus psychologiques au fil des minutes, des heures et des journées.

Une application clinique emblématique réside dans l’étude des interactions dynamiques entre les biomarqueurs du stress physiologique et les états émotionnels subjectifs. Par exemple, des mesures répétées de la variabilité de la fréquence cardiaque (VFC) associées à des auto-évaluations régulières des affects négatifs permettent de tester empiriquement l’antériorité temporelle de la dérégulation autonome sur la détérioration de l’humeur. En appliquant un test de Granger bivarié, le chercheur peut déterminer si une altération du tonus vagal à l’instant t-1 prédit de façon incrémentale un pic d’anxiété ressenti à l’instant t, ou si, à l’inverse, l’activation cognitive de ruminations émotionnelles précède et déclenche la chute des paramètres physiologiques. Cette modélisation permet de départager des modèles théoriques concurrents en psychopathologie computationnelle.

De même, l’analyse temporelle séquentielle s’applique avec une remarquable acuité à l’exploration des interactions patient-thérapeute au cours d’interventions psychothérapeutiques longitudinales. En quantifiant à intervalles réguliers l’alliance thérapeutique perçue par le patient et l’expression de résistances ou de progrès symptomatiques rapportés au fil des séances, le test de Granger offre un cadre rigoureux pour évaluer si l’affermissement du lien collaboratif précède chronologiquement la réduction des symptômes ou si l’amélioration clinique précoce du patient est le moteur principal qui consolide rétroactivement l’alliance de travail. L’identification de telles dynamiques directionnelles est cruciale pour optimiser la chronologie des interventions thérapeutiques et affiner les protocoles de soin personnalisés.

1.3 Formalisation des hypothèses nulle et alternative

La formulation statistique du test de causalité de Granger repose sur la comparaison explicite de deux structures d’hypothèses imbriquées au sein d’un modèle d’estimation linéaire. L’hypothèse nulle, traditionnellement désignée par la notation H0, stipule l’absence totale de valeur prédictive incrémentale de la variable x sur la variable dépendante y. Sous cette hypothèse nulle, la prise en compte des observations antérieures de x n’apporte aucune réduction statistiquement significative de l’erreur quadratique moyenne de prévision de y, une fois que l’on a déjà pris en considération l’ensemble des valeurs passées de y elle-même. En termes paramétriques, cela signifie que tous les coefficients associés aux décalages temporels de x dans l’équation de régression sont conjointement nuls.

À l’inverse, l’hypothèse alternative, notée H1, postule l’existence d’une causalité au sens de Granger allant de x vers y. Selon cette hypothèse, au moins un des coefficients de régression associés aux retards temporels de la variable candidate x est significativement différent de zéro, ce qui se traduit empiriquement par un gain mesurable de variance expliquée sur la série cible y. La modélisation autorégressive permet ainsi de tester si la dynamique temporelle passée de x exerce une impulsion déterministe ou stochastique sur la trajectoire contemporaine de y. Il convient de souligner que cette hypothèse alternative n’implique pas que tous les retards soient simultanément actifs, mais que leur contribution conjointe s’écarte de la simple fluctuation aléatoire.

La formalisation du protocole de décision exige la fixation préalable d’un seuil de significativité statistique alpha, généralement établi à 0,05 ou 0,01 dans les standards de la recherche académique. Ce seuil quantifie la probabilité maximale acceptable de commettre une erreur de type I, c’est-à-dire de rejeter à tort l’hypothèse nulle et de conclure à tort à l’existence d’une causalité de Granger alors que les variables sont en réalité prédictivement indépendantes. Dans le contexte des séries temporelles fortement autocorrelées, le risque d’erreur de type I peut être massivement gonflé si les propriétés fondamentales du processus sous-jacent, telles que la stationnarité et l’absence d’autocorrélation des résidus, ne sont pas scrupuleusement respectées par l’expérimentateur.

2. Cadre mathématique et modélisation autorégressive

2.1 Spécification du modèle autorégressif univarié restreint

L’assise quantitative du test repose sur la confrontation directe de deux spécifications économétriques emboîtées. La première spécification correspond au modèle univarié restreint, dans lequel la variable cible y est modélisée exclusivement à partir de sa propre histoire temporelle. Pour un horizon d’observation donné et un ordre de retard maximal p, l’équation autorégressive s’écrit de la manière suivante : la valeur de y à l’instant t est égale à une constante déterministe notée alpha, à laquelle s’ajoute la somme pondérée de ses p valeurs passées, pondérées par des coefficients autorégressifs réels gamma indice i, complétée par un terme d’erreur stochastique résiduel noté epsilon restreint à l’instant t.

L’estimation des paramètres de cette équation autorégressive est effectuée par la méthode des moindres carrés ordinaires (MCO). Sous les hypothèses classiques de stationnarité du processus univarié, les estimateurs des coefficients gamma sont convergents et asymptotiquement gaussiens. L’objectif fondamental de cette étape d’ajustement consiste à capturer l’ensemble de l’inertie propre de la série temporelle cible. En épurant le signal y de son auto-dépendance structurelle, le modèle restreint établit une ligne de base de prédiction, servant de référence neutre pour quantifier la variance non expliquée par l’historique intrinsèque de la variable dépendante.

À l’issue de l’estimation par moindres carrés ordinaires du modèle restreint, on extrait la somme des carrés des résidus, désignée par l’acronyme RSSR (pour Restricted Residual Sum of Squares). Cette somme est calculée en sommant les carrés de tous les termes d’erreur résiduels du modèle restreint sur l’ensemble de la période d’estimation effective, qui compte T observations après neutralisation des p premières valeurs initiales nécessaires au calage des retards. La variance résiduelle associée, calculée en divisant la somme des carrés des résidus par les degrés de liberté résiduels correspondants (T moins p moins 1), représente l’erreur quadratique incompressible de y lorsque seule sa propre trajectoire historique est mobilisée pour l’anticipation stochastique.

2.2 Spécification du modèle autorégressif bivarié non restreint

La seconde étape analytique introduit le modèle complet, dit non restreint, au sein duquel l’équation autorégressive univariée est enrichie par l’adjonction systématique des décalages temporels de la seconde variable stochastique x. L’équation de régression bivariée postule ainsi que la valeur de y à l’instant t est fonction linéaire de sa propre constante d’interception, de la somme des p retards de y pondérés par des coefficients gamma indice i, mais également de la somme des p retards de x pondérés par un ensemble de coefficients croisés bêta indice j, à laquelle s’ajoute un nouveau terme d’erreur stochastique résiduel noté epsilon non restreint à l’instant t.

Dans ce modèle complet, la structure matricielle des coefficients bêta indice j représente directement l’impact prédictif marginal de chaque état antérieur de x sur la réalisation présente de y. Les moindres carrés ordinaires sont de nouveau mobilisés pour résoudre simultanément le système d’équations normales et estimer l’ensemble des paramètres du modèle non restreint. Cette opération permet de réajuster simultanément les composantes autorégressives propres et les dépendances croisées, garantissant que tout effet statistique attribué aux décalages temporels de x n’est pas simplement un artefact d’une covariance linéaire non contrôlée avec le passé récent de y.

De l’estimation de cette spécification bivariée non restreinte découle une seconde somme des carrés des résidus, désignée par l’acronyme RSSU (pour Unrestricted Residual Sum of Squares). Par construction géométrique des projections orthogonales associées à la méthode des moindres carrés ordinaires, l’intégration de prédicteurs supplémentaires ne peut jamais accroître la somme des carrés des erreurs résiduelles ; ainsi, RSSU est nécessairement inférieure ou égale à RSSR. L’enjeu probabiliste du test réside donc exclusivement dans l’évaluation de l’amplitude de cette réduction : la diminution de la somme des résidus au carré est-elle suffisamment volumineuse pour être attribuée à un apport d’information substantiel plutôt qu’à l’ajustement mécanique inhérent à l’augmentation de la dimensionnalité paramétrique du modèle ?

2.3 La statistique de test F et distribution d’échantillonnage

Pour arbitrer entre le modèle restreint et le modèle non restreint, la procédure statistique mobilise le ratio de Fisher-Snedecor, formellement configuré pour tester des restrictions linéaires conjointes sur un vecteur de paramètres. La formule de la statistique de test F met directement en concurrence la variation relative des sommes des carrés résiduels et s’énonce comme le rapport entre deux quantités normalisées par leurs degrés de liberté respectifs :

F = [ (RSSR – RSSU) / p ] / [ RSSU / (T – 2p – 1) ]

Dans cette formulation algébrique rigoureuse, le numérateur quantifie la réduction de la somme des carrés des résidus permise par l’inclusion des retards de la variable x, divisée par le nombre de restrictions linéaires imposées par l’hypothèse nulle, qui correspond exactement au nombre de retards p. Le dénominateur correspond à la variance résiduelle non biaisée du modèle non restreint, calculée en divisant la somme des carrés des résidus du modèle complet (RSSU) par ses degrés de liberté d’échantillonnage, soit la taille d’échantillon effective T minorée du nombre total de paramètres estimés dans l’équation bivariée (une constante, p paramètres autorégressifs propres et p paramètres croisés, soit 2p + 1 coefficients au total).

Sous la validité de l’hypothèse nulle H0 stipulant que tous les coefficients bêta sont conjointement égaux à zéro, et sous l’hypothèse de normalité et d’indépendance sérielle des résidus, cette statistique d’échantillonnage suit exactement une loi de Fisher-Snedecor paramétrée par p degrés de liberté au numérateur et T – 2p – 1 degrés de liberté au dénominateur. Lorsque la taille de l’échantillon T devient très grande, le produit de la statistique F par le nombre de degrés de liberté du numérateur p converge asymptotiquement vers une loi du khi-deux à p degrés de liberté, ce qui correspond à la formulation asymptotique du test de Wald. La p-valeur associée à la valeur F observée représente la probabilité d’obtenir une réduction au moins aussi importante de l’erreur de prédiction sous l’effet du seul hasard d’échantillonnage.

3. Prérequis fondamentaux pour l’analyse des séries temporelles

3.1 La condition primordiale de stationnarité faible

L’application valide et légitime du test de causalité au sens de Granger repose sur une exigence statistique absolue et non négociable : les séries temporelles soumises à l’analyse doivent obligatoirement satisfaire à la condition de stationnarité faible, également appelée stationnarité au sens large ou stationnarité d’ordre deux. Un processus stochastique est qualifié de faiblement stationnaire si et seulement si ses propriétés statistiques fondamentales demeurent strictement invariantes au cours du temps. Cela implique, en premier lieu, que l’espérance mathématique de la série ne présente aucune dérive systématique et reste rigoureusement constante à chaque point de la trajectoire temporelle observée.

En second lieu, la condition de stationnarité exige que la variance de la variable demeure finie et parfaitement constante à travers les périodes, excluant tout phénomène d’explosion variance-temps ou de régimes de volatilité instables. Enfin, et il s’agit d’une dimension critique pour les processus autorégressifs, l’autocovariance entre deux observations mesurées à des instants temporels distincts ne doit dépendre exclusivement que de l’écart chronologique ou décalage temporel qui les sépare, et non du point absolu dans le temps où cette covariance est calculée. Si l’un de ces trois critères est transgressé, le processus est qualifié de non stationnaire, ce qui invalide l’ensemble des théorèmes asymptotiques sous-tendant les tests d’inférence linéaire.

L’omission de cette vérification de stationnarité expose l’analyste au piège redoutable de la régression fallacieuse, un phénomène magistralement mis en évidence par les travaux de George Udny Yule en 1926 et réactualisé de manière formelle par Clive Granger et Paul Newbold en 1974. En présence de deux séries temporelles indépendantes présentant chacune une marche aléatoire ou une tendance déterministe partagée, l’ajustement par moindres carrés ordinaires engendre artificiellement des coefficients de détermination extrêmement élevés et des statistiques de Student ou de Fisher démesurément significatives. Dans ce cas de figure pathologique, le test de Granger conclura avec une certitude statistique trompeuse à l’existence d’une relation de causalité directionnelle, alors que les deux signaux ne partagent qu’une non-stationnarité commune sans le moindre couplage informationnel authentique.

3.2 Homoscédasticité et indépendance sérielle des perturbations

Au-delà de la stationnarité globale des processus observés, la mécanique inférentielle du test de causalité impose des contraintes structurales sévères sur le comportement des termes d’erreur résiduels du modèle autorégressif bivarié. L’hypothèse fondamentale requiert que les résidus epsilon se comportent comme un bruit blanc gaussien pur. Cette caractéristique implique que l’espérance mathématique conditionnelle des résidus est rigoureusement nulle, que leur distribution présente une variance constante au fil du temps (homoscédasticité) et, de façon primordiale, qu’aucune dépendance temporelle ne persiste entre les résidus d’époques différentes (absence d’autocorrélation sérielle).

Lorsque la structure de retard p retenue dans la modélisation autorégressive est insuffisante pour épurer la totalité de la mémoire du système stochastique, les perturbations résiduelles conservent une corrélation sérielle interne. Cette autocorrélation des résidus contrevient frontalement aux hypothèses du théorème de Gauss-Markov. Si les estimateurs des coefficients de régression demeurent formellement sans biais, la matrice de variance-covariance des coefficients estimés devient quant à elle complètement erronée et généralement sous-estimée. Il en résulte que le calcul de la statistique F de Fisher-Snedecor devient mathématiquement caduc, provoquant une distorsion majeure des probabilités critiques et une inflation drastique du risque de fausse découverte.

De surcroît, la survenue d’hétéroscédasticité conditionnelle, phénomène caractérisé par une fluctuation temporelle structurée de la variance des chocs où des périodes de calme succèdent à des phases de forte volatilité, altère profondément l’exactitude du test standard. Dans de tels contextes, la distribution exacte du ratio F s’écarte des tables théoriques conventionnelles de Snedecor. Il devient alors impératif d’identifier la présence potentielle d’effets de type autorégressif à hétéroscédasticité conditionnelle (effets ARCH) sur les séries de résidus et d’appliquer, le cas échéant, des corrections d’erreurs-types hétéroscédastiques et autocorrélées (matrices robustes de type Newey-West) pour préserver la validité de l’inférence causale.

3.3 Linéarité de la relation dynamique sous-jacente

La formulation conventionnelle du test de Granger repose sur le postulat d’une relation de dépendance strictement linéaire entre les observations contemporaines de la variable cible et les réalisations décalées dans le temps des variables explicatives. Le modèle postule que l’amplitude de la réponse de y à une variation unitaire passée de x demeure strictement proportionnelle et identique, quelle que soit la position du système dans son espace d’état et quelle que soit la magnitude ou la polarité positive ou négative du choc initialement induit par x. Cette restriction mathématique assure une tractabilité analytique optimale par le recours aux moindres carrés ordinaires.

Néanmoins, ce cadre linéaire constitue une approximation simplificatrice qui atteint rapidement ses limites épistémologiques lorsqu’il est confronté aux dynamiques comportementales, physiologiques ou écologiques complexes. Dans les systèmes vivants et psychologiques, les relations fonctionnelles manifestent fréquemment des effets de seuil, des saturations sigmoïdes, de l’hystérèse ou des bifurcations de régimes qualitatifs. Par exemple, une élévation modérée du niveau de stress peut faciliter la performance attentionnelle, tandis qu’un dépassement d’un seuil critique entraîne un effondrement brutal des capacités cognitives. Face à de tels processus non linéaires, le test de causalité linéaire de Granger s’avère structurellement aveugle : il peut échouer à rejeter l’hypothèse nulle d’indépendance prédictive alors même qu’un lien déterministe non linéaire puissant régit la dynamique conjointe des séries.

Enfin, le modèle autorégressif standard est caractérisé par une vulnérabilité aiguë face aux observations aberrantes isolées ainsi qu’aux ruptures structurelles temporelles. L’apparition d’un choc exogène soudain ou d’une transition de régime au sein de l’environnement expérimental engendre des perturbations locales démesurées qui pèsent de manière disproportionnée sur le calcul de la somme des carrés des résidus. Une unique rupture structurelle non modélisée dans la moyenne ou la variance d’une série chronologique peut suffire à créer l’illusion statistique d’une causalité de Granger ou, à l’inverse, masquer complètement une relation prédictive stable en faussant les estimations des moindres carrés.

4. Configuration de l’écosystème R et bibliothèques requises

4.1 Présentation du package central lmtest

Dans l’environnement informatique R, l’infrastructure dédiée à l’évaluation formelle des hypothèses linéaires en économétrie temporelle s’articule autour du package fondamental lmtest (développé par Torsten Hothorn, Achim Zeileis, Richard W. Farebrother et Christian Cummins). Cette bibliothèque constitue la référence méthodologique standard pour la réalisation de tests diagnostiques appliqués aux modèles de régression linéaire, offrant un ensemble d’outils statistiques optimisés pour interroger la validité des spécifications empiriques face aux violations classiques des hypothèses de modélisation.

Au cœur du package lmtest figure la fonction native grangertest(), spécifiquement conçue pour exécuter avec une efficacité computationnelle maximale le test bivarié d’antériorité prédictive selon le formalisme de Fisher-Snedecor. Cette routine encapsule de manière transparente l’ajustement du modèle autorégressif bivarié complet, l’estimation conjointe du modèle univarié restreint, l’extraction systématique des résidus d’estimation et le calcul précis du ratio F associé à sa probabilité critique sous l’hypothèse nulle d’absence de causalité informationnelle.

L’un des atouts déterminants de l’architecture logicielle de grangertest() réside dans son interopérabilité universelle avec les structures de données séquentielles de R. La fonction est capable de traiter sans altération structurelle aussi bien les objets chronologiques de base de type ts que les structures temporelles avancées manipulées via les bibliothèques zoo (spécialisée dans les séries ordonnées régulières et irrégulières) et xts (structure de séries chronologiques extensibles). Cette flexibilité facilite grandement l’incorporation de l’évaluation causale au sein de pipelines de traitement de données préexistants.

4.2 Outils complémentaires pour la manipulation et la modélisation

Bien que la fonction issue du package lmtest fournisse l’outil direct d’exécution du test, la conduite méthodique d’une recherche causale complète nécessite la mobilisation d’un ensemble cohérent de packages spécialisés. Pour l’évaluation préliminaire et incontournable des propriétés de stationnarité des données, la bibliothèque tseries s’avère indispensable. Elle met à la disposition de l’utilisateur les routines de test de racine unitaire, en particulier l’implémentation standardisée du test augmenté de Dickey-Fuller et du test de stationnarité KPSS.

Pour dépasser le cadre bivarié simplificateur et aborder l’analyse dynamique multivariée, l’intégration du package vars s’impose comme une nécessité technique incontournable. Conçu par Bernhard Pfaff, ce module constitue le socle le plus puissant de modélisation vectorielle autorégressive sous R. Il permet non seulement de sélectionner automatiquement la profondeur optimale des retards temporels à l’aide de critères d’information stochastiques multivariés, mais offre également les outils nécessaires pour estimer des systèmes VAR complets, calculer des fonctions de réponse impulsionnelle et exécuter des décompositions de variance de l’erreur de prévision.

Enfin, la manipulation ergonomique, le nettoyage rigoureux et l’exploration graphique des séries temporelles modernes tirent un profit considérable de la suite tidyverse, incluant notamment dplyr, tidyr et le moteur de visualisation graphique ggplot2. L’association de ces paradigmes modernes de manipulation fonctionnelle avec des bibliothèques de séries temporelles spécialisées comme tibbletime ou timetk garantit une traçabilité totale du flux de travail statistique, de la phase de transformation des variables brutes jusqu’à la production de figures de qualité éditoriale destinées aux publications scientifiques.

4.3 Script d’installation et de chargement reproductible

La reproductibilité computationnelle exige l’élaboration d’un script d’amorce transparent et déterministe, garantissant que l’ensemble des modules requis est identifié, installé si nécessaire et chargé dans l’espace de noms de la session de travail sans déclencher d’incohérence logicielle. Dans l’écosystème R, une bonne pratique de programmation consiste à définir un vecteur de chaînes de caractères énumérant les packages essentiels, puis à interroger le registre des bibliothèques locales pour n’exécuter l’installation que pour les dépendances absentes de l’environnement de l’utilisateur.

L’exécution séquentielle des commandes R s’articule autour de la syntaxe suivante : on recense les bibliothèques indispensables que sont lmtest, tseries, vars et tidyverse. À l’aide de la primitive installed.packages(), le script isole les modules manquants et déclenche leur acquisition via la fonction install.packages(), en pointant explicitement vers un miroir sécurisé du réseau CRAN (Comprehensive R Archive Network). Une fois l’intégrité de l’environnement validée, chaque bibliothèque est successivement appelée dans l’espace de travail au moyen de l’instruction canonique library().

Pour prévenir les collisions d’espaces de noms, fréquentes lors de l’utilisation simultanée de modules de séries temporelles et du tidyverse, il est particulièrement conseillé d’utiliser les préfixes explicites d’appels de fonction lors de l’accès aux routines critiques. L’analyste averti prendra soin de documenter la version exacte du moteur d’exécution de R ainsi que les numéros de version de chaque package mobilisé en consignant l’état global du système grâce à la fonction sessionInfo(). Cette démarche constitue le fondement d’une démarche de recherche transparente, garantissant que les calculs de causalité pourront être répliqués à l’identique à travers différentes plateformes informatiques et versions logicielles.

5. Préparation, nettoyage et structuration des séries temporelles

5.1 Traitement des observations manquantes et espacement temporel

Avant d’engager toute modélisation autorégressive sous R, il est impératif d’évaluer la régularité structurelle de l’échantillonnage temporel des données recueillies. Le formalisme du test de causalité de Granger repose de manière axiomatique sur l’hypothèse d’un intervalle temporel strictement constant séparant chaque point de mesure successif (delta t constant). Dans les études observationnelles, en psychologie écologique momentanée ou lors de recueil de biomarqueurs, les ruptures de protocole, les non-réponses accidentelles ou les pannes techniques génèrent fréquemment des trous d’échantillonnage et des discontinuités temporelles.

La présence d’observations manquantes (valeurs notées NA sous R) brise la contrainte de retard temporel régulier. L’application brute d’un décalage temporel d’ordre 1 sur une série lacunaire associerait arbitrairement deux mesures séparées en réalité par plusieurs heures ou plusieurs jours, introduisant un biais majeur dans l’estimation des coefficients de mémoire autorégressive. L’analyste doit donc impérativement vérifier l’absence d’intervalles temporels erratiques. Si des valeurs manquantes isolées parsèment la série, des méthodes d’imputation adaptées aux séries temporelles, telles que l’interpolation par splines cubiques, le lissage par modèles d’espace d’état ou l’algorithme du filtre de Kalman disponibles dans le package imputeTS, doivent être mises en œuvre avec prudence.

Néanmoins, l’imputation artificielle de séries chronologiques comporte des dangers statistiques considérables. Toute procédure d’interpolation mathématique génère par nature une forme de lissage artificiel qui atténue la variance intrinsèque des signaux et induit une autocorrélation de substitution au sein des résidus d’échantillonnage. Si la proportion de données manquantes dépasse un seuil critique, généralement estimé à 5 ou 10 pour cent de l’échantillon global, les estimations autorégressives deviennent tributaires des hypothèses algorithmiques de l’imputation plutôt que de la dynamique réelle des processus étudiés. Dans une telle situation, il convient d’éviter l’interpolation et de privilégier des découpages par fenêtres stationnaires continues ou des modèles autorégressifs en temps continu adaptés.

5.2 Transformation des objets de données au format ts

Pour que les routines de modélisation chronologique puissent interpréter correctement les dépendances temporelles, les données brutes, classiquement importées sous forme de data frame ou de structure tibble rectangulaire, doivent être formellement encapsulées au sein d’objets temporels natifs. Dans le langage R de base, la fonction canonique ts() est dédiée à cette conversion typologique. Elle associe au vecteur de mesures quantitatives un vecteur d’attributs temporels indispensables, comprenant l’indice temporel initial, l’indice final et la fréquence d’échantillonnage calendaire.

Le paramétrage explicite de la fonction ts() s’effectue en renseignant minutieusement l’argument numérique start et l’argument d’échelle frequency. Si les observations sont recueillies sur un rythme journalier avec une dynamique hebdomadaire récurrente, la fréquence temporelle sera paramétrée à 7. Pour des collectes horaires répétées sur une journée complète, elle prendra la valeur de 24. L’initialisation précise du calendrier via un vecteur numérique à deux dimensions spécifiant l’unité principale et la subdivision temporelle garantit que l’alignement des signaux bivariés sera respecté de manière rigoureuse lors de la construction des retards.

Une fois l’objet de classe ts instancié dans la session R, il est fondamental de procéder à une inspection diagnostique de ses attributs métrologiques en utilisant les primitives start(), end(), frequency() et deltat(). Ces commandes permettent de certifier que l’infrastructure temporelle de l’objet est intacte, qu’aucun décalage d’échantillonnage n’est apparu lors de la phase de conversion typologique et que le vecteur de données est parfaitement aligné pour accueillir sans ambiguïté les opérations d’estimation de décalages croisés.

5.3 Visualisation exploratoire des trajectoires temporelles

L’exploration visuelle constitue l’étape préliminaire incontournable avant tout déploiement de formalismes quantitatifs. L’œil humain et l’analyse graphique permettent de détecter instantanément des anomalies structurelles profondes que les tests numériques globaux peuvent occulter. Sous R, la mobilisation combinée des instructions plot.ts() ou l’usage perfectionné de ggplot2 (via l’extension de formatage de données temporelles tidyr::pivot_longer()) permet de tracer simultanément les trajectoires longitudinales des deux variables explorées sur un repère orthonormé partagé.

L’examen minutieux de ces graphiques chronologiques vise en premier lieu à détecter l’existence manifeste de tendances déterministes. Une série temporelle qui s’élève ou décline continuellement sur l’ensemble de la fenêtre d’observation ne possède pas une moyenne constante et viole donc d’emblée la condition de stationnarité faible. De même, l’inspection permet de déceler la présence de composantes saisonnières cycliques, de patterns périodiques ou d’oscillations synchronisées susceptibles de générer des covariances mécaniques illusoires si elles ne sont pas désaisonnalisées au préalable.

Enfin, le tracé chronologique bivarié sert à repérer d’éventuels changements de régime structurel, tels que des transitions brutales de niveau moyen, des décalages d’amplitude ou des chocs ponctuels hors normes. Dans les sciences du comportement, de telles discontinuités traduisent souvent des événements contextuels majeurs vécus par les sujets (prise médicamenteuse soudaine, choc émotionnel aigu, modification de l’environnement de vie). L’identification graphique de ces ruptures incite le modélisateur à délimiter des sous-périodes d’analyse homogènes ou à planifier l’introduction de variables indicatrices muettes (variables muettes d’intervention) au sein du modèle économétrique.

6. Évaluation formelle de la stationnarité sous R

6.1 Application du test augmenté de Dickey-Fuller (ADF)

L’inspection visuelle des trajectoires temporelles doit impérativement être corroborée par des tests d’hypothèse formels destinés à trancher la question de la présence ou de l’absence d’une racine unitaire dans le polynôme autorégressif des séries étudiées. Le test le plus largement déployé en économétrie appliquée est le test augmenté de Dickey-Fuller (ADF test). Sous R, sa mise en œuvre standardisée s’effectue aisément par l’appel à la fonction adf.test() hébergée au sein de la bibliothèque tseries.

Le principe mathématique du test ADF consiste à régresser les variations premières de la variable sur son niveau décalé d’une période, en incorporant une constante, éventuellement une tendance linéaire, ainsi qu’un nombre suffisant de variations décalées pour épurer complètement l’autocorrélation des résidus. L’hypothèse nulle fondamentale H0 de ce test stipule que le processus stochastique possède une racine unitaire, ce qui équivaut à affirmer que la série est non stationnaire (processus intégré d’ordre 1, noté I(1)). L’hypothèse alternative H1 correspond à la stationnarité stricte du processus autour d’une moyenne déterministe ou d’une tendance linéaire.

La règle de décision statistique repose sur l’analyse rigoureuse de la p-valeur générée dans la console de sortie de R. Si la probabilité critique rapportée par adf.test() est inférieure au seuil de risque choisi de 0,05, l’expérimentateur est fondé à rejeter formellement l’hypothèse nulle de racine unitaire, concluant ainsi que la série ne présente pas de non-stationnarité stochastique. En revanche, si la p-valeur s’avère supérieure à 0,05, le non-rejet de l’hypothèse nulle impose d’admettre la présence d’une non-stationnarité qui interdit l’application immédiate du test de causalité de Granger dans les niveaux bruts des données.

6.2 Test complémentaire de stationnarité KPSS

En matière d’inférence temporelle, s’en remettre exclusivement au test ADF expose l’analyste à un biais méthodologique bien documenté : la faible puissance statistique du test de Dickey-Fuller en présence de petits échantillons ou face à des racines autorégressives proches du cercle unité (coefficients autorégressifs proches de 0,95 ou 0,99). Pour remédier à cette vulnérabilité, la démarche économétrique moderne préconise une approche de validation croisée par l’adjonction systématique du test de Kwiatkowski, Phillips, Schmidt et Shin (test KPSS), accessible via la commande kpss.test() du package tseries.

La distinction conceptuelle majeure du test KPSS réside dans l’inversion délibérée de la structure des hypothèses stochastiques. Dans le cadre KPSS, l’hypothèse nulle H0 postule formellement que la série temporelle est stationnaire (autour d’un niveau moyen ou d’une tendance déterministe), tandis que l’hypothèse alternative H1 stipule la présence d’une racine unitaire non stationnaire. Ce renversement de perspective modifie substantiellement le fardeau de la preuve statistique, permettant de vérifier si les données apportent une évidence empirique suffisante pour rejeter formellement l’hypothèse de stationnarité.

L’exploitation conjointe des deux procédures de test débouche sur une grille de décision robuste :

• Si le test ADF rejette la racine unitaire (p < 0,05) et que le test KPSS ne rejette pas la stationnarité (p > 0,05), la stationnarité de la série est confirmée de façon concordante.
• Si le test ADF échoue à rejeter la racine unitaire (p > 0,05) et que le test KPSS rejette l’hypothèse de stationnarité (p < 0,05), la présence d’une non-stationnarité stochastique est indiscutablement établie.
• Si les deux tests aboutissent à des conclusions contradictoires, cela indique que l’information contenue dans les données est insuffisante ou parasitée par des anomalies dynamiques (présence de racines fractionnaires ou de ruptures de tendance), justifiant une prudence maximale et le recours systématique à des transformations stabilisatrices.

6.3 Procédures de différenciation pour stationnariser les séries

Lorsque le protocole de vérification croisée ADF-KPSS établit formellement le caractère non stationnaire d’une ou des deux séries chronologiques sous revue, une transformation mathématique préalable s’avère incontournable afin de rétablir les conditions d’application de la théorie asymptotique des moindres carrés. L’outil privilégié de stationnarisation face aux processus possédant une racine unitaire stochastique consiste à appliquer l’opérateur de différenciation première, matérialisé sous l’environnement R par la fonction native diff().

La différenciation première consiste à soustraire à chaque observation temporelle la valeur prise par la série à l’instant immédiatement antérieur. Cette opération filtre efficacement la mémoire persistante d’une marche aléatoire et neutralise toute tendance linéaire déterministe présente dans le niveau du signal. Si une série brute nécessite une seule différenciation pour que son comportement devienne stationnaire au sens large, elle est classée dans la catégorie des processus intégrés d’ordre un, désignée par la nomenclature I(1). La série transformée issue de diff() est quant à elle qualifiée d’intégrée d’ordre zéro, notée I(0).

Après application de la transformation par l’opérateur diff(), il est impératif d’exécuter de nouveau les routines de diagnostic adf.test() et kpss.test() sur les séries différenciées afin de certifier formellement leur retour à la stationnarité stochastique. Seules ces séries épurées de toute racine unitaire pourront alimenter légitimement le test de causalité au sens de Granger. Il est à noter que l’interprétation substantielle des résultats portera alors non plus sur le niveau absolu des variables d’intérêt, mais sur leurs incréments ou leurs taux de variation séquentiels d’une période à l’autre.

7. Sélection optimale du nombre de retards (order)

7.1 Importance théorique du calibrage des décalages temporels

Le calibrage du paramètre désigné par la variable order ou retard p dans la modélisation autorégressive représente la décision méthodologique la plus stratégique et la plus délicate qui incombe à l’analyste. Ce paramètre quantifie l’horizon temporel sur lequel la mémoire passée du système est interrogée pour prédire les réalisations futures. Un choix inadéquat de cette profondeur temporelle entraîne des dégradations sévères de la validité inférentielle du test de causalité.

D’une part, la sous-paramétrisation du retard, qui consiste à retenir une valeur de p inférieure à la durée réelle de l’effet dynamique sous-jacent, induit un biais de variable omise. En tronquant arbitrairement la mémoire stochastique du processus, le modèle autorégressif échoue à absorber la totalité de l’autocorrélation interne du signal. Cette mémoire résiduelle non modélisée se déverse directement dans les termes d’erreur, provoquant une autocorrélation sérielle des résidus qui détruit la validité de la statistique F de Fisher et majore artificiellement le risque de détection de causalités fallacieuses.

D’autre part, la sur-paramétrisation, qui se caractérise par la sélection d’un nombre de retards temporellement excessif, pénalise lourdement la parcimonie et la précision statistique du modèle. L’introduction de chaque retard supplémentaire requiert l’estimation d’un coefficient supplémentaire pour y et d’un coefficient supplémentaire pour x, consommant deux degrés de liberté au sein de l’échantillon fini disponible. Cette dispersion inutile de la puissance d’échantillonnage gonfle mécaniquement la variance des estimateurs par moindres carrés ordinaires, diluant la significativité des coefficients croisés pertinents et augmentant considérablement le risque de commettre une erreur de type II (incapacité à détecter une causalité pourtant bien réelle au niveau populationnel).

7.2 Exploitation des critères d’information (AIC, BIC, HQ)

Pour dépasser l’arbitraire d’un choix subjectif de la profondeur de décalage, la théorie statistique moderne s’appuie sur l’optimisation mathématique de critères d’information stochastiques. Ces métriques évaluent le compromis quantitatif entre la qualité d’ajustement empirique du modèle aux données observées (mesurée par la maximisation de la vraisemblance ou la minimisation de la somme des carrés des résidus) et une pénalité formelle indexée sur la complexité paramétrique du modèle (le nombre total de coefficients estimés).

Le critère le plus largement documenté est le critère d’information d’Akaike (AIC). Fondé sur la minimisation de la divergence de Kullback-Leibler, l’AIC applique une pénalité modérée proportionnelle à deux fois le nombre de paramètres estimés. En conséquence, l’AIC tend à privilégier la performance de prédiction pure et tolère des structures de retards relativement amples, au risque d’une légère sur-paramétrisation dans les contextes d’échantillons volumineux. À l’opposé, le critère d’information bayésien de Schwarz (BIC ou SC) applique une pénalisation logarithmique rigoureuse indexée sur la taille d’échantillon T (pénalité proportionnelle au produit du nombre de paramètres par le logarithme népérien de T). Le critère BIC est asymptotiquement convergent et favorise des modèles extrêmement parcimonieux, convenant particulièrement aux échantillons larges.

Entre ces deux approches figure le critère de Hannan-Quinn (HQ), dont la pénalité s’établit selon une formule logarithmique double, offrant un compromis robuste entre l’optimisme relatif de l’AIC et l’austérité contraignante du BIC. Dans les échantillons de taille modeste fréquemment rencontrés dans la recherche comportementale ou les études longitudinales intensives, il est classique que ces critères proposent des ordres de retards légèrement divergents. L’analyste rigoureux arbitre généralement en faveur de la structure la plus parcimonieuse (guidée par le BIC) dès lors que les résidus associés sont exempts d’autocorrélation, ou adopte une approche d’analyse de sensibilité explorant la stabilité du test de Granger sous les différents ordres recommandés.

7.3 Utilisation pratique de la fonction VARselect() du package vars

Dans l’environnement computationnel de R, l’identification automatisée de la structure optimale de retard pour un système multivarié ou bivarié est grandement facilitée par la commande VARselect(), intégrée au package vars. Cette fonction élimine la contrainte de programmer manuellement une boucle d’ajustement de régressions multiples imbriquées et standardise le calcul simultané des différents critères d’information sur une plage continue de décalages temporels.

La mise en œuvre pratique débute par l’assemblage matriciel des deux séries chronologiques stationnarisées, réalisé au moyen de l’instruction cbind(). Cet objet matriciel bivarié est ensuite passé en premier argument à la routine VARselect(). L’argument crucial lag.max doit être paramétré pour fixer la borne supérieure de l’espace de recherche des retards, par exemple un maximum de 5 ou 10 périodes, en veillant à ce que cette valeur demeure cohérente avec la taille globale de l’échantillon disponible et la fréquence d’échantillonnage théorique du phénomène mesuré.

L’exécution de cette fonction renvoie une structure de données contenant deux éléments d’une grande clarté diagnostique : une matrice détaillée affichant les valeurs numériques prises par les critères AIC, HQ, SC (BIC) et FPE (erreur de prédiction finale) pour chaque valeur de décalage explorée, et un vecteur synthétique désigné par $selection qui isole directement le retard optimal correspondant au minimum global de chaque critère d’information. Dès lors, le chercheur dispose d’une justification algorithmique objective et reproductible pour sélectionner la valeur d’ordre p qui sera injectée dans le test de Granger formel.

8. Exécution pas à pas du test avec la fonction grangertest()

8.1 Syntaxe et anatomie des paramètres de grangertest()

La fonction grangertest() issue de la bibliothèque lmtest constitue le vecteur d’exécution le plus épuré et le plus polyvalent pour tester l’absence de causalité au sens de Granger sous R. Sa syntaxe générale s’appuie sur le modèle universel de formulation des relations sous R, utilisant l’opérateur tilde pour découpler la variable dépendante (située à gauche) des prédicteurs (situés à droite). L’expression standard grangertest(y ~ x, order = p, data = df) formalise l’interrogation suivante : les valeurs passées de la série x améliorent-elles la prévision de la série y au-delà de sa propre dynamique, en considérant un horizon de p décalages temporels ?

Alternativement, la fonction accepte une syntaxe directe sous forme de passage de deux arguments vectoriels indépendants : grangertest(x, y, order = p). Il convient d’attirer l’attention sur une inversion syntaxique fréquente qui constitue une source majeure d’erreur d’interprétation chez les praticiens : dans la syntaxe par formule y ~ x, la variable y est la cible (l’effet potentiel) et x est la source causale postulée. En revanche, dans la syntaxe à arguments positionnels grangertest(x, y), le premier argument spécifié x correspond à la variable source et le second argument y représente la variable cible dont la prévisibilité est interrogée. Une vigilance absolue s’impose donc lors de la rédaction du script pour préserver la directionnalité voulue du test.

L’argument numérique order reçoit l’entier scalaire p déterminé préalablement par l’arbitrage des critères d’information stochastiques issus de VARselect(). Enfin, la fonction propose un argument de traitement des valeurs manquantes, désigné par na.action, paramétré par défaut sur na.omit. Ce paramètre exclut mécaniquement les premières lignes du jeu de données qui comportent des valeurs indéterminées en raison de la génération des décalages temporels, garantissant que l’estimation des deux modèles (restreint et non restreint) s’opère sur un sous-ensemble strict et rigoureusement superposable d’observations temporelles synchronisées.

8.2 Exemple concret de code R : du chargement à l’exécution

Pour illustrer de manière totalement transparente et déterministe le déroulement de la séquence computationnelle complète, considérons la génération d’un système dynamique simulé au sein duquel une dépendance causale unidirectionnelle d’ordre 2 est explicitement injectée. L’initialisation du générateur congruentiel de nombres pseudo-aléatoires via l’instruction set.seed(123) garantit la reproductibilité intégrale des données simulées à travers toute machine hôte exécutant R.

Considérons la séquence d’instructions suivante, décrite ici de manière séquentielle et rigoureuse :

En premier lieu, l’expérimentateur génère un vecteur d’innovations stochastiques gaussiennes indépendantes pour une variable x sur une longueur de 300 pas de temps. La variable x est configurée pour suivre un processus autorégressif univarié stable d’ordre 1, où la valeur contemporaine de x dépend à hauteur de 0,6 de sa valeur immédiatement antérieure, augmentée d’un terme d’aléa pur. En second lieu, une variable y est générée en définissant sa valeur à l’instant t comme une fonction linéaire de ses propres retards (auto-dépendance de 0,4 au retard 1), à laquelle s’ajoute de manière causale l’influence explicite de la variable x au retard 2 (avec un coefficient de transmission de 0,55), augmentée de son propre terme de bruit blanc indépendant.

Les deux séries stationnaires ainsi générées sont regroupées au sein d’une matrice ou d’un data frame standard. L’analyste invoque ensuite les packages lmtest et vars. Après avoir exécuté VARselect() pour confirmer que l’horizon d’ordre 2 minimise les critères d’information conjoints, l’analyste lance la commande canonique grangertest(y ~ x, order = 2). Cette ligne de code orchestre l’ensemble des opérations sous-jacentes : alignement des décalages temporels, partitionnement matriciel, inversion des équations de moindres carrés ordinaires pour les modèles avec et sans la variable x, et calcul de la statistique de Fisher finale, le tout s’exécutant en une fraction de seconde dans la mémoire vive de l’ordinateur.

8.3 Lecture technique de la sortie console de R

L’exécution de la commande grangertest() génère instantanément dans la console R un tableau synthétique d’analyse de variance (ANOVA) standardisé, dont la lecture requiert une compréhension précise de l’architecture de modélisation déployée par la routine. La console affiche un en-tête mentionnant formellement le test de causalité au sens de Granger (Granger causality test), rappelant explicitement la formule testée ainsi que l’ordre de retard p mobilisé pour l’ajustement.

Le corps du résultat se structure sous la forme d’un tableau comparatif à deux lignes, confrontant directement le Modèle 1 (Model 1) et le Modèle 2 (Model 2) :

• Le Modèle 1 représente la spécification bivariée non restreinte, intégrant l’ensemble des retards propres de la variable cible y ainsi que les retards croisés de la variable prédictive x.
• Le Modèle 2 correspond au modèle univarié restreint, épuré de tous les retards de x, conservant exclusivement les composantes autorégressives de y.

Pour chacune de ces deux spécifications, le tableau rapporte deux grandeurs métrologiques fondamentales :

• La colonne étiquetée Res.Df indique les degrés de liberté résiduels (Residual Degrees of Freedom), calculés comme la taille d’échantillon T tronquée des décalages initiaux, diminuée du nombre de coefficients de régression estimés dans le modèle correspondant. On observe systématiquement que le Modèle 2 possède p degrés de liberté supplémentaires par rapport au Modèle 1, ce qui traduit directement l’omission des p paramètres associés à la variable x.
• La colonne intitulée Df met en évidence cette différence stricte de paramétrisation, matérialisée par la valeur négative ou positive de p.
• La colonne majeure F consigne la valeur empirique du ratio de Fisher-Snedecor, issue de la confrontation des sommes des carrés des résidus des deux équations.
• Enfin, la colonne finale désignée par Pr(>F) présente la p-valeur exacte associée à la statistique F observée, assortie de codes de significativité visuels standardisés (symbolisés par des astérisques). C’est cette valeur critique qui scelle la conclusion du test statistique.

9. Interprétation statistique et biométrique des résultats

9.1 Décision statistique face au seuil de significativité

La formulation de la conclusion statistique repose sur l’application rigoureuse du cadre formel d’inférence de Neyman-Pearson. L’analyste confronte la p-valeur (colonne Pr(>F)) au seuil de significativité nominal alpha préalablement fixé (classiquement 0,05). Si la p-valeur observée est inférieure ou égale à 0,05, la probabilité d’observer un tel gain de prédictibilité sous l’hypothèse de pure indépendance stochastique est infinitésimale. Le statisticien prend alors la décision formelle de rejeter l’hypothèse nulle H0 au profit de l’hypothèse alternative H1.

Le rejet de l’hypothèse nulle autorise à énoncer que la variable candidate x « cause au sens de Granger » la variable dépendante y au seuil de risque considéré. Cette décision atteste que la contribution informationnelle des p valeurs passées de x surpasse les simples fluctuations aléatoires propres à l’échantillonnage, fournissant un avantage systématique pour réduire la variance des résidus d’anticipation de la série cible.

À l’inverse, si la probabilité critique rapportée est strictement supérieure à 0,05, la décision méthodologique s’énonce en termes d’échec au rejet de l’hypothèse nulle. Il est épistémologiquement crucial de rappeler que l’échec au rejet de H0 ne constitue en aucun cas une démonstration formelle de l’absence totale d’influence de x sur y. Conformément à l’adage fondamental de Carl Sagan, « l’absence de preuve n’est pas la preuve de l’absence ». Cet échec au rejet signifie simplement que, dans le cadre strict du modèle linéaire autorégressif testé, au niveau de retard p sélectionné et compte tenu de la puissance de l’échantillon fini disponible, les données ne permettent pas de mettre en évidence un gain de prévision statistiquement concluant.

9.2 Traduction substantielle dans le domaine d’application

Une fois la sentence statistique arrêtée, l’analyste doit impérativement traduire ce résultat dans le champ conceptuel de la discipline scientifique dans laquelle s’ancre la recherche. Cette contextualisation substantielle exige une discipline sémantique sans faille. Il est impératif d’éviter le raccourci intellectuel consistant à proclamer sans nuance que « x cause y ». La formulation scientifiquement exacte impose d’affirmer que « les trajectoires antérieures de la variable x procurent une valeur prédictive incrémentale statistiquement significative quant aux réalisations futures de la variable y, au-delà de l’information déjà encapsulée dans l’histoire propre de y ».

Par ailleurs, la démarche d’interprétation moderne ne saurait se cantonner à la seule évaluation de la p-valeur, dont la petitesse peut être un simple sous-produit mécanique d’une taille d’échantillon massive. L’analyste méticuleux évalue systématiquement la taille de l’effet stochastique sous-jacent. Cette évaluation s’appuie sur la quantification du gain relatif de variance résiduelle, mesuré par la variation proportionnelle de la somme des carrés des résidus entre le modèle restreint et le modèle complet, assimilable à un coefficient de détermination partiel (R-deux partiel). Un test de Granger peut afficher une significativité extrême (p < 0,001) dans un échantillon de 50 000 observations tout en ne réduisant l’erreur de prévision que d’une fraction marginale de 0,2 %, ce qui indique une portée pratique ou clinique quasi négligeable du phénomène prédictif.

Cette distinction est particulièrement prégnante en psychologie et en psychiatrie computationnelle. Mettre en évidence que des oscillations du niveau de cortisol salivaire causent au sens de Granger des fluctuations de l’humeur dépressive ne prouve aucunement que le cortisol constitue la cause moléculaire unique ou dominante du syndrome dépressif. Ce résultat atteste simplement d’un couplage temporel directionnel observable au sein du flux de données, constituant un indice empirique solide pour orienter des recherches expérimentales plus interventionnistes (par exemple des protocoles d’administration pharmacologique contrôlée avec randomisation contrefactuelle).

9.3 Rédaction académique des résultats aux normes professionnelles

La communication des résultats au sein d’articles soumis à des revues académiques à comité de lecture doit se conformer aux standards de transparence et de rigueur formelle, notamment ceux édictés par l’American Psychological Association (normes APA, 7e édition). La rédaction doit expliciter sans équivoque l’ensemble des paramètres métrologiques ayant gouverné l’inférence : l’indication explicite de la stationnarisation préalable des séries, l’ordre de retard retenu (p), la statistique F exacte, les degrés de liberté résiduels des deux modèles comparés, la valeur précise de la p-valeur ainsi qu’une mesure de l’ampleur de l’effet.

Un modèle de rédaction académique standardisé peut être formulé de la façon suivante :

« Afin d’examiner la directionnalité temporelle liant la variabilité de la fréquence cardiaque (VFC) et le niveau d’anxiété subjective auto-rapporté, un test de causalité au sens de Granger a été appliqué sur les séries chronologiques stationnarisées en différences premières. L’ordre optimal de décalage temporel a été préalablement calibré à p = 2 sur la base de la minimisation conjointe des critères d’information d’Akaike (AIC) et bayésien de Schwarz (BIC). L’analyse de variance comparant le modèle autorégressif univarié restreint et le modèle bivarié complet révèle un apport prédictif incrémental significatif de la VFC passée sur la trajectoire contemporaine de l’anxiété, F(2, 293) = 8,42, p < 0,001, delta RSS = 5,44 %. Les fluctuations antérieures du tonus vagal permettent ainsi d’anticiper de manière statistiquement significative l’aggravation ou l’atténuation subséquente des affects anxieux, confirmant l’hypothèse d’une antériorité temporelle de la régulation autonome sur l’expérience émotionnelle consciente. »

Une telle formalisation textuelle procure à la communauté scientifique une visibilité totale sur la chaîne de traitement inférentielle, permettant une réplication immédiate des constats empiriques et une intégration rigoureuse des données au sein de futures méta-analyses.

10. Analyse bidirectionnelle et modélisation vectorielle autorégressive (VAR)

10.1 Test de rétroaction et causalité réciproque

L’une des limites fondamentales de l’exécution isolée de la commande grangertest(y ~ x) réside dans son caractère asymétrique unilatéral. Dans la vaste majorité des systèmes naturels, biologiques, économiques et psychologiques, les variables n’interagissent pas au sein de structures causales à sens unique. Il est donc méthodologiquement impératif d’exécuter de façon systématique le test réciproque symétrique, en évaluant dans un second temps si la variable y cause également la variable x au sens de Granger via l’instruction grangertest(x ~ y, order = p).

Cette démarche bidirectionnelle exhaustive débouche sur une typologie quadruple des configurations stochastiques :

• Indépendance prédictive complète : Ni x ne cause y, ni y ne cause x au sens de Granger. Les séries évoluent selon des dynamiques temporelles orthogonales ou découplées.
• Causalité unidirectionnelle directe : La variable x cause la variable y, sans que l’influence inverse ne présente de significativité statistique. L’asymétrie temporelle est totale.
• Causalité unidirectionnelle rétrograde : La variable y cause la variable x, tandis que la trajectoire de x est prédictivement stérile pour anticiper y.
• Boucle de rétroaction dynamique (feedback loop) : Les deux tests conjoints rejettent simultanément l’hypothèse nulle d’absence de causalité. Les variables x et y se prédisent mutuellement à travers le temps.

L’identification de boucles de rétroaction temporelle revêt une importance conceptuelle majeure dans le champ de la psychopathologie computationnelle et de la théorie des réseaux de symptômes. De nombreux troubles mentaux, tels que le cercle vicieux de l’insomnie et de l’anxiété, s’auto-entretiennent précisément sous la forme de causalités réciproques circulaires. Dans de telles architectures dynamiques, l’analyse unilatérale est réductrice, et seule la caractérisation formelle de la double causalité permet de capturer l’essence homéostatique ou déstabilisatrice du système clinique observé.

10.2 Transition vers le modèle VAR global avec la fonction VAR()

Lorsque la présence de causalités réciproques ou l’ambition de modéliser plus de deux variables temporelles simultanément est avérée, l’approche séquentielle par équations isolées atteint ses limites. Il devient indispensable de basculer vers le formalisme unifié de la modélisation vectorielle autorégressive (modèle VAR). Sous R, ce saut qualitatif est orchestré par la fonction VAR() de la bibliothèque vars, qui permet d’estimer simultanément l’intégralité d’un système dynamique multivarié.

La fonction VAR() prend en entrée la matrice regroupant l’ensemble des séries chronologiques stationnarisées et ajuste de manière concomitante un système d’équations autorégressives où chaque variable du vecteur devient successivement la cible de ses propres retards et des retards de l’ensemble des autres composantes du système. La commande synthétique mon_modele_var <- VAR(donnees_stationnaires, p = p_opt, type = "const") estime ainsi par moindres carrés ordinaires multivariés l’intégralité de la matrice de transition dynamique reliant les états temporels du système.

Une fois le modèle vectoriel ajusté, le package vars permet d’interroger la causalité de Granger multivariée globale grâce à la méthode spécialisée causality(mon_modele_var, cause = "x"). Cette routine exécute un test de Wald multivarié évaluant simultanément si l’ensemble des retards de la variable indicée comme cause peut être exclu du bloc d’équations prédisant l’ensemble des autres variables du système. Ce passage à l’échelle supérieure dote le chercheur d’un cadre analytique d’une grande puissance pour sonder les interactions dynamiques au sein de systèmes complexes hautement multidimensionnels.

10.3 Fonctions de réponse impulsionnelle (IRF) et décomposition de variance

Si le test de causalité de Granger indique avec précision l’existence d’une influence prédictive significative et en identifie la directionnalité, il reste silencieux sur deux questions empiriques majeures : quel est le signe qualitatif (positif ou négatif) de cette influence, et quel est son profil d’amortissement ou de propagation à travers le temps ? Pour répondre à ces interrogations, la modélisation VAR s’enrichit de deux outils post-estimatoires essentiels : les fonctions de réponse impulsionnelle (Impulse Response Functions, IRF) et la décomposition de la variance de l’erreur de prévision (FEVD).

Sous le package vars, la commande irf(mon_modele_var, impulse = "x", response = "y", n.ahead = 10, boot = TRUE) simule l’impact dynamique qu’engendre l’application d’un choc stochastique exogène unitaire (une impulsion équivalente à un écart-type) sur la variable x à l’instant présent sur la trajectoire future de la variable y au cours des 10 périodes temporelles suivantes. L’option de rééchantillonnage par bootstrap (boot = TRUE) permet de construire des intervalles de confiance empiriques à 95 % autour de la trajectoire de réponse. La visualisation graphique résultante (via la méthode plot() appliquée à l’objet IRF) permet d’observer visuellement si le choc sur x induit une augmentation transitoire, une chute dépressive immédiate ou une oscillation amortie de la variable cible y, éclairant la dynamique intime de l’effet causal.

En complément, la décomposition de variance de l’erreur de prévision, mobilisée via l’instruction fevd(mon_modele_var, n.ahead = 10), calcule le pourcentage de la variance résiduelle de prévision de y imputable aux chocs structurels spécifiques de x comparativement aux chocs propres à y, pour différents horizons temporels futurs. Ce tableau de pourcentages offre une mesure directe et intuitive de l’importance relative de la variable causale au sein du système stochastique, montrant par exemple que si x n’explique que 2 % de la variance de y à court terme (horizon de 1 pas de temps), sa part contributive monte à 35 % à moyen terme (horizon de 6 pas de temps), soulignant un effet prédictif à latence d’accumulation différée.

11. Diagnostics post-estimation et vérification de la robustesse

11.1 Diagnostic de l’autocorrélation résiduelle

L’estimation d’un modèle autorégressif bivarié ou multivarié ne peut être considérée comme valide sur le plan inférentiel tant que ses termes d’erreurs résiduels n’ont pas subi une batterie complète d’examens diagnostiques post-estimation. Le diagnostic le plus critique porte sur l’indépendance sérielle des résidus. Si les erreurs d’ajustement manifestent une corrélation temporelle interne, les degrés de liberté réels sont violés et les valeurs p associées au test de Granger deviennent artificiellement optimistes.

Au sein du package lmtest, le diagnostic d’autocorrélation s’opère au moyen de la routine de Breusch-Godfrey (Breusch-Godfrey test), accessible via la commande bgtest(). Cette méthode surpasse le classique test de Durbin-Watson car elle demeure rigoureusement exacte en présence de variables dépendantes décalées dans le membre de droite de l’équation de régression, ce qui caractérise précisément les modèles autorégressifs. L’hypothèse nulle H0 du test de Breusch-Godfrey stipule l’absence totale d’autocorrélation résiduelle sérielle jusqu’à un ordre de décalage donné.

Si la p-valeur générée par bgtest() est inférieure au seuil conventionnel de 0,05, l’expérimentateur rejette l’hypothèse d’indépendance sérielle et doit en tirer les conséquences méthodologiques immédiates : la structure de retard retenue (l’ordre p) est insuffisante pour épurer la dynamique temporelle du processus. La solution économétrique s’impose d’elle-même : il est nécessaire d’incrémenter le paramètre order d’une ou plusieurs unités et de réestimer le test de Granger jusqu’à ce que les résidus du modèle complet satisfassent pleinement au critère de bruit blanc séquentiel.

11.2 Détection de l’hétéroscédasticité conditionnelle (effets ARCH)

Un deuxième angle mort diagnostique fréquemment négligé concerne l’homoscédasticité de la variance des résidus d’échantillonnage. Dans de nombreuses séries chronologiques, en particulier celles mesurant des réactivités physiologiques ou des tensions affectives, les résidus présentent fréquemment des grappes de volatilité où la variance instantanée des chocs dépend de l’amplitude des chocs passés. Cette structure caractérise la présence d’hétéroscédasticité autorégressive conditionnelle (effets ARCH).

La détection de ces effets sous R s’exécute à l’aide de la fonction arch.test() présente dans le package vars, qui réalise un test du multiplicateur de Lagrange (test ARCH-LM) sur les résidus du modèle vectoriel autorégressif. L’hypothèse nulle H0 stipule que les carrés des résidus sont sériellement indépendants, traduisant une variance conditionnelle rigoureusement constante. Si l’évaluation statistique débouche sur une p-valeur significative (p < 0,05), l’homoscédasticité classique est réfutée, signalant que la distribution théorique du ratio F est déformée par l’instabilité de la variance.

Face à une confirmation d’effets ARCH au sein des résidus, l’analyste doit impérativement consolider son inférence causale. Dans le cadre univarié sous lmtest, il est préconisé de recalculer la statistique de Wald du test de Granger en appliquant une matrice de variance-covariance robuste aux hétéroscédasticités et aux autocorrélations sérielles (matrice HAC de type Newey-West ou Andrews), mobilisable via la fonction coeftest() combinée avec le package sandwich. Cette correction ajuste à la hausse les erreurs-types des estimateurs, prévenant l’analyste contre le risque de faux positifs causaux induits par les bouffées de variabilité non modélisées.

11.3 Stabilité structurelle des paramètres du modèle

L’estimation autorégressive repose implicitement sur l’hypothèse forte de stationnarité structurelle : les coefficients de dépendance dynamique (les paramètres gamma et bêta) sont supposés constants et immuables d’un bout à l’autre de la période temporelle couverte par l’échantillon. Si un événement exogène, une réorganisation institutionnelle ou un virage comportemental modifie en profondeur la dynamique relationnelle liant les séries au cours de l’étude, l’estimation globale d’un unique jeu de paramètres devient invalide et produit une moyenne spatio-temporelle artificielle dépourvue de sens physique.

Pour certifier la permanence de la relation causale à travers le temps, l’environnement R offre les outils de la bibliothèque spécialisée strucchange (développée par Achim Zeileis et collègues). L’une des méthodes les plus robustes repose sur le test des sommes cumulées des résidus récursifs, universellement désigné sous l’acronyme de test CUSUM. Accessible via la commande efp() (pour Empirical Fluctuation Process), cette approche suit le cumul séquentiel des erreurs de prévision normalisées au fur et à mesure que l’échantillon s’accroît dans le temps.

L’exploitation graphique de l’objet CUSUM via l’instruction plot() trace la dérive empirique des résidus au sein d’un couloir délimité par des frontières critiques correspondant au seuil alpha de 5 %. Tant que la trajectoire des résidus cumulés oscille à l’intérieur de ces bornes statistiques, l’hypothèse de stabilité structurelle des coefficients du modèle autorégressif ne peut être rejetée. À l’inverse, toute sortie de trajectoire hors du faisceau critique atteste de la survenue d’une cassure structurelle majeure, contraignant l’analyste à identifier la date de rupture et à segmenter son test de Granger en sous-périodes distinctes pour documenter la façon dont le lien de causalité temporelle est apparu, a muté ou s’est éteint au fil du temps.

12. Pièges méthodologiques courants et recommandations en recherche

12.1 Le problème de la variable tierce omise (fausse causalité)

Le risque épistémologique le plus dévastateur qui pèse sur l’interprétation d’un test de causalité au sens de Granger bivarié réside dans le problème de la variable tierce omise, fréquemment désigné sous le terme de causalité fallacieuse par facteur de confusion dynamique (spurious causality). Dans un système bivarié restreint aux seules séries x et y, la présence d’une prédictibilité temporelle apparente de x sur y peut être entièrement provoquée par l’action asynchrone d’une troisième variable sous-jacente non mesurée, notée z, qui pilote simultanément les deux signaux avec des constantes de temps différentes.

Considérons une illustration comportementale et chronobiologique éloquente : supposons que l’on étudie à très haute fréquence temporelle le niveau d’éveil cognitif ou attentionnel (variable y) et la température corporelle centrale (variable x). Un test de Granger bivarié révélera indubitablement que les variations thermiques corporelles devancent chronologiquement et causent au sens de Granger les variations de l’attention. Pourtant, la température ne constitue aucunement le moteur causal direct de la cognition ; c’est le rythme circadien sous-jacent piloté par l’horloge biologique suprachiasmatique (variable tierce z) qui synchronise les deux paramètres, en modulant la thermorégulation avec une latence plus courte que l’activation du cortex préfrontal.

Pour prémunir l’analyse contre cette fausse causalité systémique, la recommandation méthodologique moderne impose d’élargir autant que possible le cadre d’analyse vers des modélisations vectorielles autorégressives (VAR) multivariées intégrant l’ensemble des covariables contextuelles majeures identifiées par la théorie substantive du domaine. En contrôlant simultanément les antériorités temporelles des facteurs d’ambiance, des cycles circadiens, des médications et des événements environnementaux, le chercheur « conditionne » le test de Granger sur un ensemble informationnel Omega bien plus dense, limitant substantiellement le risque d’ériger une simple avance de phase descriptive en un couplage causal authentique.

12.2 Taux d’échantillonnage et agrégation temporelle inadéquate

La validité des conclusions issues d’un test de causalité de Granger dépend de façon intrinsèque et étroite de l’adéquation existant entre la vitesse cinétique réelle du phénomène investigué et la fréquence d’échantillonnage empirique à laquelle les mesures sont relevées dans le temps. L’asymétrie temporelle selon Granger n’est mesurable que si la discrétisation du temps opérée par le dispositif de mesure est capable de capter la propagation séquentielle de l’onde causale.

Le premier écueil majeur réside dans le biais d’agrégation temporelle. Lorsqu’un processus causal opère à une échelle temporelle fine (par exemple, une réaction neurophysiologique qui se déploie en quelques fractions de seconde ou quelques minutes) alors que le protocole de recueil agrège ou moyenne les observations à une maille trop large (par exemple des agrégats journaliers ou hebdomadaires), la dynamique séquentielle sous-jacente est mathématiquement effacée. À l’échelle macroscopique résultante, la causalité réelle à haute fréquence se dissout et se transforme artificiellement en une corrélation contemporaine instantanée au sein de la même période d’observation, interdisant tout rejet de l’hypothèse nulle du test de Granger classique.

Le second péril, encore plus insidieux, est celui de l’inversion fallacieuse du sens causal par sous-échantillonnage (phénomène de repliement de spectre ou d’aliasing causal). Comme l’ont démontré plusieurs travaux théoriques en économétrie temporelle, échantillonner un système dynamique continu à des intervalles trop distendus peut engendrer une distorsion de phase telle que la série source apparaît mathématiquement en retard sur la série cible dans les données discrétisées, inversant de manière spectaculaire le sens de la causalité de Granger dans la console du logiciel statistique. La définition de l’intervalle de recueil delta t doit par conséquent être méticuleusement fondée sur des considérations théoriques et biomécaniques rigoureuses préalablement à toute campagne d’acquisition de données.

12.3 Synthèse des bonnes pratiques pour une recherche reproductible dans R

Pour garantir la rigueur, l’intégrité et la reproductibilité des analyses de causalité de Granger au sein de la communauté scientifique, l’ensemble de la démarche empirique sous R doit se conformer à un protocole standardisé strict et séquentiel. Ce protocole peut être synthétisé sous la forme d’une chaîne de décision analytique intégrée :

• Vérification de la régularité temporelle : S’assurer que le pas d’échantillonnage est constant, traiter les observations manquantes avec circonspection sans induire de sur-lissage artificiel, et formater rigoureusement les objets chronologiques dans les structures de séries appropriées sous R (ts, zoo, xts).
• Évaluation formelle de la stationnarité : Ne jamais exécuter de test sur des séries en niveaux sans avoir préalablement validé leur stationnarité faible par une approche de validation croisée combinant le test ADF (recherche de racine unitaire) et le test KPSS (recherche de stationnarité). Différencier à l’ordre requis si une non-stationnarité stochastique est attestée.
• Sélection algorithmique des retards : Calibrer le paramètre d’ordre p au moyen des critères d’information stochastiques (AIC, BIC, HQ) via VARselect(), en privilégiant la parcimonie tout en veillant à la couverture de la dynamique théorique attendue.
• Exécution bidirectionnelle symétrique : Conduire systématiquement le test dans les deux directions (x causant y, puis y causant x) pour discriminer les causalités unilatérales, les indépendances stochastiques et les boucles de rétroaction dynamique. Élargir vers des modèles VAR multivariés si des variables de confusion sont identifiées.
• Diagnostics post-estimation exhaustifs : Certifier formellement la pureté des résidus en vérifiant l’absence d’autocorrélation sérielle (test de Breusch-Godfrey), l’absence d’hétéroscédasticité conditionnelle (test ARCH-LM) et la stabilité temporelle structurelle des paramètres (test CUSUM). Mobiliser des corrections de covariance robustes de type HAC si des violations résiduelles mineures sont constatées.
• Reproductibilité computationnelle intégrale : Figer l’aléa de calcul par des graines de reproductibilité explicites (set.seed()), encapsuler les environnements logiciels et les versions de packages via des gestionnaires de dépendances hermétiques tels que renv, et publier l’intégralité du code source annoté parallèlement aux jeux de données pseudonymisés dans des dépôts ouverts pérennes.

En observant scrupuleusement ces impératifs méthodologiques, l’analyste transforme le test de causalité au sens de Granger sous R en un instrument d’investigation d’une remarquable puissance descriptive, capable d’éclairer la dynamique temporelle fine des phénomènes complexes tout en conservant une stricte lucidité épistémologique sur la nature prédictive et non déterministe de la causalité statistique ainsi caractérisée.

Références

• Akaike, H. (1974). A new look at the statistical model identification. IEEE Transactions on Automatic Control, 19(6), 716–723. https://doi.org/10.1109/TAC.1974.1100705
• Breusch, T. S. (1978). Testing for autocorrelation in dynamic linear models. Australian Economic Papers, 17(31), 334–355. https://doi.org/10.1111/j.1467-8454.1978.tb00635.x
• Dickey, D. A., & Fuller, W. A. (1979). Distribution of the estimators for autoregressive time series with a unit root. Journal of the American Statistical Association, 74(366), 427–431. https://doi.org/10.1080/01621459.1979.10482531
• Engle, R. F. (1982). Autoregressive conditional heteroscedasticity with estimates of the variance of United Kingdom inflation. Econometrica, 50(4), 987–1007. https://doi.org/10.2307/1912773
• Godfrey, L. G. (1978). Testing against general autoregressive and moving average error models when the regressors include lagged dependent variables. Econometrica, 46(6), 1293–1301. https://doi.org/10.2307/1913829
• Granger, C. W. J. (1969). Investigating causal relations by econometric models and cross-spectral methods. Econometrica, 37(3), 424–438. https://doi.org/10.2307/1912791
• Granger, C. W. J., & Newbold, P. (1974). Spurious regressions in econometrics. Journal of Econometrics, 2(2), 111–120. https://doi.org/10.1016/0304-4076(74)90034-7
• Hannan, E. J., & Quinn, B. G. (1979). The determination of the order of an autoregression. Journal of the Royal Statistical Society: Series B (Methodological), 41(2), 190–195. https://doi.org/10.1111/j.2517-6161.1979.tb01072.x
• Kwiatkowski, D., Phillips, P. C. B., Schmidt, P., & Shin, Y. (1992). Testing the null hypothesis of stationarity against the alternative of a unit root. Journal of Econometrics, 54(1–3), 159–178. https://doi.org/10.1016/0304-4076(92)90104-Y
• Lütkepohl, H. (2005). New Introduction to Multiple Time Series Analysis. Springer-Verlag. https://doi.org/10.1007/978-3-540-27752-1
• Newey, W. K., & West, K. D. (1987). A simple, positive semi-definite, heteroskedasticity and autocorrelation consistent covariance matrix. Econometrica, 55(3), 703–708. https://doi.org/10.2307/1913610
• Pfaff, B. (2008). Analysis of Integrated and Cointegrated Time Series with R (2nd ed.). Springer. https://doi.org/10.1007/978-0-387-75967-8
• Schwarz, G. (1978). Estimating the dimension of a model. The Annals of Statistics, 6(2), 461–464. https://doi.org/10.1214/aos/1176344136
• Wiener, N. (1956). The theory of prediction. In E. F. Beckenbach (Ed.), Modern Mathematics for Engineers (pp. 165–190). McGraw-Hill.
• Zeileis, A., & Hothorn, T. (2002). Diagnostic checking in regression relationships. R News, 2(3), 7–10. https://cran.r-project.org/doc/Rnews/Rnews_2002-3.pdf
• Zeileis, A., Leisch, F., Hornik, K., & Kleiber, C. (2002). strucchange: An R package for testing for structural change in linear regression models. Journal of Statistical Software, 7(2), 1–38. https://doi.org/10.18637/jss.v007.i02