L’estimation en économétrie spatiale est essentielle pour comprendre les interactions dans les modèles de régression. Cet article explore les défis d’estimation liés aux modèles d’auto-corrélation spatiale et propose une application empirique sur les déterminants du chômage en Tunisie.
L’estimation
Les problèmes d’estimation associés aux modèles de regression spatiale sont distincts de part et d’autre des deux catégories de modèles spatiaux. En effet, les modèles d’auto- corrélation spatiale des erreurs représentent des cas particuliers de spécification avec des termes d’erreurs non-sphériques tandis que l’insertion d’une variable endogène décalée conduit à une forme d’endogénéité.
Ces complications peuvent être abordées moyennant les méthodes économétriques usuelles. Néanmoins, la nature bidirectionnelle (feed-back) des relations spatiales et le recours aux poids spatiaux exigent des techniques dédiées au traitement des problèmes associés à l’inadaptation des estimateurs MCO. Pour trouver des estimateurs convergents et asymptotiquement normaux et efficients, la méthode la plus couramment utilisée est celle du maximum de vraisemblance.
Il est également pos- sible d’avoir recours à la méthode des variables instrumentales ou à celle des moments généralisés.
Dans cette section, une étude de ces approches sera présentée précédée d’un bref examen des propriétés asymptotiques associées à la dimension spatiale.
Propriétés asymptotiques et Espace :
Les propriétés des estimateurs en économétrie spatiale sont fondées sur des consi- dérations asymptotiques qui utilisent des lois des grands nombres et des théorèmes
centraux limites en vue d’établir la convergence et la normalité asymptotique. Selon Anselin[2010], il ne s’agit pas de simples extentions à deux dimensions de résultats des séries temporelles car il existe des différences sur trois aspects.
Premièrement, de nombreux processus spatiaux (par exemple, SAR et SMA) in- duisent une hétéroscedasticité, nécessitant un traitement conjoint de l’autocorrélation et de l’hétérogénéité. Plus précisément, cette non stationnarité exclut le recours aux théorèmes centraux limites relatifs aux processus stationnaires, qui ont été utilisées pour obtenir les propriétés des estimateurs dans les modèles où l’interdépendance spa- tiale est basée sur la distance métrique. Le traitement des processus spatiaux refermant des poids spatiaux est compliqué par la nécessité de considérer des TCL et des LGN spécifiques. Ceci est causé par le fait que les poids dépendent de la taille de l’échan- tillon, empêchant les résultats standards (par exemple, l’estimation du maximum de vraisemblance) d’être directement applicable.
Un deuxième aspect distinctif concerne la façon dont l’échantillonnage spatial est conceptualisée et extrait. Les données spatiales peuvent être considérées comme consti- tuant soit une surface continue ou un ensemble d’objet discret. Le premier point de vue est répandu dans les sciences physiques et les approches géostatistiques. Le second est plus courant dans l’économétrie spatiale.
En effet, lors de l’examen d’une population d’objets spatiaux, le mécanisme par lequel les observations ou les lieux sont choisis n’est pas arbitraire et peut être déterministe ou stochastique. Dans le cas stochastique, l’ensemble d’indices aléatoires utilisés pour prélever des échantillons de la population correspond à un processus directeur (directing process).
En conséquence, le processus spatial aléatoire lui même devient subordonné au processus directeur. Cette complica- tion n’est généralement pas rencontré dans le domaine temporel.
La troisième distinction se rapporte à la manière dont l’échantillon augmente pour
atteindre la limite asymptotique n → ∞. Les propriétés asymptotiques des estima- teurs peuvent être déterminées de deux manières. La première correspond aux modèles d’échantillonnage usuels : c’est-à-dire aux modèles d’échantillonnage extensifs (increa- sing domain asymptotics). L’extension de l’échantillon résulte alors de l’addition de nouvelles observations à celles qui sont déjà disponibles.
Dans le cas de données spa- tiales, un modèle d’échantillonnage extensif conduit à ajouter indéfiniment de nouvelles unités à celles qui ont été sélectionnées antérieurement sans que les unités déjà présentes soient modifiées. Dans le deuxième cas, l’espace étudié est borné et non extensible, le modèle d’échantillonnage intensif (infill asymptotics) est plus approprié.
Ces modèles correspondent à des échantillons dont la taille croît suite à des découpages de plus en plus fins de la zone étudiée : les unités spatiales sont redécoupées et donc modifiées.
Cette procédure donne une surface d’échantillonnage de plus en plus dense, avec une distance minimale entre les lieux d’échantillonnage proche de zéro quand n → ∞. Les applications en économétrie spatiale font davantage appel au premier type de modèle. En plus des restrictions habituelles sur la variance et les moments d’ordre supérieur des variables du modèle, ces conditions se traduisent par des contraintes sur la structure des matrices de poids. Hormis les modèles à spécifications complexes, ces conditions sont
satisfaites en construisant des matrices de contiguïté.
Le maximum de vraisemblance :
L’estimation par le maximum de vraisemblance des modèles de décalage spatial et des modèles d’autocorrélation spatiale sont décrits dans Anselin [1988].
Le modèle de décalage spatial :
Le point de départ de l’estimation par maximum de vraissemblance des modèles de régression spatiale est l’hypothèse de normalité des termes d’erreurs. Afin d’étudie le cas le plus général d’hétéroscédasticité et/ou d’autocorrélation, on confère au vecteur des termes d’erreurs de taille n × 1 une distribution normale multivariée de telle sorte que ε ∼ N (0, Σθ). L’indice θ dénote que Σ pourrait être une fonction d’un vecteur de paramètres θ de taille p × 1. Néanmoins, on admet communémment des erreurs i.i.d. permettant d’écrire ε ∼ N (0, σ2I) avec θ = σ2.
Pour passer de la vraissemblance du vecteur des erreurs à la vraissemblance des variables dépendantes observées, une transformation Jacobienne doit être introdtuite. Dans la cas du modèle de décalage spatial [1.29], ceci correspond à :
|∂ε/∂y| = |∂(y−ρWy−Xβ)/∂y| = |I − ρW | (1.48)
le déterminant d’une matrice de taille n × n. Il est à noter que le même Jaco- bien ce réduit au scalaire 1 dans le cadre d’un modèle de régression standad puisque
|∂(y−Xβ)/∂y| = |I| = 1. La présence du Jacobien constitue une complication majeure au niveau du calcul.
L’utilisation du résultat standard pour une distribution normale multivariée, permet d’écrire la log-vraisemblance pour le modèle décalage spatial comme suit :
L = − (ln 2π) − ln Σ +ln |I − ρW | − 1 (y − ρWy − Xβ)𝘫Σ−1(y − ρWy − Xβ) (1.49)
1 1
2 2 θ 2 θ
La maximisation de la log-vraisemblance n’est pas équivalente à la minimisation des moindres carrés pondérés (le dernier terme de L) relative au modèle de régression
linéaire standard. La principale différence réside dans la présence du terme log-Jacobien
|I − ρW |, montrant implicitement qu’un recours aux moindres carrés pondérés débou- chera sur des estimateurs non-convergeants du fait de l’endogénéité dans le terme Wy. Par ailleurs, la maximisation de la fonction de log-vraisemblance nécessite l’évaluation du jacobien pour chaque nouvelle valeur de ρ. Même pour des échantillons de taille modérée, cela peut être une opération lourde puisqu’il s’agit de calculer le déterminant d’une matrice carrée de dimension n.
Plusieurs solutions sont proposées dans la littérature dont la plus ancienne et encore la plus utilisée proposée par Cliff et Ord [1981]. Les estimateurs de β,ρet θ est obtenus à partir des conditions de premier ordre usuelles :
∂L/∂β = e𝘫Σ−θ 1X = 0 (1.50)
∂L/∂θ
∂L/∂ρ = −Tr[W (I − ρW )−1] + e𝘫Σθ−1Wy = 0 (1.51)
= − 1 Tr[Σ−1(∂Σθ/∂θ )] + 1 e𝘫Σ−1(∂Σθ/∂θ )Σ−1e = 0, (1.52)
i 2 θ i 2 θ i θ
pour i = 1, . . . , p, avec e = y − ρWy − Xβ et Tr l’opérateur trace. Les solutions de ces conditions découlent d’une optimisation numérique. L’inférence est basée sur une
matrice de variance asymtotique (asymptotic variance matrix) qui représente l’inverse de la matrice d’information. Dans le cas général, cette dernière prend la forme suivante pour [ρ, β, θ] :
X𝘫Σ−1W y X𝘫Σ−1X 0
ˆ θ θ
Tr(Wρ)2 + Tr[ΣθWρ𝘫 Σ−θ 1Wρ] + (Wˆy)𝘫Σθ−1(Wˆy) (X𝘫Σθ−1Wˆy)𝘫 ψ𝘫
ψ 0 1 Ψ
2
(1.53)
où Wρ = W (I −ρW )−1, W y = W (I −ρW )Xβ, et ψ un vecteur de traces de matrices de taille n × 1, ψi relatif au paramètre θi, et prend la forme :
ˆ
ψi = Tr[Σθ−1(∂Σθ/∂θi)Wρ] (1.54)
Il en est de même pour Ψ ; matrice de taille n × n contenant des traces de matrice et s’écrivant :
Ψi,j = Tr[Σθ−1(∂Σθ/∂θi)[Σ−θ 1(∂Σθ/∂θj )] (1.55)
La présence des vecteurs ψ implique que la matrice d’information n’est pas block- diagonale entre les paramètres du modèle [ρ, β] et les paramètres d’erreur θ. Il s’agit
d’une spécificité des modèles de décalage spatial entrainant un impact sur la nature des tests de spécification comme indiqué dans Anselin et Bera [1998].
L’échantillon peut contenir des sous-groupes : ce cas particulier est intéressant en pratique. En effet, la matrice variance-covariance Σ des erreurs est serait diagonale avec g = 1, . . . , p ses sous-diagonales, chacune correspondant à un sous-groupe :
σ2I1 0 · · · 0
1
Σ =
2
0 σ2I2 . . . 0
. . . . . 0
0 0 . . . σ2Ip
p
(1.56)
avec Ig,g = 1, . . . , p une matrice identité de ngéléments correspondant au nombre d’observations de chaque sous-groupe et Σgng = n.
En plaçant ces expressions ainsi formulées dans les conditions de premier ordre [1.50]-[1.52], des solutions analytiques assujetties aux valeurs de ρ sont obtenues pour
βMLet σg,ML :
2
σ
g,MV
βMV = (X𝘫Σ−1X)−1X𝘫Σ−1(y − ρWy) (1.57)
= (y − ρWy − Xβ)−1Hg(y − ρWy − Xβ)/ng (1.58)
pour g = 1, . . . , p et Hg une matrice identité avec tous les éléments de la diagonale nuls à l’exception de de ceux relatifs au sous-groupe g, faisant de [1.58] une somme moyenne des carrés des résidus du groupe respectif. L’estimateur de βMV est une somme pondérée, ayant pour poids −ρ de X par y et de X par Wy (X estimé par la méthode Moindres Carrés Généralisés Réalisables).
Contrairement à la situation d’erreurs i.i.d. (Anselin[1980]), la procédure précédente ne conduit pas à une simple vraissemblance concentrée ; les MCGR dépendent des termes σ2, qui eux dépendent du paramètre inconnu ρ . Les estimateurs peuvent être obtenus de façon itérative en substituant les valeurs de β et Σ d’une itération précédente dans [1.51] résolvant cette condition sur ρ .
Ainsi, de nouvelles valeurs de β et Σ sont obtenues à travers [1.57] et [1.58]. Par ailleurs, Σ−1 est une matrice diagonale avec des éléments 1/σ2 correspondant à chaque sous- groupe. L’inférence asymptotique peut être basée sur la matrice d’information [1.53], avec Σ−1 dans sa plus simple expression dérivée de la forme générale avec ∂Σ/∂σ2 = Hg.
g
g
g
Le modèle d’autocorrélation des erreurs :
L’estimation par maximum de vraissemblance des modèles d’autocorrélation spatiale des erreurs est un cas spécial des résultats de Magnus [1978]. Pour des termes d’erreurs
non-sphériques Σθ, dont le paramètre est θ, l’estimateur par ML de β est l’expression usuelle par les MCG :
βˆM V = (X𝘫Σ−1X)−1X𝘫Σ−1y (1.59)
θ θ
Cet estimateur [1.59] représente la solution des conditions de premier ordre de la vraissemblance :
L = − n ln 2π − 1 ln |Σ | − (y − Xβ)−1Σ−1(y − Xβ) (1.60)
2 2 θ θ
Quant aux estimateurs de θi, la solution s’obtient à partir des conditions de premier ordre suivantes :
Tr[(∂Σ−1/∂θ )Σ ] = e𝘫(∂Σ−1/∂θ )e (1.61)
θ i θ θ i
avec e = y − XβˆM L. L’estimateur convergent de β peut ainsi être obtenu par les MCGR sur les estimateurs convergents de θi. Similairement aux modèle de décalage spatial, l’inférence asymptotique est basée sur l’inverse de la matrice d’information block-diagonale entre β et les paramètres de variance-covariance des erreurs. La variance asymptotique de βˆM V prend la forme MCG :
AsyV (βˆM V ) = (X𝘫Σ−1X)−1 (1.62)
θ
Avec le block correspondant aux paramètres de variance-covariance de la forme 2Ψ−1, les éléments de Ψ étant ceux de [1.55].
Le rang du processus des erreurs considéré dans la section (1.2.1) va déboucher sur des formes spéciales de Σθ de manière à simplifier [1.59]-[1.61].
En guise d’illustration, le vecteur des paramètres relatif à un processus SAR (avec erreurs homoscédastiques) est : θ = [σ2, λ]. La matrice de varicance-covariance est [1.18]. Ainsi, l’estimateur par MCGF de ce modèle est simplifié en :
βˆM V = [X𝘫(I − λˆW )𝘫(I − λˆW )X]−1X𝘫(I − λˆW )𝘫(I − λˆW )y (1.63)
Cette simplification peut être considérée comme la regression de la variable avec filtrage spatial XL = X − λˆW X sur la variable, également avec filtrage spatial, yL = y − λˆW y . Il s’agit d’une forme de moindres carrés avec filtrage spatial. La dérivée partielle à utiliser dans la condition [1.61] est la suivante :
∂Σ−1/∂λ = −W − W 𝘫 + λW 𝘫W (1.64)
Contrairement à son homologue des séries temporelles, une estimation convergente de λ ne peut pas être obtenue à partir d’une régression auxiliaire simple. Ainsi, la condi- tion de premier ordre [1.61] doit être résolue par des moyens numériques. L’inférence
asymptotique est basée sur l’inverse de la matrice d’information comme pour le modèle de décalage spatial (Anselin[1988b]).
La régression avec filtrage spatial [1.63] aboutit à une expression simple de βˆM L. Ce résultat n’est pas valable pour d’autres processus spatiaux. Par exemple, dans le processus SMA de la section (1.2.1),la matrice de variance-covariance est donnée par [1.23] et :
Σ−1 = (1/σ2)[I + γ(W + W 𝘫) + γ2WW 𝘫]−1 (1.65)
nécessitant l’inverse d’une matrice de taille n × n pour effectuer des MCGF. De même, les dérivées partielles nécessaires pour la condition (73) impliquent inverses de cet ordre. Les modèles de représentation directe ou les modèles avec des erreurs compsées souffrent du même problème. En outre, les formes fonctionnelles et de distance métriques utilisées dans l’approche par la représentation directe peuvent causer des problèmes avec l’optimasation de la fonction log-vraisemblance (Anselin [2006]).
Les problèmes de Calcul :
L’estimation par le maximum de vraisemblance dans les modèles de régression spa- tiales implique l’application des techniques d’optimisation non-linéaire à la fonction de log-vraisemblance. Le problème principal de l’estimation par le maximum de vrai- semblance est la présence du log-jacobien dans la fonction de log-vraisemblance, de la forme ln|I − ρW |. De plus, les conditions de premier ordre et la matrice d’informa- tion contiennent des traces de produits de matrices dont W (I − ρW )−1. Le calcul de ces termes, même pour des échantillons de taille modérée, de manière habituelle est difficile.
La solution la plus ancienne et encore la plus utilisée est proposée par Ord [1975]. Elle consiste à exploiter la décomposition du jacobien en termes des n valeurs propres de la matrice de poids W :
soit :
n
I ρW = (1 ρωi) (1.66)
| − | −
i=1
n
| − | −
ln I ρW = ln(1 ρωi) (1.67)
i=1
L’apport de cette forme apparaît lors de la maximisation de la fonction de log- vraisemblance complète ou de la procédure de balayage de la fonction de log-vraisemblance
concentrée. Au lieu de calculer le déterminant de I − ρW à chaque étape, il n’est besoin
de le calculer qu’une seule fois pour toutes les valeurs propres de W et l’évaluation à chaque étape devient plus commode. Cette propriété permet alors d’écrire la fonction
de log-vraisemblance complète en une somme d’éléments correspondant aux observa- tions individuelles. Lorsque le calcul des valeurs propres pour une matrice de poids W asymétrique pose problème, Ord [1975] propose une propriété intéressante dans le cas où la matrice de poids W est le résultat de la standardisation d’une matrice de poids symétrique, notée W ∗.
Cet auteur démontre que les valeurs propres de W se trouvent en calculant les valeurs propres de la matrice symétrique,D1/2W ∗D1/2 où D est une matrice diagonale dans laquelle chaque élément de la diagonale principale est égal à l’inverse de la somme des éléments de la ligne correspondante de W .
Cette méthode est aujour- d’hui encore très largement utilisée. Pourtant, pour de très grands échantillons, le calcul des valeurs propres peut devenir numériquement instable. Une variante de cette pro- cédure est suggérée par Smirnov et Anselin [2001]. Elle consiste à évaluer directement les coefficients de la fonction caractéristique. Les simulations effectuées par les auteurs indiquent que cette méthode est la seule méthode « directe » capable de calculer les estimations du maximum de vraisemblance pour de très grands échantillons (plus d’un million d’observations).
D’autres techniques sont proposées pour de très grands échan- tillons. Certaines d’entre elles exploitent la structure particulière des matrices de poids qui contiennent en général beaucoup d’éléments nuls. Il s’agit des méthodes de facto- risation de matrices qui s’avèrent très puissantes pour évaluer rapidement le jacobien : la décomposition de Cholesky pour une matrice symétrique et la décomposition de LU sinon (Anselin[2006]).