B.3 Variance de l’estimateur des moindres carrés ordinaires
B.3.1 Cas homoscédastique
Faire l’hypothèse d’homoscédasticité revient à supposer que \(\mathbb{E}[\epsilon^2 \mid X]\) est constante égale à un certain réel positif \(\sigma^2\).
Dans ce cas, le calcul en est grandement simplifié puisque la matrice de variance-covariance \(\mathbb{E}[XX']^{-1}\mathbb{E}[XX'\epsilon^2]\mathbb{E}[XX']^{-1}\) peut se simplifier en \(\sigma^2 \mathbb{E}[XX']^{-1}\). Un estimateur du premier terme est \(\frac{1}{n} \sum_{i=1}^n \hat{\epsilon_i}^2\). Un estimateur du second terme est \((\frac{1}{n}\mathbf{X}'\mathbf{X})^{-1}\) que l’on a de toute façon déjà calculé pour évaluer \(\hat{\beta}\). Ainsi, un estimateur de la matrice de variance-covariance de \(\sqrt{n}(\hat{\beta}-\beta)\) est donné par \((\frac{1}{n}\mathbf{X}'\mathbf{X})^{-1} \left(\frac{1}{n}\sum_{i=1}^n \hat{\epsilon_i}^2\right)\), et donc un estimateur de la matrice de variance-covariance de \(\hat{\beta}\) est donné par \(\frac{1}{n}(\frac{1}{n}\mathbf{X}'\mathbf{X})^{-1} \left(\frac{1}{n}\sum_{i=1}^n \hat{\epsilon_i}^2\right)\).
Le code suivant permet de vérifier que c’est bien ainsi que la fonction pré-implémentée lm estime la matrice de variance-covariance par défaut.
library(AER,
quietly=TRUE)
library(data.table,
quietly=TRUE)
#On charge dans un premier temps les données du CPS 1985
data("CPS1985")
CPS<-data.table(CPS1985)
#On estime la régression du salaire horaire sur l'éducation et l'âge
reg_wage_educ_age<-lm(wage~education + age,
data=CPS)
#On récupère les résidus estimés
CPS$estimated_residuals<-CPS$wage-
reg_wage_educ_age$coefficients["(Intercept)"]-
reg_wage_educ_age$coefficients["education"]*CPS$education-
reg_wage_educ_age$coefficients["age"]*CPS$age
#On estime sigma2 comme la moyenne empirique du carré des résidus
estimated_sigma2<-mean(CPS$estimated_residuals^2)
#On a aussi besoin de X'X
#On construit la matrice X
matrice_covariables<-as.matrix(cbind(rep(1,
nrow(CPS)),
CPS[,
c("education",
"age")]))
#On calcule X'X : c'est bien une matrice 3*3
XprimeX<-t(matrice_covariables)%*%matrice_covariables
XprimeX## V1 education age
## V1 534 6952 19669
## education 6952 94152 253613
## age 19669 253613 797769#L'estimateur de la matrice de variance-covariance : n(X'X)-1 multiplié par
# l'estimateur de sigma2 multiplié par 1/n
estimated_vcov<-(1)/(nrow(CPS)-nrow(XprimeX))*#On divise par n-d et pas par n pour
#avoir de meilleures propriétés à distance finie
solve((1/nrow(CPS))*XprimeX)*
estimated_sigma2
estimated_vcov## V1 education age
## V1 1.63677072 -0.0845952551 -0.0134615244
## education -0.08459526 0.0059360405 0.0001986127
## age -0.01346152 0.0001986127 0.0002952717#On vérifie l'égalité entre cet estimateur et la matrice de variance-covariance
# estimée par lm
all.equal(as.numeric(estimated_vcov),
as.numeric(vcov(reg_wage_educ_age)))## [1] TRUEB.3.2 Cas hétéroscédastique
Faire l’hypothèse d’hétéroscédasticité revient à autoriser la dispersion de \(\epsilon\), mesurée par l’espérance conditionnelle \(\mathbb{E}[\epsilon^2 \mid X]\), à varier avec les variables indépendantes.
Dans ce cas, il faut estimer la matrice \(\mathbb{E}[XX'\epsilon^2]\). La solution usuelle est de considérer comme estimateur de \(\mathbb{E}[XX']^{-1}\mathbb{E}[XX'\epsilon^2]\mathbb{E}[XX']^{-1}\) la variable aléatoire \((\frac{1}{n}\mathbf{X}'\mathbf{X})^{-1} (\frac{1}{n} \mathbf{X}' \hat{\Sigma} \mathbf{X}) (\frac{1}{n}\mathbf{X}'\mathbf{X})^{-1}\), où \(\hat{\Sigma}\) est la matrice carrée de taille \(n \times n\), dont les termes diagonaux sont égaux à \(\hat{\epsilon_i}^2\) et les termes non-diagonaux sont nuls. En effet, en décomposant le terme central :
\[\begin{align} \mathbf{X}' \hat{\Sigma} \mathbf{X} & = \left( \begin{array}{ccc} X_1^1 & \cdots & X_n^1 \\ \vdots & \ddots & \vdots \\ X_1^d & \cdots & X_n^d \\ \end{array}\right) \left( \begin{array}{ccc} \hat{\epsilon_1}^2 & \cdots & 0 \\ \vdots & \ddots \vdots \\ 0 & \cdots & \hat{\epsilon_n}^2 \\ \end{array}\right) \left(\begin{array}{ccc} X_1^1 & \cdots & X_i^d \\ \vdots & \ddots & \vdots \\ X_n^1 & \cdots & X_n^d \\ \end{array}\right) \nonumber \\ &= \left( \begin{array}{ccc} X_1^1 \hat{\epsilon_1}^2 & \cdots & X_1^n \hat{\epsilon_1}^2 \\ \vdots & \ddots & \vdots \\ X_1^d \hat{\epsilon_1}^2 & \cdots & X_n^d \hat{\epsilon_n}^2 \\ \end{array} \right) \left(\begin{array}{ccc} X_1^1 & \cdots & X_i^d \\ \vdots & \ddots & \vdots \\ X_n^1 & \cdots & X_n^d \\ \end{array}\right) \nonumber \\ &= \left( \begin{array}{ccc} \sum_{i=1}^n X_i^1 X_i^1 \hat{\epsilon_i}^2 & \cdots & \sum_{i=1}^n X_i^1 X_i^d \hat{\epsilon_i}^2 \\ \vdots & \ddots & \vdots \\ \sum_{i=1}^n X_i^d X_i^1 \hat{\epsilon_i}^2 & \cdots & \sum_{i=1}^n X_i^d X_i^d \hat{\epsilon_i}^2 \\ \end{array} \right) \end{align}\]
Le code suivant permet de vérifier que c’est bien cette approche qui sert à R pour calculer la matrice de variance-covariance robuste à l’hétéroscédasticité.
library(AER,
quietly=TRUE)
library(data.table,
quietly=TRUE)
#On charge dans un premier temps les données du CPS 1985
data("CPS1985")
CPS<-data.table(CPS1985)
#On estime la régression du salaire horaire sur l'éducation et l'âge
reg_wage_educ_age<-lm(wage~education + age,
data=CPS)
#On récupère les résidus estimés
CPS$estimated_residuals<-CPS$wage-
reg_wage_educ_age$coefficients["(Intercept)"]-
reg_wage_educ_age$coefficients["education"]*CPS$education-
reg_wage_educ_age$coefficients["age"]*CPS$age
#On construit la matrice X
matrice_covariables<-as.matrix(cbind(rep(1,
nrow(CPS)),
CPS[,
c("education",
"age")]))
#On construit la matrice diagonale Sigma
Sigma<-diag(x=CPS$estimated_residuals^2,
nrow=nrow(CPS),
ncol=nrow(CPS))
#On estime le terme central du sandwich : X'SigmaX
XprimeSigmaX<-1/nrow(CPS)*t(matrice_covariables)%*%Sigma%*%matrice_covariables
#On estime les termes des côtés
XprimeX<-t(matrice_covariables)%*%matrice_covariables
#On estime la matrice de variance-covariance
estimated_vcovHC<-1/(nrow(CPS))*
solve(1/nrow(CPS)*XprimeX)%*%
XprimeSigmaX%*%
solve(1/nrow(CPS)*XprimeX)
estimated_vcovHC## V1 education age
## V1 1.74565165 -0.0935543358 -0.0156346157
## education -0.09355434 0.0067706327 0.0003106483
## age -0.01563462 0.0003106483 0.0003213295#On vérifie que la matrice de variance-covariance robuste à l'hétéroscédaisticté
# estimée par R est égale à celle que l'on vient d'estimer
all.equal(as.numeric(estimated_vcovHC),
as.numeric(vcovHC(reg_wage_educ_age,
type="HC0")))## [1] TRUE#En fait le résultat par défaut prend en plus en compte une correction qui
# améliore la précision à distance finie
vcovHC(reg_wage_educ_age)## (Intercept) education age
## (Intercept) 1.78606223 -0.095847356 -0.0159540103
## education -0.09584736 0.006927378 0.0003197900
## age -0.01595401 0.000319790 0.0003272044B.3.3 En relâchant l’hypothèse d’indépendance
Grâce aux résultats de l’Annexe A.12, on sait que lorsque le nombre \(n_C\) de clusters devient très grand, la loi de la variable aléatoire \(\sqrt{n_C}(\hat{\beta}-\beta)\) s’identifie à la loi normale centrée de matrice de variance-covariance \(\mathbb{E}[\sum_{i \in c} X_iX_i']^{-1} \mathcal{V}(\tilde{\epsilon_c}) \mathbb{E}[\sum_{i \in c} X_iX_i']^{-1}\). Les termes latéraux de ce produit matriciel ne posent pas intrisèquement de difficulté : en effet, en reprenant les notations de la preuve en question, on voit que \(\frac{1}{n_C}\mathbf{X}'\mathbf{X}=\frac{1}{n_C}\sum_{c=1}^{n_C} \mathbf{X}_c'\mathbf{X}_c\) doit s’identifier quand \(n_C\) tend vers l’infini à cette matrice.
La difficulté est donc d’estimer le terme central \(\mathcal{V}(\tilde{\epsilon_c})\) avec : \[\tilde{\epsilon_c}:=\sum_{\begin{array}{c}i=1 \\ i \in c\end{array}}^n X_i \epsilon_i\] Pour ce faire, on va simplement remplacer dans la somme \(\epsilon_i\) qui est inconnu, par \(\hat{\epsilon_i}\). En définitive, le terme central de l’estimateur retenu de la matrice de variance-covariance de \(\sqrt{n_C}(\hat{\beta}-\beta)\) s’écrit : \[\frac{1}{n_C} \sum_{c=1}^{n_C} \left\{\left(\sum_{\begin{array}{c}i=1 \\ i \in c\end{array}}^n X_i \hat{\epsilon_i}\right)\left(\sum_{\begin{array}{c}i=1 \\ i \in c\end{array}}^n X_i \hat{\epsilon_i}\right)'\right\}\] Lorsque chaque observation correspond à un cluster, on retombe bien sur l’estimateur de la matrice de variance-covariance robuste à l’hétéroscédasticité examiné plus haut.
Le fragment de code suivant permet de vérifier sur l’exemple proposé en 2.9.5 que c’est bien de cette façon que R prend en compte les clusters dans le calcul des écarts-types.
library(AER,
quietly=TRUE)
library(data.table,
quietly=TRUE)
library(ggplot2,
quietly=TRUE)
#On charge dans un premier temps les données nécessaires
data("TeachingRatings")
TeachingRatings<-data.table(TeachingRatings)
#On réplique Hamermesh and Parker, 2005, Table 3
reg_eval<-lm(eval ~ beauty +
gender +
minority +
native +
tenure +
division +
credits,
weights = students,
data = TeachingRatings)
#On récupère les résidus estimés
TeachingRatings$residuals<-reg_eval$residuals
#On construit la matrice X
#Une difficulté ici : on a des variables qualitatives qu'il faut au préalable
# convertir en variables dichotomiques
#On récupère le nom des variables indépendantes incluses dans la régression
noms_covariables<-gsub(" ",
"",
strsplit(
strsplit(
as.character(reg_eval$call)[2],
"~")[[1]][2],
"\\+")[[1]])
#On récupère d'abord les colonnes des variables incluses dans la régression mais
# qui ne sont pas dans un format numérique
nonnumeric_variables<-names(sapply(TeachingRatings[,..noms_covariables],
is.numeric)[sapply(TeachingRatings
[,..noms_covariables],
is.numeric)==FALSE])
#Pour chaque colonne non-numérique il faut ensuite récupérer les niveaux de
# cette variable
nonnumeric_variables_levels<-lapply(TeachingRatings[,..nonnumeric_variables],
levels)
#Lorsqu'il y a strictement plus de 1 niveau il faut créer une variable
# dichotomique d'appartenance à chaque niveau
dichotomique<-function(variable,
niveau){
as.numeric(TeachingRatings[,..variable]==niveau)
}
#On boucle sur les niveaux à l'intérieur de chaque variable non-numérique
dichotomique_list<-function(variable){
if(variable %in% nonnumeric_variables){
if(length(nonnumeric_variables_levels[[variable]])>1){
dichot<-sapply(nonnumeric_variables_levels[[
variable]][
2:length(nonnumeric_variables_levels[[variable]])],
function(x){
dichotomique(variable,
x)
})
}
}
else{
TeachingRatings[,..variable]
}
}
matrice_covariables<-as.matrix(cbind(rep(1,
nrow(TeachingRatings)),
do.call(cbind,
sapply(noms_covariables,
dichotomique_list))))
colnames(matrice_covariables)[1]<-"(Intercept)"
#On est dans une régression pondérée, on va aussi avoir besoin du vecteur des
# poids
if(!is.null(reg_eval$weights)){
vecteur_poids<-reg_eval$weights
}
if(is.null(reg_eval$weights)){
vecteur_poids<-rep(1,
nrow(TeachingRatings))
}
#Une quantité utile : le nombre d'étudiants (le poids total des observations)
sum_weight<-as.numeric(sum(vecteur_poids))
#On définit le niveau de clustering
clustering<-"prof"
#Une quantité utile : n_C le nombre de clusters
nb_clusters<-nrow(unique(TeachingRatings[,..clustering]))
#On estime les termes des côtés
XprimeX<-
1/sum_weight*
t(sqrt(vecteur_poids)*matrice_covariables)%*%
(matrice_covariables*sqrt(vecteur_poids))
#La vraie difficulté : on veut estimer le terme central
#Il faut : sommer toutes les variables indépendantes*residus estimé à
# l'intérieur de chaque cluster = le groupe des observations correspondant à un
# enseignant
Xepsilon<-cbind(matrice_covariables*TeachingRatings$residuals,
TeachingRatings[,..clustering],
vecteur_poids)
colnames(Xepsilon)[ncol(Xepsilon)]<-"poids"
sum_Xepsilon_cluster<-
Xepsilon[,
lapply(X=.SD,
FUN=function(x){
sum(x*as.numeric(poids)/sum_weight)
}),
.SDcols=colnames(matrice_covariables),
by=clustering]
#Ensuite on multiplie la matrice obtenue par sa transposée
terme_central<-
1/nb_clusters*
t(as.matrix(
sum_Xepsilon_cluster[,
.SD,
.SDcols=names(sum_Xepsilon_cluster)!=clustering]))%*%
as.matrix(
sum_Xepsilon_cluster[,
.SD,
.SDcols=names(sum_Xepsilon_cluster)!=clustering])
#Enfin on multiplie par X'X^-1 des deux côtés
estimated_vcov_cluster<-
nb_clusters/(nb_clusters-1)*
1/nb_clusters*
solve(1/nb_clusters*XprimeX)%*%
terme_central%*%
solve(1/nb_clusters*XprimeX)
estimated_vcov_cluster## (Intercept) beauty.beauty gender.female minority.yes
## (Intercept) 0.0093537390 3.274401e-04 -0.0018700974 0.0008650045
## beauty.beauty 0.0003274401 3.396539e-03 -0.0010481125 -0.0008142540
## gender.female -0.0018700974 -1.048113e-03 0.0071363205 -0.0023394829
## minority.yes 0.0008650045 -8.142540e-04 -0.0023394829 0.0122835749
## native.no -0.0009147969 5.952765e-05 -0.0001297671 -0.0044308522
## tenure.yes -0.0074473070 2.740039e-04 -0.0002957024 -0.0016536652
## division.lower -0.0039995046 6.002078e-04 -0.0014972379 -0.0016677874
## credits.single 0.0004566126 1.006060e-03 0.0015614848 -0.0020578474
## native.no tenure.yes division.lower credits.single
## (Intercept) -9.147969e-04 -7.447307e-03 -0.0039995046 4.566126e-04
## beauty.beauty 5.952765e-05 2.740039e-04 0.0006002078 1.006060e-03
## gender.female -1.297671e-04 -2.957024e-04 -0.0014972379 1.561485e-03
## minority.yes -4.430852e-03 -1.653665e-03 -0.0016677874 -2.057847e-03
## native.no 1.768224e-02 -5.910749e-04 0.0024595592 -2.059399e-03
## tenure.yes -5.910749e-04 8.744043e-03 0.0025752894 6.397171e-05
## division.lower 2.459559e-03 2.575289e-03 0.0099946976 -5.861908e-03
## credits.single -2.059399e-03 6.397171e-05 -0.0058619076 2.698866e-02#On vérifie que la matrice de variance-covariance robuste à la dépendance
# estimée par R est égale à celle que l'on vient d'estimer
all.equal(as.numeric(estimated_vcov_cluster),
as.numeric(vcovCL(reg_eval,
type="HC0",
cluster=TeachingRatings[,..clustering])))## [1] TRUE#En fait l'estimateur par défaut utilise encore une correction supplémentaire
# qui prend en compte la dimension des variables indépendendantes (qui réduit
# les degrés de liberté)
all.equal((nrow(TeachingRatings)-1)/
(nrow(TeachingRatings)-ncol(matrice_covariables))*
as.numeric(estimated_vcov_cluster),
as.numeric(vcovCL(reg_eval,
cluster=TeachingRatings[,..clustering])))## [1] TRUE