# Advanced Econometrics: Quantiles and Expectiles

Next Thursday, I will give the fourth lecture of the PhD course on advanced tools for econometrics, on quantile and expectile regressions. Slides are available online.

# Inference and autoregressive processes

Consider a (stationary) autoregressive process, say of order 2,

for some white noise  with variance . Here is a code to generate such a process,

```> phi1=.5
> phi2=-.4
> sigma=1.5
> set.seed(1)
> n=240
> WN=rnorm(n,sd=sigma)
> Z=rep(NA,n)
> Z[1:2]=rnorm(2,0,1)
> for(t in 3:n){Z[t]=phi1*Z[t-1]+phi2*Z[t-2]+WN[t]}```

Here, we have to estimate two sets of parameters: the autoregressive coefficients, and the variance of the innovation process . There are (at least) three techniques to estimate those parameters.

• using least square regression

A natural idea is to see here a regression model, and thus, if we consider a matrix formulation,

Here we can run (conditional) ordinary least squares estimation,

```> base=data.frame(Y=Z[3:n],X1=Z[2:(n-1)],X2=Z[1:(n-2)])
> regression=lm(Y~0+X1+X2,data=base)
> summary(regression)

Call:
lm(formula = Y ~ 0 + X1 + X2, data = base)

Residuals:
Min      1Q  Median      3Q     Max
-4.3491 -0.8890 -0.0762  0.9601  3.6105

Coefficients:
Estimate Std. Error t value Pr(>|t|)
X1  0.45107    0.05924   7.615 6.34e-13 ***
X2 -0.41454    0.05924  -6.998 2.67e-11 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.449 on 236 degrees of freedom
Multiple R-squared: 0.2561,	Adjusted R-squared: 0.2497
F-statistic: 40.61 on 2 and 236 DF,  p-value: 6.949e-16

> regression\$coefficients
X1         X2
0.4510703 -0.4145365
> summary(regression)\$sigma
[1] 1.449276```
• using Yule-Walker equations

As we’ve seen in class, we can easily get the following equations for the autocovariance functions,

which can also be written

So we just have to solve a simple linear system of equations. Note that if we divide by the variance, those equations can be written in terms of the autocorrelation functions

The code is the following

```> rho1=cor(Z[1:(n-1)],Z[2:n])
> rho2=cor(Z[1:(n-2)],Z[3:n])
> A=matrix(c(1,rho1,rho1,1),2,2)
> b=matrix(c(rho1,rho2),2,1)
> (PHI=solve(A,b))
[,1]
[1,]  0.4517579
[2,] -0.4155920```

Now, we need to extract the estimated innovation process, from this set of parameters (note that it could be possible to include the variance term in Yule-Walker equations, to get a three dimensional linear equation)

```> estWN=base\$Y-(PHI[1]*base\$X1+PHI[2]*base\$X2)
> sd(estWN)
[1] 1.445706```

This estimator is probably not the best one (we can take into account that we’ve lost two degrees of freedom), but as a starting point, let us consider this one.

• using (conditional) likelihood estimators

Finally, we can assume some distribution for the innovation process. Thestandard model is a Gaussian model, i.e.

In that case, the conditional log likelihood (conditional since we set the first two observations here) is

```> CondLogLik=function(A,TS){
+ phi1=A[1];  phi2=A[2]
+ sigma=A[3]	; L=0
+ for(t in 3:length(TS)){
+ L=L+dnorm(TS[t],mean=phi1*TS[t-1]+
+ phi2*TS[t-2],sd=sigma,log=TRUE)}
+ return(-L)}```

Now, we can run standard optimization procedures,

```> LogL=function(A) CondLogLik(A,TS=Z)
> optim(c(0,0,1),LogL)
\$par
[1]  0.4509685 -0.4144938  1.4430930

\$value
[1] 425.0164

\$counts
88       NA

\$convergence
[1] 0

\$message
NULL```

Here, our three estimators are rather close. Actually, if we generate 1,000 time series (of size 240), those are the Box-plots of our three estimators, for the first order autoregressive coefficient

for the second one,

and finally for the standard deviation of the innovation process

All those estimators behave nicely, and are rather close. Note that they all might be biased, but they are consistent (see Davidson and MacKinnon for instance, in their book, for more details).

# Régression, variables explicatives et géométrie

La régression (comme tout calcul d’espérance conditionnelle) est un problème de projections

Un théorème intéressant est le théorème dit de Frisch-Waugh, permettant de comprendre la différence fondamentale entre un modèle de régression multiple, et les modèles de régression simple. On reviendra sur ce point en cours en évoquant rapidement le paradoxe de Simpson. La formulation est la suivante: on veut régresser sur et , deux ensembles (a priori disjoints) de variables explicatives,

Les équations normales (associées au problème de minimisation de la somme des carrés des erreurs) sont

de telle sorte qu’à l’optimum on peut relier les estimateurs des deux jeux de paramètres par

La première partie correspond à la régression de sur , mais il reste un second terme dès lors que et ne sont pas orthogonales. Notons la matrice de projection (orthogonale) sur

Alors on peut écrire

En posant

et ,

on retrouve un modèle linéaire classique (à condition de travailler sur la projection des variables sur le sous-espace engendré par , i.e. en transformant les variables),

Pour aller plus loin sur la géométrie des moindres carrées, et sur le théorème de Frisch-Waugh, je peux renvoyer à des notes de cours, et à quelques transparents

Le graphique ci-dessous correspond au cas où les variables explicatives sont orthogonales, et dans ce cas, la régression multiple est équivalente à deux régressions simples

Le graphique ci-dessous au cas non orthogonal,

Le code qui permet de vérifier ces histoires de projections successives est relativement simple. Tout d’abord on importe les données, et on regarde le modèle global,

```> chicago=read.table(
+ "http://freakonometrics.free.fr/chicago.txt",
> Y=chicago\$Fire
> X1=chicago\$X_1
> X2=chicago\$X_2
> X3=chicago\$X_3
> base=data.frame(Y,X1,X2,X3)
> tail(base)
Y    X1 X2     X3
42  4.8 0.152 19 13.323
43 10.4 0.408 25 12.960
44 15.6 0.578 28 11.260
45  7.0 0.114  3 10.080
46  7.1 0.492 23 11.428
47  4.9 0.466 27 13.731
> regression=lm(Y~X1+X2+X3)
> summary(regression)

Call:
lm(formula = Y ~ X1 + X2 + X3)

Residuals:
Min     1Q Median     3Q    Max
-9.737 -4.565 -1.479  3.751 16.079

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 22.07525    6.19447   3.564 0.000910 ***
X1          -0.62764    5.28130  -0.119 0.905953
X2           0.22378    0.06161   3.632 0.000744 ***
X3          -1.55059    0.38195  -4.060 0.000204 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.527 on 43 degrees of freedom
Multiple R-squared: 0.4417,	Adjusted R-squared: 0.4027
F-statistic: 11.34 on 3 and 43 DF,  p-value: 1.314e-05

> n=length(Y)
> X=matrix(c(rep(1,n),X1,X2,X3),n,4)
> X[1:5,]
[,1]  [,2] [,3]   [,4]
[1,]    1 0.604   29 11.744
[2,]    1 0.765   44  9.323
[3,]    1 0.735   36  9.948
[4,]    1 0.669   37 10.656
[5,]    1 0.814   53  9.730```

Ensuite, on va projeter seulement sur les deux premières variables (et la constante)

```> FWX1=X[,1:3]
> regression12=lm(Y~X1+X2)
> solve(t(FWX1)%*%FWX1)%*%t(FWX1)%*%Y
[,1]
[1,]  0.08069764
[2,] 11.56913900
[3,]  0.15108490
> summary(regression12)\$coefficients
Estimate Std. Error    t value   Pr(>|t|)
(Intercept)  0.08069764 3.49182154 0.02311047 0.98166664
X1          11.56913900 5.05014782 2.29085156 0.02681947
X2           0.15108490 0.06854408 2.20420043 0.03278652
>
>
> FWX2=X[,4]
> H1=FWX1%*%solve(t(FWX1)%*%FWX1)%*%t(FWX1)
> M1=diag(rep(1,n))-H1
> FWX2s=M1%*%FWX2
> FWYs =M1%*%Y
> (beta2=solve(t(FWX2s)%*%FWX2s)%*%t(FWX2s)%*%FWYs)
[,1]
[1,] -1.550594
> summary(regression)\$coefficients[4]
[1] -1.550594
>
> (beta1=solve(t(FWX1)%*%FWX1)%*%t(FWX1)%*%Y-
+        solve(t(FWX1)%*%FWX1)%*%t(FWX1)%*%FWX2%*%beta2)
[,1]
[1,] 22.0752495
[2,] -0.6276442
[3,]  0.2237765
> summary(regression)\$coefficients[1:3]
[1] 22.0752495 -0.6276442  0.2237765```

On retrouve bien l’estimateur du dernier paramètre à l’aide du théorème de Frisch-Waugh, puis en utilisant les équations normales, on en déduit les trois premiers.