Category Archives: 2012-H

2012-H, ACT6420, Time Series

Unit root, or not ? is it a big deal ?

10/09/2012 Arthur Charpentier 1 Comment

Consider a time series, generated using

set.seed(1)
E=rnorm(240)
X=rep(NA,240)
rho=0.8
X[1]=0
for(t in 2:240){X[t]=rho*X[t-1]+E[t]}

The idea is to assume that an autoregressive model can be considered, but we don’t know the value of the parameter. More precisely, we can’t choose if the parameter is either one (and the series is integrated), or some value strictly smaller than 1 (and the series is stationary). Based on past observations, the higher the autocorrelation, the lower the variance of the noise.

rhoest=0.9; H=260
u=241:(240+H)
P=X[240]*rhoest^(1:H)
s=sqrt(1/(sum((rhoest^(2*(1:300))))))*sd(X)

Now that we have a model, consider the following forecast, including a confidence interval,

 
plot(1:240,X,xlab="",xlim=c(0,240+H),
ylim=c(-9.25,9),ylab="",type="l")
V=cumsum(rhoest^(2*(1:H)))*s
polygon(c(u,rev(u)),c(P+1.96*sqrt(V),
rev(P-1.96*sqrt(V))),col="yellow",border=NA)
polygon(c(u,rev(u)),c(P+1.64*sqrt(V),
rev(P-1.64*sqrt(V))),col="orange",border=NA)
lines(u,P,col="red")

Here, forecasts can be derived, with any kind of possible autoregressive coefficient, from 0.7 to 1. I.e. we can chose to model the time series either with a stationary, or an integrated series,

As we can see above, assuming either that the series is stationary (parameter lower – strictly – than 1) or integrated (parameter equal to 1), the shape of the prediction can be quite different. So yes, assuming an integrated model is a big deal, since it has a strong impact on predictions.

2012-H, Time Series

Les tests de non-stationnarité (racine unité)

10/09/2012 Arthur Charpentier 3 Comments

Pour commencer la modélisation d’une série temporelle, la première étape est de savoir si on peut la considérer comme stationnaire. L’alternative étant qu’elle soit intégrée. Le test le plus classique est le test de Dickey-Fuller.De manière générale, considérons un modèle autorégressive à l’ordre 1,

On peut éventuellement réécrire ce processus sous la forme

où . On aura un processus non-stationnaire dès lors que ou encore . Si on n’a pas de tendance linéaire, tester ou peut se faire avec un test de Student. L’idée de Dickey-Fuller est de généraliser le test de Student, en rajoutant cette tendance linéaire (pour commencer).

Afin d’illustrer l’utilisation de ce test, commençons par un cas “simple”: on va intégrer un bruit blanc simulé.

>  set.seed(1)
>  E=rnorm(240)
>  X=cumsum(E)
>  plot(X,type="l")

On peut regarder s’il y a une constante non nulle (on parlera de modèle avec drift) voire une tendance linéaire (on parlera de trend),

La première étape pourrait être de tenter un modèle sans rien.

> df=ur.df(X,type="none",lags=1)
> summary(df)

############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression none

Call:
lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)

Residuals:
Min       1Q   Median       3Q      Max
-2.87492 -0.53977 -0.00688  0.64481  2.47556

Coefficients:
Estimate Std. Error t value Pr(>|t|)
z.lag.1    -0.005394   0.007361  -0.733    0.464
z.diff.lag -0.028972   0.065113  -0.445    0.657

Residual standard error: 0.9666 on 236 degrees of freedom
Multiple R-squared: 0.003292,	Adjusted R-squared: -0.005155
F-statistic: 0.3898 on 2 and 236 DF,  p-value: 0.6777

Value of test-statistic is: -0.7328

Critical values for test statistics:
1pct  5pct 10pct
tau1 -2.58 -1.95 -1.62

La région critique du test est ici (pour un seuil à 5%) l’ensemble des valeurs inférieures à -1.95. Or ici la statistique de test est ici -0.73, on est alors dans la région d’acceptation du test, et on va retenir l’hypothèse de racine unité, i.e. la série n’est pas stationnaire. Mais peut-être faudrait-il juste prendre en compte une constante ?

> df=ur.df(X,type="drift",lags=1)
> summary(df)

############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression drift

Call:
lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)

Residuals:
Min       1Q   Median       3Q      Max
-2.91930 -0.56731 -0.00548  0.62932  2.45178

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.29175    0.13153   2.218   0.0275 *
z.lag.1     -0.03559    0.01545  -2.304   0.0221 *
z.diff.lag  -0.01976    0.06471  -0.305   0.7603
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.9586 on 235 degrees of freedom
Multiple R-squared: 0.02313,	Adjusted R-squared: 0.01482
F-statistic: 2.782 on 2 and 235 DF,  p-value: 0.06393

Value of test-statistic is: -2.3039 2.7329

Critical values for test statistics:
1pct  5pct 10pct
tau2 -3.46 -2.88 -2.57
phi1  6.52  4.63  3.81

Deux statistiques de test sont calculées, ici: la première relative à la racine unité, la seconde relative à la constante. On observe ici que la statistique de test pour le test de racine unité (-2.3) est ici supérieure à toutes les valeurs critiques associées (données sur la première ligne). On va encore accepter l’hypothèse de racine unité. Mais le modèle était peut-être faux, et peut-être avait-on en fait en tendance linéaire ?

> df=ur.df(X,type="trend",lags=1)
> summary(df)

############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression trend

Call:
lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)

Residuals:
Min       1Q   Median       3Q      Max
-2.87727 -0.58802 -0.00175  0.60359  2.47789

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.3227245  0.1502083   2.149   0.0327 *
z.lag.1     -0.0329780  0.0166319  -1.983   0.0486 *
tt          -0.0004194  0.0009767  -0.429   0.6680
z.diff.lag  -0.0230547  0.0652767  -0.353   0.7243
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.9603 on 234 degrees of freedom
Multiple R-squared: 0.0239,	Adjusted R-squared: 0.01139
F-statistic:  1.91 on 3 and 234 DF,  p-value: 0.1287

Value of test-statistic is: -1.9828 1.8771 2.7371

Critical values for test statistics:
1pct  5pct 10pct
tau3 -3.99 -3.43 -3.13
phi2  6.22  4.75  4.07
phi3  8.43  6.49  5.47

On obtient cette fois trois statistique, la première relative au test de racine unité, et les deux suivantes aux tests sur la constante, et sur la tendance (la pente de l’ajustement linéaire). Là encore, la valeur de test (-1.98) excède les valeurs critiques: la p-value serait ici supérieure à 10%. On va là encore accepter l’hypothèse de racine unité.
Mais la modélisation autorégressive à l’ordre 1 était peut-être elle aussi fausse. Aussi, il existe un test de Dickey-Fuller augmenté. L’idée est de considérer, de manière beaucoup plus générale, un processus autorégressif à un ordre supérieur,

Comme auparavant, on peut réécrire ce modèle en fonction des variations (on parlera de modèle à correction d’erreur)

où on construit les nouveaux coefficients par une relation de récurrence, du genre . Bref, à nouveau, on est ramené à tester ou . Et là encore, on a le choix sur la tendance.
Si on commence par supposer qu’il n’y a pas de tendance,

> df=ur.df(X,type="none",lags=2)
> summary(df)

############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression none

Call:
lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)

Residuals:
Min       1Q   Median       3Q      Max
-2.86738 -0.53887 -0.02009  0.67058  2.45970

Coefficients:
Estimate Std. Error t value Pr(>|t|)
z.lag.1     -0.005168   0.007389  -0.699    0.485
z.diff.lag1 -0.029619   0.065289  -0.454    0.650
z.diff.lag2 -0.037856   0.065436  -0.579    0.563

Residual standard error: 0.9685 on 234 degrees of freedom
Multiple R-squared: 0.004704,	Adjusted R-squared: -0.008056
F-statistic: 0.3687 on 3 and 234 DF,  p-value: 0.7757

Value of test-statistic is: -0.6994

Critical values for test statistics:
1pct  5pct 10pct
tau1 -2.58 -1.95 -1.62

la valeur de la statistique de test excède là encore les valeurs critiques, i.e. la p-value dépasse 10%. Si on rajoute une constante,

> df=ur.df(X,type="drift",lags=2)
> summary(df)

############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression drift

Call:
lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)

Residuals:
Min       1Q   Median       3Q      Max
-2.91609 -0.56702  0.01025  0.62109  2.43970

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.31105    0.13332   2.333   0.0205 *
z.lag.1     -0.03744    0.01565  -2.392   0.0175 *
z.diff.lag1 -0.01917    0.06483  -0.296   0.7677
z.diff.lag2 -0.02794    0.06496  -0.430   0.6676
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.9594 on 233 degrees of freedom
Multiple R-squared: 0.02664,	Adjusted R-squared: 0.0141
F-statistic: 2.125 on 3 and 233 DF,  p-value: 0.09774

Value of test-statistic is: -2.3924 2.9709

Critical values for test statistics:
1pct  5pct 10pct
tau2 -3.46 -2.88 -2.57
phi1  6.52  4.63  3.81

on valide encore l’hypothèse de racine unité et de même avec une tendance linéaire,

> df=ur.df(X,type="trend",lags=2)
> summary(df)

############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression trend

Call:
lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)

Residuals:
Min       1Q   Median       3Q      Max
-2.85835 -0.58826 -0.00867  0.61407  2.47280

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.3530591  0.1522373   2.319   0.0213 *
z.lag.1     -0.0338831  0.0168523  -2.011   0.0455 *
tt          -0.0005674  0.0009879  -0.574   0.5663
z.diff.lag1 -0.0237595  0.0654158  -0.363   0.7168
z.diff.lag2 -0.0328215  0.0656102  -0.500   0.6174
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.9608 on 232 degrees of freedom
Multiple R-squared: 0.02802,	Adjusted R-squared: 0.01126
F-statistic: 1.672 on 4 and 232 DF,  p-value: 0.1573

Value of test-statistic is: -2.0106 2.0849 3.0185

Critical values for test statistics:
1pct  5pct 10pct
tau3 -3.99 -3.43 -3.13
phi2  6.22  4.75  4.07
phi3  8.43  6.49  5.47

On notera toutefois que ces trois modèles nous suggèrent de ne pas retenir autant de retard, qui ne semblent pas significatifs.
On peut bien entendu faire ces tests sur de vraies données par exemple le niveau du Nil

> library(datasets)
> NILE=Nile

ou encore les taux d’intérêt américains,

> base=read.table(
+ "http://freakonometrics.free.fr/basedata.txt",
+ header=TRUE)
> Y=base[,"R"]
> Y=Y[(base$yr>=1960)&(base$yr<=1996.25)]
> TAUX=ts(Y,frequency = 4, start = c(1960, 1))

Par exemple, sur cette dernière, on rejette l’hypothèse de stationnarité

>  df=ur.df(y=TAUX,lags=3,type="drift")
>  summary(df)

############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression drift

Call:
lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)

Residuals:
Min      1Q  Median      3Q     Max
-3.1982 -0.2947 -0.0629  0.3524  3.1899

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.40609    0.16706   2.431 0.016362 *
z.lag.1     -0.06339    0.02500  -2.535 0.012354 *
z.diff.lag1  0.32613    0.08145   4.004 0.000102 ***
z.diff.lag2 -0.31102    0.08027  -3.875 0.000165 ***
z.diff.lag3  0.28712    0.08103   3.543 0.000541 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.7805 on 137 degrees of freedom
Multiple R-squared: 0.2051,	Adjusted R-squared: 0.1819
F-statistic: 8.839 on 4 and 137 DF,  p-value: 2.227e-06

Value of test-statistic is: -2.5354 3.2457

Critical values for test statistics:
1pct  5pct 10pct
tau2 -3.46 -2.88 -2.57
phi1  6.52  4.63  3.81

Mais toutes sortes d’autres tests (plus robustes) peuvent être faits. Les plus connus sont le test de Philips-Perron et le test dit KPSS. Pour ce dernier, il faut spécifier s’il l’on suppose que la série est de moyenne constante, ou si une tendance doit être prise en compte. Si l’on suppose une constante non-nulle, mais pas de tendance, on utilise

> summary(ur.kpss(X,type="mu"))

####################### 
# KPSS Unit Root Test # 
####################### 

Test is of type: mu with 4 lags.

Value of test-statistic is: 0.972

Critical value for a significance level of:
10pct  5pct 2.5pct  1pct
critical values 0.347 0.463  0.574 0.739

alors que pour prendre en compte une tendance

> summary(ur.kpss(X,type="tau"))

####################### 
# KPSS Unit Root Test # 
####################### 

Test is of type: tau with 4 lags.

Value of test-statistic is: 0.5057

Critical value for a significance level of:
10pct  5pct 2.5pct  1pct
critical values 0.119 0.146  0.176 0.216

Derrière se cache un teste du multiplicateur de Lagrange. L’hypothèse nulle correspond à l’absence de racine unité: plus la statistique de test est grande, plus on s’éloigne de la stationnarité (hypothèse nulle). Avec ces deux tests, on rejette là encore l’hypothèse de stationnarité de notre marche aléatoire simulée.

Pour le test de Philips-Perron, on a un test de type Dickey-Fuller,

> PP.test(X)

Phillips-Perron Unit Root Test

data:  X
Dickey-Fuller = -2.0116, Trunc lag parameter = 4, p-value = 0.571

Là encore, la p-value nous recommande de valider l’hypothèse de non-stationnarité. A titre de comparaison, si on avait travaillé sur la série différenciée, on aurait accepté l’hypothèse de stationnarité

> PP.test(diff(X))

Phillips-Perron Unit Root Test

data:  diff(X)
Dickey-Fuller = -15.9522, Trunc lag parameter = 4, p-value = 0.01

2012-H, Time Series

Simulation de séries temporelles

06/09/2012 Arthur Charpentier 1 Comment

Un billet rapide pour reprendre le code tapé en cours, la semaine passée. Considérons un processus autorégressif d’ordre 1, où est un bruit blanc, stationnaire, i.e. appartient à l’intervalle . Le code pour simuler un tel processus est

n=1000
bruit=rnorm(n)
phi1= .85
X=rep(NA,n)
X[1]=0
for(t in 2:n){X[t]=phi1*X[t-1]+bruit[t]}
plot(acf(X),lwd=5,col='blue')
plot(pacf(X),lwd=5,col='blue')

ou avec un autocorrélation au premier ordre négative,

phi1= -0.7

On peut aussi regarder un processus autorégressif au second ordre,

sur la figure ci-dessous (avec en haut à gauche le triangle de stationnarité du couple de paramètres).

phi1=  0.3
phi2=  0.5
X=rep(NA,n)
X[1:2]=0
for(t in 3:n){
X[t]=phi1*X[t-1]+phi2*X[t-2]+bruit[t]}

Histoire de changer un peu, on peut regarder un processus moyenne mobile au premier ordre, où est un paramètre dans .

theta1=  .8
X=rep(NA,n)
X[1]=0
for(t in 2:n){
X[t]=bruit[t]+theta1*bruit[t-1]}

ou une moyenne mobile du second ordre,

theta1= -.6
theta2=  .5
X=rep(NA,n)
X[1:2]=0
for(t in 3:n){
X[t]=bruit[t]+theta1*bruit[t-1]+
theta2*bruit[t-2]}

2012-H, ACT6420, Time Series

Lissage exponentiel et séries temporelles

20/02/2012 Arthur Charpentier Leave a comment

L’idée du lissage exponentiel est de supposer que , où on suppose que le bruit est ici un bruit indépendant (imprévisible). On va alors supposer que

http://freakonometrics.hypotheses.org/files/2015/12/lex03.gif

où désigne un poids attribué à la nouvelle observation dans la fonction de lissage .

On parle de lissage pour moyenne mobile à poids exponentiel (ewma, exponentially weighted moving average, évoqué ici ou là sur ce blog). En effet

i.e., en itérant

(même s’il faut faire attention à cette décomposition qui, en pratique, s’arrête avec la première observation).

La méthode de Holt-Winters propose de généraliser ce lissage, en introduisant deux (voire trois) composantes.

sera interprété comme un niveau
correspondant à la variation du processus
qui correspond à un cycle (éventuel)

On suppose désormais, dans cette méthode (dite de Holt-Winters), que

D’un point de vue plus pratique, la programmation se fait comme suit,

> HW=HoltWinters(X,alpha=.2,beta=0)
> plot(HW)

Cette trois composantes seront utilisées pour faire une prévision de la série: la prédiction sera de la forme , à la date , pour un horizon (avec les notations du précédant billet).

> P=predict(HW,24,prediction.interval=TRUE)
> plot(HW,xlim=range(c(time(X),time(P))))
> polygon(c(time(P),rev(time(P))),c(P[,2],rev(P[,3])),
+ col="yellow",border=NA)
> lines(P[,1],col="red",lwd=3)

2012-H

Example of linear regression model

06/02/2012 Arthur Charpentier 1 Comment

(via http://www.afc-cca.com/ page 158)

2012-H, ACT6420, Regression

On linear models with no constant and R2

02/02/2012 Arthur Charpentier 3 Comments

In econometrics course we always say to our students that “if you fit a linear model with no constant, then you might have trouble. For instance, you might have a negative R-squared”. So I tried to find databases on the internet such that, when we compute a linear regression, we actually obtain a negative R squared. I have generated hundreds to random databases that should exhibit such a property, in R. With no success. Perhaps to be more specific, I should explain what might happen if we do not include a constant in a linear model. Consider the following dataset, where points are on a straight line, with a negative slope, far from the origin, symmetric with respect to the first diagonal.

> x=1:3
> y=3:1
> plot(x,y)

Points are on a straight line, so it is actually possible to get a perfect linear model. But only if we integrate a constant in our model. This is related to the fact that the correlation between our two variates is -1,

> cor(x,y)
[1] -1

The least-square program is here

i.e. the estimate of the slope is

Numerically, we obtain

> sum(x*y)/sum(x^2)
[1] 0.7142857

which is the actual slope on the illustration above. If we compute the sum of squares of errors (as a function of the slope), we have here

> ssr=function(b){sum((y-b*x)^2)}
> SSR=Vectorize(ssr)
> B=seq(-1,3,by=.1)
> plot(B,SSR(B),ylim=c(0,ssr(3)),cex=.6,type="b")

so the value we have computed is actually the minimum of the sum of squares of errors. But note that the sum of squares always exceeds the total sum of squares in red on the graph above

> ssr(b)
[1] 6.857143
> sum((y-mean(y))^2)
[1] 2

Recall that the total “coefficient of variation“, denoted , is defined as

i.e.

> 1-ssr(b)/sum((y-mean(y))^2)
[1] -2.428571

which is negative. It is also sometimes defined as “the square of the sample correlation coefficient between the outcomes and their predicted values“. Here it would be related to

> cor(b*x,y)
[1] -1

so we would have a unit . So obviously, using the in a model without a constant would give odd results. But the weird part is that if we run that regression with R, we get

> summary(lm(y~0+x))

Call:
lm(formula = y ~ 0 + x)

Residuals:
1       2       3
2.2857  0.5714 -1.1429

Coefficients:
Estimate Std. Error t value Pr(>|t|)
x   0.7143     0.4949   1.443    0.286

Residual standard error: 1.852 on 2 degrees of freedom
Multiple R-squared: 0.5102,	Adjusted R-squared: 0.2653
F-statistic: 2.083 on 1 and 2 DF,  p-value: 0.2857

Here, the estimation is correct. But the we obtain tells us that the model is not that bad… ~~So if anyone knows what R computes, I’d be glad to know~~. The value given by R (thanks Vincent for asking me to look carefully at the R source code) is obtained using Pythagoras’s theorem to compute the total sum of square,

> sum((b*x)^2)/(sum((b*x)^2)+sum((y-b*x)^2))
[1] 0.5102041

So be careful, the might look good, but meaningless !

2012-H, ACT6420, Regression, Statistics

Normalité des estimateurs dans une régression

18/01/2012 Arthur Charpentier Leave a comment

Cette semaine en cours, on a évoqué la normalité des estimateurs dans une régression linéaire. En supposant la normalité des résidus, on a vu en cours que était un estimateur Gaussien. En particulier, chacun des estimateurs est alors Gaussien, au sens où , ce qui peut se visualiser sur le graphique suivant (la constante est en abscisse, et la pente en ordonnée), avec un intervalle de confiance à 95%,

En fait, la variance étant inconnue (mais pouvant être estimée), si on remplace la variance des résidus par leur estimateur, la loi de l’estimateur est une loi de Student. De même pour l’estimateur de la pente de la droite de régression,

Le fait que le couple soit également Gaussien fait que l’on peut construire non plus un intervalle de confiance (on n’est plus en dimension 1) mais une ellipse de confiance. On retient une forme elliptique car c’est la plus petite région dans laquelle on se trouvera avec une probabilité de 95% (comme discuté dans un vieux billet).

Mais pour mieux comprendre cette notion d’ellipse de confiance, le plus simple est d’avoir recours à du rééchantillonnage dans la base de données. En effet, comme on ne dispose que d’un jeu de données, on a un estimateur. Et la discussion sur la loi de notre estimateur est purement théorique. Pour visualiser la loi de notre estimateur, il faudrait des centaines, voire des milliers de bases similaires. Que l’on n’a pas. La solution est alors de tirer au hasard des points de notre échantillon (avec remise), i.e. de faire du bootstrap,

set.seed(1)
COEF=matrix(NA,10000,2)
for(s in 1:nrow(COEF)){
I=sample(1:nrow(cars),nrow(cars),replace=TRUE)
COEF[s,]=lm(dist~speed,data=cars[I,])$coefficients
}

Si on regarde ce qui se passe, tirage après tirage, on génère un paquet d’échantillon, et pour chaque échantillon, on ajuste une droite de régression.

Les distributions – sur tous ces échantillons – de la constante et de la pente semblent effectivement Gaussiennes,

hist(COEF[,1],col="light blue",prob=TRUE)
u=seq(min(COEF[,1]),max(COEF[,1]),length=500)
v=dnorm(u,mean(COEF[,1]),sd(COEF[,1]))
lines(u,v,lwd=3,col="red")
hist(COEF[,2],col="light blue",prob=TRUE)
u=seq(min(COEF[,2]),max(COEF[,2]),length=500)
v=dnorm(u,mean(COEF[,2]),sd(COEF[,2]))
lines(u,v,lwd=3,col="red")

pour la distribution de l’estimateur de la pente, alors que pour l’estimateur de la constante, on a la distribution suivante,

Si on regarde maintenant la loi jointe,

on retrouve une forme elliptique pour le nuage des points (et une forte corrélation négative entre les deux estimateurs). Mais si on creuse un peu

library(ellipse)
reg=lm(dist~speed,data=cars)
e=ellipse(reg)
plot(e,type="l",lwd=2)
polygon(e,col="light blue")
points(COEF,cex=.5)

l’ellipse ne coïncide pas (tout à fait) avec celle obtenu théoriquement. Ce qui laisse à penser que l’hypothèse de normalité des résidus est peut-être à revoir..,

Freakonometrics