# Probabilistic Foundations of Econometrics, part 4

This post is the fourth one of our series on the history and foundations of econometric and machine learning models. Part 3 is online here.

## Goodness of Fit, and Model

In the Gaussian linear model, the determination coefficient – noted $R^2$ – is often used as a measure of fit quality. It is based on the variance decomposition formula $$\underbrace{\frac{1}{n}\sum_{i=1}^n (y_i-\bar{y})^2}_{\text{total variance}}=\underbrace{\frac{1}{n}\sum_{i=1}^n (y_i-\widehat{y}_i)^2}_{\text{residual variance}}+\underbrace{\frac{1}{n}\sum_{i=1}^n (\widehat{y}_i-\bar{y})^2}_{\text{explained variance}}$$ The $R^2$ is defined as the ratio of explained variance and total variance, another interpretation of the coefficient that we had introduced from the geometry of the least squares $$R^2= \frac{\sum_{i=1}^n (y_i-\bar{y})^2-\sum_{i=1}^n (y_i-\widehat{y}_i)^2}{\sum_{i=1}^n (y_i-\bar{y})^2}$$The sums of the error squares in this writing can be rewritten as a log-likelihood. However, it should be remembered that, up to one additive constant (obtained with a saturated model) in generalized linear models, deviance is defined by ${Deviance}(\widehat{\beta}) = -2\log[\mathcal{L}]$ which can also be noted $Deviance(\widehat{\mathbf{y}})$. A null deviance can be defined as the one obtained without using the explanatory variables $\mathbf{x}$, so that $\widehat{y}_i=\overline{y}$. It is then possible to define, in a more general context (with a non-Gaussian distribution for $y$)$$R^2=\frac{{Deviance}(\overline{y})-{Deviance}(\widehat{\mathbf{y}})}{{Deviance}(\overline{y})}=1-\frac{{Deviance}(\widehat{\mathbf{y}})}{{Deviance}(\overline{y})}$$However, this measure cannot be used to choose a model, if one wishes to have a relatively simple model in the end, because it increases artificially with the addition of explanatory variables without significant effect. We will then tend to prefer the adjusted $R^2$,$$\bar R^2 = {1-(1-R^{2})\cdot{n-1 \over n-p}} = R^{2}-\underbrace{(1-R^{2})\cdot{p-1 \over n-p}}_{\text{penalty}}$$where $p$ is the number of parameters of the model. Measuring the quality of fit will penalize overly complex models.

This idea will be found in the Akaike criterion, where $AIC=Deviance+2\cdot p$ or in the Schwarz criterion, $BIC=Deviance+log(n)\cdot p$. In large dimensions (typically $p>\sqrt{n}$), we will tend to use a corrected $AIC$, defined by $AIC_c=Deviance+2⋅p⋅n/(n-p-1)$.

These criterias are used in so-called “stepwise” methods, introducing the set methods. In the “forward” method, we start by regressing to the constant, then we add one variable at a time, retaining the one that lowers the $AIC$ criterion the most, until adding a variable increases the AIC criterion of the model. In the “backward” method, we start by regressing on all variables, then we remove one variable at a time, removing the one that lowers the $AIC$ criterion the most, until removing a variable increases the $AIC$ criterion from the model.

Another justification for this notion of penalty (we will come back to this idea in machine learning) can be the following. Let us consider an estimator in the class of linear predictors,$$\mathcal{M}=\big\lbrace m:~m(\mathbf{x})=s_h(\mathbf{x})^T\mathbf{y} \text{ where }S=(s(\mathbf{x}_1),\cdots,s(\mathbf{x}_n))^T\text{ is some smoothing matrix}\big\rbrace$$ and assume that $y=m_0 (x)+\varepsilon$, with $\mathbb{E}[\varepsilon]=0$ and $Var[\varepsilon]=\sigma^2\mathbb{I}$, so that $m_0 (x)=\mathbb{E}[Y|X=x]$. From a theoretical point of view, the quadratic risk, associated with an estimated model $\widehat{m}$, $\mathbb{E}\big[(Y-\widehat{m}(\mathbf{X}))^2\big]$, is written$$\mathcal{R}(\widehat{m})=\underbrace{\mathbb{E}\big[(Y-m_0(\mathbf{X}))^2\big]}_{\text{error}}+\underbrace{\mathbb{E}\big[(m_0(\mathbf {X})-\mathbb{E}[\widehat{m}(\mathbf{X})])^2\big]}_{\text{bias}^2}+\underbrace{\mathbb{E}\big[(\mathbb{E}[\widehat{m}(\mathbf{X})]-\widehat{m}(\mathbf{X}))^2\big]}_{\text{variance}}$$if $m_0$ is the true model. The first term is sometimes called “Bayes error”, and does not depend on the estimator selected, $\widehat{m}$.

The empirical quadratic risk, associated with a model $m$, is here:$$\widehat{\mathcal{R}}_n(m)=\frac{1}{n}\sum_{i=1}^n (y_i-m(\mathbf{x}_i))^2$$ (by convention). We recognize here the mean square error, “mse”, which will more generally give the “risk” of the model $m$ when using another loss function (as we will discuss later on). It should be noted that:$$\displaystyle{\mathbb{E}[\widehat{\mathcal{R}}_n(m)]=\frac{1}{n}\|m_0(\mathbf{x})-m(\mathbf{x})\|^2+\frac{1}{n}\mathbb{E}\big(\|{Y}-m_0(\mathbf{X})\|^2\big)}$$We can show that:$$n\mathbb{E}\big[\widehat{\mathcal{R}}_n(\widehat{m})\big]=\mathbb{E}\big(\|Y-\widehat{m}(\mathbf{x})\|^2\big)=\|(\mathbb{I}-\mathbf{S})m_0\|^2+\sigma^2\|\mathbb{I}-\mathbf{S}\|^2$$so that the (real) risk of $\widehat{m}$ is:$${\mathcal{R}}_n(\widehat{m})=\mathbb{E}\big[\widehat{\mathcal{R}}_n(\widehat{m})\big]+2\frac{\sigma^2}{n}\text{trace}(\boldsymbol{S})$$So, if $\text{trace}(\boldsymbol{S})\geq0$ (which is not a too strong assumption), the empirical risk underestimates the true risk of the estimator. Actually, we recognize here the number of degrees of freedom of the model, the right-hand term corresponding to Mallow’s $C_p$, introduced in Mallows (1973) using not deviance but $R^2$.

## Statistical Tests

The most traditional test in econometrics is probably the significance test, corresponding to the nullity of a coefficient in a linear regression model. Formally, it is the test of $H_0:\beta_k=0$ against $H_1:\beta_k\neq 0$. The so-called Student test, based on the statistics $t_k=\widehat{\beta}_k/se_{\widehat{β}_k}$, allows to decide between the two alternatives, using the test $p$-value, defined by $\mathbb{P}[|T|>|t_k|]$ avec $T\overset{\mathcal{L}}{\sim} Std_\nu$, where $\nu$ is the number of degrees of freedom of the model ($\nu=p+1$ for the standard linear model). In large dimension, however, this statistic is of very limited interest, given a significant FDR (“False Discovery Ratio”). Classically, with a level of significance $\alpha=0.05$, 5% of the variables are falsely significant. Suppose that we have $p$=100 explanatory variables, but that 5 (only) are really significant. We can hope that these 5 variables will pass the Student test, but we can also expect that 5 additional variables (false positive test) will emerge. We will then have 10 variables perceived as significant, while only half are significant, i.e. an FDR ratio of 50%. In order to avoid this recurrent pitfall in multiple tests, it is natural to use the procedure of Benjamini & Hochberg (1995).

## From a correlation to some causal effect

Econometric models are used to implement public policy evaluations. It is therefore essential to fully understand the underlying mechanisms in order to know which variables actually make it possible to act on a variable of interest. But then we move on to another important dimension of econometrics. Jerry Neyman was responsible for the first work on the identification of causal mechanisms, and then Rubin (1974) formalized the test, called the “Rubin causal model” in Holland (1986). The first approaches to the notion of causality in econometrics were based on the use of instrumental variables, models with discontinuity of regression, analysis of differences in differences, and natural or unnatural experiments. Causality is usually inferred by comparing the effect of a policy – or more generally of a treatment – with its counterfactual, ideally given by a random control group. The causal effect of the treatment is then defined as $\Delta=y_1-y_0$, i.e. the difference between what the situation would be with treatment (noted t=1) and without treatment (noted t=0). The concern is that only $y=t\cdot y_1+(1-t)\cdot y_0$ and $t$ are observed. In other words, the causal effect of variable $t$  on $t$  is not observed (since only one of the two potential variables – $y_0$ or $y_1$  is observed for each individual), but it is also individual, and therefore a function of x-covariates. Generally, by making assumptions about the distribution of the triplet $(Y_0,Y_1,T)$, some parameters of the causal effect distribution become identifiable, based on the density of the observable variables $(Y,T)$. Classically, we will be interested in the moments of this distribution, in particular the average effect of treatment in the population, $\mathbb{E}[\Delta]$, or even just the average effect of treatment in the case of treatment $\mathbb{E}[\Delta|T=1]$. If the result $(Y_0,Y_1)$ is independent of the processing access variable $T$, it can be shown that $\mathbb{E}[\Delta]=\mathbb{E}[Y|T=1]- \mathbb{E} [Y|T=0]$. But if this independence hypothesis is not verified, there is a selection bias, often associated with $\mathbb{E}[Y_0|T=1]- \mathbb{E} [Y_0|T=0]$. Rosenbaum & Rubin (1983) propose to use a propensity to be treated score, $p(x)=\mathbb{P}[T=1|X=x]$, noting that if variable $Y_0$\ is independent of access to treatment $T$ conditionally to the explanatory variables $X$, then it is independent of $T$  conditionally to the score $p(X)$ : it is sufficient to match them using their propensity score. Heckman et al (2003) thus proposes a kernel estimator on the propensity score, which simply provides an estimator of the effect of the treatment, provided that it is treated.

To be continued next time, we’ll introduce “machine learning techniques” (references mentioned above are online here)

# Régression sur une variable qualitative et ANOVA

Ce matin, pour le cours STT5100, on évoquait la régression sur une variable catégorielle. En particulier, on avait commencé par regarder ce que donnerait la régression sans la constante, et son interprétation. On s’était appuyé sur la base des poids et des tailles des élèves, et de la variable de genre.

Davis=read.table( "http://socserv.socsci.mcmaster.ca/jfox/Books/Applied-Regression-2E/datasets/Davis.txt") Davis[12,c(2,3)]=Davis[12,c(3,2)] Davis=data.frame(Y=Davis$weight * 2.204622, X1=Davis$sex)

On voulait estimer le modèle $$y_i =\beta_F\boldsymbol{1}_F(x_i)+\beta_H\boldsymbol{1}_H(x_i)+\varepsilon_i$$On avait vu que l’on pouvait passer par l’écriture matricielle

 X=cbind(Davis$X1=='F',Davis$X1=='M') Y=Davis$Y car la matrice $\mathbf{X}^T\mathbf{X}$ est inversible (une fois que l’on enlève la constante)  solve(t(X)%*%X) [,1] [,2] [1,] 0.008928571 0.00000000 [2,] 0.000000000 0.01136364 et donc l’estimateur par moindres carrés est (classiquement)$$\widehat{\mathbf{\beta}} = (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{y}$$  solve(t(X)%*%X) %*% (t(X)%*%Y) [,1] [1,] 125.4272 [2,] 167.3258 ce qui correspond effectivement à la sortie de R,  reg=lm(Y~0+X1,data=Davis) summary(reg) Coefficients: Estimate Std. Error t value Pr(&gt;|t|) X1F 125.427 1.960 64.00 &lt;2e-16 *** X1M 167.326 2.211 75.68 &lt;2e-16 *** Considérons maintenant les deux sous-populations, avec le poids des femmes, et le poids des hommes x=Y[X[,1]==1] y=Y[X[,2]==1] nx=length(x) ny=length(y) On avait vu en cours que les $\widehat{\mathbf{\beta}}$ avaient une interprétation très simple, puisque$$\widehat{{\beta}}_M = \frac{1}{n_M}\sum_{i:x_i=M} y_i$$autrement dit $\widehat{{\beta}}_M$ est le poids moyen des hommes. Et en effet  mean(y)  167.3258 C’est finalement très naturel, ou intuitif. On peut maintenant s’interroger sur l’écart-type de l’estimateur de $\widehat{{\beta}}_M$. Intuitivement, on aurait envie d’avoir la variance de l’estimateur de la moyenne, soit ici  sqrt(var(y)/ny)  2.794391 sqrt(1/(ny-1)*sum( (y-mean(y))^2 )/ny)  2.794391 car pour rappel$$\text{Var}[\overline{y}]=\frac{\text{Var}(y)}{n}$$Comme on l’a vue dans le modèle de régression multiple, la variance de l’estimateur de $\mathbf{\beta}$ est proportionnel à $\sigma^2$, la variance globale des résidus (c’est l’hypothèse d’homoscédasticité ! les deux groupes doivent avoir la même variance). On va donc calculer l’estimateur naturel de $\sigma^2$  s2=1/(nx+ny-2)*(sum( (x-mean(x))^2 )+sum( (y-mean(y))^2)) sqrt(s2/ny)  2.210863 et en effet, on retombe sur la valeur donnée dans le tableau de régresion  sqrt(s2/nx)  1.959721 (pareil pour l’autre coefficient). On avait ensuite regardé la régression telle qu’elle faite classiquement, sous R : on garde la constante, et on enlève une des variables indicatrices (qui devient alors la “modalité de référence”).  X=cbind(1,Davis$X1=='M')

Là encore, le modèle devient identifiable, et on obtient ici

 solve(t(X)%*%X) %*% (t(X)%*%Y) [,1] [1,] 125.42724 [2,] 41.89855

On avait noté qu’il y avait un interprétation de cette seconde valeur, comme un différentiel par rapport à la modalité de référence

mean(y)-mean(x)  41.89855

La sortie de régression devient ici

 reg2=lm(Y~X1,data=Davis) summary(reg2)   Coefficients: Estimate Std. Error t value Pr(&gt;|t|) (Intercept) 125.427 1.960 64.00 &lt;2e-16 *** X1M 41.899 2.954 14.18 &lt;2e-16 ***

Et comme je l’avais dit, le test de Student correspond ici à un test d’égalité entre la taille moyenne des hommes et celle des femmes. Et en effet, si on fait le test, on voit que la différence est significative, comme attendu (pour la même raison qu’au dessus, on suppose la même variance dans les deux groupes)

 t.test(Y[X[,1]==1],Y[X[,2]==1],var.equal=TRUE)   Two Sample t-test   data: Y[X[, 1] == 1] and Y[X[, 2] == 1] t = -6.4475, df = 286, p-value = 4.826e-10 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -30.62603 -16.30035 sample estimates: mean of x mean of y 143.8626 167.3258

Je suis par contre un peu surpris que les p-values soient différente. Mon interprétation est que les p-values sont (de toutes façons) très faibles, et donc ça a peu d’importance. En fait, si on rend les deux variables indépendantes (par exemple en mélangeant la variable $\mathbf{y}$), ça marche ! Posons

 Davis$Y=sample(Davis$Y)

ce qui revient à permuter toutes les observations de la variable dépendante (mais pas les autres !). La régression donne ici

 reg2=lm(Y~X1,data=Davis) summary(reg2)   Call: lm(formula = Y ~ X1, data = Davis)   Residuals: Min 1Q Median 3Q Max -57.458 -22.184 -5.512 17.809 118.912   Coefficients: Estimate Std. Error t value Pr(&gt;|t|) (Intercept) 143.4382 2.7820 51.56 &lt;2e-16 *** X1M 0.9645 4.1940 0.23 0.818 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1   Residual standard error: 29.44 on 198 degrees of freedom Multiple R-squared: 0.000267, Adjusted R-squared: -0.004782 F-statistic: 0.05289 on 1 and 198 DF, p-value: 0.8183

autrement dit, le genre n’est plus significatif, avec une p-value de 81.8%. Ce qui est bien au dessus de 5%. Si on fait maintenant le test de comparaison de moyenne, sur les deux sous-groupes, on obtient

 Y=Davis$Y t.test(Y[X[,1]==1],Y[X[,2]==1],var.equal=TRUE) Two Sample t-test data: Y[X[, 1] == 1] and Y[X[, 2] == 1] t = -0.22998, df = 198, p-value = 0.8183 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -9.235209 7.306165 sample estimates: mean of x mean of y 143.4382 144.4027 et le test a ici également une p-value de 81.8%. Les deux tests sont donc rigoureusement équivalents. # Copulas and Financial Time Series I was recently asked to write a survey on copulas for financial time series. The paper is, so far, unfortunately, in French, and is available on https://hal.archives-ouvertes.fr/. There is a description of various models, including some graphs and statistical outputs, obtained from read data. To illustrate, I’ve been using weekly log-returns of (crude) oil prices, Brent, Dubaï and Maya. The dataset is available from an excel file, oil.xls (I thought it was possible to load it direclty from the internet, but it did not work… so I suggest to download the file first, and then load it) > library(xlsx) > temp <- tempfile() > download.file( + "http://freakonometrics.free.fr/oil.xls",temp) trying URL 'http://freakonometrics.free.fr/oil.xls' Content type 'application/vnd.ms-excel' length 99328 bytes (97 KB) downloaded 97 KB > oil=read.xlsx(temp,sheetName="DATA",dec=",") Error in .jcall("RJavaTools", "Ljava/lang/Object;", "invokeMethod", cl, : java.io.IOException: block[ 0 ] already removed - does your POIFS have circular or duplicate block references? > oil=read.xlsx("D:\\home\\acharpen\\mes documents\\oil.xls",sheetName="DATA") Then we can plot those three time series > head(oil) Date WTI brent Dubai Maya 1 1997-01-10 2.73672 2.25465 3.3673 1.5400 2 1997-01-17 -3.40326 -6.01433 -3.8249 -4.1076 3 1997-01-24 -4.09531 -1.43076 -6.6375 -4.6166 4 1997-01-31 -0.65789 0.34873 0.7326 -1.5122 5 1997-02-07 -3.14293 -1.97765 -0.7326 -1.8798 6 1997-02-14 -5.60321 -7.84534 -7.6372 -11.0549 > Time=as.Date(oil$Date,"%Y-%m-%d")
> plot(Time,oil[,3],type="l",ylab="Brent, weekly log returns",ylim=range(oil[,3:5])) The idea is to use some multivariate ARMA-GARCH processes here. The heuristics here is that the first part is used to model the dynamics of the average value of the time series, and the second part is used to model the dynamics of the variance of the time series. Two kinds of models are considered in the paper

• a mutivariate GARCH process (or a model on the dynamics of the variance matrix) on the residuals from the ARMA models
• a multivariate model (based on copulas) on the residuals of the ARMA-GARCH process

# Conditional Distributions from some Elliptical Vectors

This winter, in my ACT8595 course, I asked my students (that was some homework) to prove that it was possible to derive the conditional distribution when we have a Student-t random vector (and to get the analytical expression of the later). But before, let us recall a standard result about the Gaussian vector. If  is a Gaussian random vector, i.e.

then  has a Gaussian distribution. More precisely, it is a  distribution, with

and  is the Schur complement of the block  of the matrix ,

Observe that  is also related to well known quantity: in the bivariate case, where  and  are univariate Gaussian variables,

which is the slope in the linear regression of  on .

In the case of the Student-t distribution, the conditional distrubution will not be a Student-t distribution anymore, but it will still be an elliptical distribution, and some interpretations of various quantities can actually be obtained.

The density of the multivariate centred Student-t distribution, with unit variance, and parameters  and  is

If we consider the following blocks,

then we can get that marginal distributions have a centred Student-t distribution, with unit variance, and parameters  and ,

Then, to derive the conditional density, we can use Bayes formula,

One can write (as in Section 9.1 in Tong, 1990, The Multivariate Normal Distribution, but other expressions can be found in Section 2.5 in Fang, Ng and Kotz, 1989, Symmetric multivariate and related distributions, or in Section 1.11 in Kotz and Nadarajah, 2004, Multivariate t distributions and their applications) this conditional density as

with

and

This conditional distribution is elliptical, but it is not a Student-t distribution, except in the case where , or when the correlation matrix  is the identity.

[June 2016] actually, as shown in Ding (2016), this is a Student-t distribution. “Kotz & Nadarajah (2004) and Nadarajah & Kotz (2005) failed to recognize that the conditional distribution of the MVT distribution is also a MVT distribution due to the complexity of the conditional density function […] Conditional distributions of elliptically contoured distributions are also elliptically contoured distributions. But this does not immediately guarantee that conditional distributions of the MVT distributions are also MVT distributions without some further algebra.

Now, if we look at the components of this density, we can observe that we have

which was mentioned previously, in the Gaussian case: the term on the right is the conditional mean,

and the bloc that appears at several places is the conditional variance,

Now, if we want to visualize that conditional density, let us plot it. The code below is based on Bayes formula

> library(mnormt)
> r=.6
> R=matrix(c(1,r,r,1),2,2)
> nu=4
> f2=function(x2) dt(x2,df=nu)
> f =function(x) dmt(x,S=R,df=nu)
> f1.2=function(x1,x2) f(c(x1,x2))/f2(x2)

In order to compare that conditional density with a Student-t one, let us define the density of a non-centred Student-t random variable,

> dstd=function(x,mu,s,nu) gamma((nu+1)/2)/
+ (gamma(nu/2)*s*sqrt(pi*nu))*
+ (1+1/nu*(x-mu)^2/(s^2))^(-(nu+1)/2)

Here is the function we can use to plot those two densities,

> graphdensity=function(x2=-1.5){
+ vectx1=seq(-3,3,length=251)
+ y=Vectorize(function(x) f1.2(x,x2))(vectx1)
+ plot(vectx1,y,type="l",col="red",ylim=c(0,.5),
+ xlab="",ylab="")
+ abline(v=r*x2,lty=2)
+ lines(vectx1,dstd(vectx1,x2*r,sqrt(1-r^2),nu),col="blue",lty=2)}
> graphdensity(-1.5) In the case where , the two lines are rather close (the difference migth come from computational issues)

> graphdensity(-1) and just to conclude, a last one

> graphdensity(0) # Bivariate Densities with N(0,1) Margins

This Monday, in the ACT8595 course, we came back on elliptical distributions and conditional independence (here is an old post on de Finetti’s theorem, and the extension to Hewitt-Savage’s). I have shown simulations, to illustrate those two concepts of dependent variables, but I wanted to spend some time to visualize densities. More specifically what could be the joint density is we assume that margins are  distributions.

• The Bivariate Gaussian distribution

Here, we consider a Gaussian random vector, with margins , and with correlation . This is the standard graph, with elliptical isodensity curves

r=.5
library(mnormt)
S=matrix(c(1,r,r,1),2,2)
f=function(x,y) dmnorm(cbind(x,y),varcov=S)
vx=seq(-3,3,length=201)
vy=seq(-3,3,length=201)
z=outer(vx,vy,f)
set.seed(1)
X=rmnorm(1500,varcov=S)
xhist <- hist(X[,1], plot=FALSE)
yhist <- hist(X[,2], plot=FALSE)
top <- max(c(xhist$density, yhist$density,dnorm(0)))
nf <- layout(matrix(c(2,0,1,3),2,2,byrow=TRUE), c(3,1), c(1,3), TRUE)
par(mar=c(3,3,1,1))
image(vx,vy,z,col=rev(heat.colors(101)))
points(X,cex=.2)
par(mar=c(0,3,1,1))
barplot(xhist$density, axes=FALSE, ylim=c(0, top), space=0,col="light green") lines((density(X[,1])$x-xhist$breaks)/diff(xhist$breaks),
dnorm(density(X[,1])$x),col="red") par(mar=c(3,0,1,1)) barplot(yhist$density, axes=FALSE, xlim=c(0, top), space=0,
horiz=TRUE,col="light green")
lines(dnorm(density(X[,2])$x),(density(X[,2])$x-yhist$breaks)/ diff(yhist$breaks),col="red") That was the simple part.

• The Bivariate Student-t distribution

Consider now another elliptical distribution. But we want here to normalize the margins. Thus, instead of a pair , we would like to consider the pair , so that the marginal distributions are . The new density is obtained simply since the transformation is a one-to-one increasing transformation. Here, we have

k=3
r=.5
G=function(x) qnorm(pt(x,df=k))
dg=function(x) dt(x,df=k)/dnorm(qnorm(pt(x,df=k)))
Ginv=function(x) qt(pnorm(x),df=k)
S=matrix(c(1,r,r,1),2,2)
f=function(x,y) dmt(cbind(Ginv(x),Ginv(y)),S=S,df=k)/(dg(x)*dg(y))
vx=seq(-3,3,length=201)
vy=seq(-3,3,length=201)
z=outer(vx,vy,f)
set.seed(1)
Z=rmt(1500,S=S,df=k)
X=G(Z) Because we considered a nonlinear transformation of the margins, the level curves are no longer elliptical. But there is still some kind of symmetry.

• The Exchangeable Case with Conditionally Independent Random Variables

We did consider the case where  and  with independent random variables, given , and that both variables are exponentially distributed, with parameter . As we’ve seen in class, it might be difficult to visualize that sample, unless we have log scales on both axis. But instead of a log transformation, why not consider a transformation so that margins will be . The only technical problem is that we do not have the (nonconditional) distributions of the margins. Well, we have them, but they are integral based. From a computational point of view, that’s not a bit deal… Computations might take a while, but we can visualize the density using the following code (here, we assume that  is Gamma distributed)

a=.6
b=1
h=.0001
G=function(x) qnorm(ifelse(x<0,0,integrate(function(z) pexp(x,z)*
dgamma(z,a,b),lower=0,upper=Inf)$value)) Ginv=function(x) uniroot(function(z) G(z)-x,lower=-40,upper=1e5)$root
dg=function(x) (Ginv(x+h)-Ginv(x-h))/2/h
H=function(xy) integrate(function(z) dexp(xy,z)*dexp(xy,z)*
dgamma(z,a,b),lower=0,upper=Inf)$value f=function(x,y) H(c(Ginv(x),Ginv(y)))*(dg(x)*dg(y)) vx=seq(-3,3,length=151) vy=seq(-3,3,length=151) z=matrix(NA,length(vx),length(vy)) for(i in 1:length(vx)){ for(j in 1:length(vy)){ z[i,j]=f(vx[i],vy[j])}} set.seed(1) Theta=rgamma(1500,a,b) Z=cbind(rexp(1500,Theta),rexp(1500,Theta)) X=cbind(Vectorize(G)(Z[,1]),Vectorize(G)(Z[,2])) There is a small technical problem, but no big deal. Here, the joint distribution is quite different. Margins are – one more time – standard Gaussian, but the shape of the joint distribution is quite different, with an asymmetry from the lower (left) tail to the upper (right) tail. More details when we’ll introduce copulas. The only difference will be that the margins will be uniform on the unit interval, and not standard Gaussian. # Happy St Patrick’s Day I love Saint Patrick’s Day for, at least, two reasons. The first one is that, on March 17th, you can play out loud The Pogues, the second one is that it’s the only day in the year when I really enjoy getting a Guiness in a pub. And Guiness is important in statistical science (I did mention a couple of hours ago – on this blog – that beers were important for social reasons in the academic world, but that was for other reasons…) > theta=seq(0,pi/2,length=101) > leaf=sin(2*theta)+.25*sin(6*theta) > for(k in 0:3) + polygon(leaf*cos(theta+k*pi/2),leaf*sin(theta+k*pi/2),col="green") As mentioned in all my statistics and econometrics courses, the history of statistics (I mean here mathematical statistics) is closely related to Guinness. A long time ago, there was a Guinness Brewing Company of Dublin, which – as its name suggests – was an Irish brewing company. And the boss, who was to inherit the family business, decided to attract young students, trained in chemistry at Cambridge or Oxford. In 1899, William Sealy Gosset, who had obtained a double degree in math and chemistry, left Oxford to Dublin. And to be quite honest, being graduate in maths meant when he had studied differential equations and astronomy. Basically, mathematics were useless for Guinness, and he got there with his expertise in chemistry. In fact, William turned out to be also a very good administrator, but this has nothing to do with our story. William had good memories of his studies in math, and he wondered if he could find a problem to look at. He started studies on workmanship, noting that conditions vary so much (temperature, from hops, malt, manufacturing conditions …) that there were only few consistent data. The “law of errors” (the central limit theorem) can not apply under these conditions. In short, Bill (now we know each other a little, we’ll call him Bill) took many measurements, and noticed that the Poisson distribution could be an interesting model to work with. To make the story short, Bill managed to use statistical techniques to control the variance of the production, meaning that he was able to lower losses in the production of beer. A nice application like this one deserved publication in a scientific journal … Well, of course the Poisson distribution has long been known (it was 1904 and a few months before, Von Bortkiewicz found elegant applications of this law, as discussed in a post a few weeks ago). But there was a disclosure issue there: Bill’s contract prohibited him from disclosing secrets to the competitors. Meanwhile, Bill had met Karl Pearson, who was then editor of Biometrika, and encouraged him to publish his results. In 1906, Bill who had helped Guiness to gain a lot of money – doing applied mathematics can be usefull – managed to take a sabbatical to work with Pearson to Galton Laboratory biometrics. Bill and Karl decided to publish the work under a pseudonym “Student.” The legend claims that they had hesitated to use “pupil.” And for almost 30 years, “Mr Gosset” honorable employee Guinness led a dissolute life by publishing in statistical journals (after work in the brewery) always under the pseudonym “Student”. Of course, it might not be that simple. I mean, Bill had a family life, too. And his wife was the captain of the national Hockey team. So I hardly imagine Bill playing the smart ass and doing mathematical computations, when it was time to wash the dishes or iron his shirt… In 1908, he wrote a remarkable “the probable error of the mean” remarked, at least, by Ronald Fisher. In fact, Bill found that there was a interesting law, but – as the normal – it was difficult to manipulate to obtain confidence intervals. Without a computer, he had the idea of ​​using monte carlo methods to tabulate quantiles and construct its tables. And he was probably the first one to look carefully at the problem of small samples, unlike Karl Pearson, who always put focus on the asymptotic case. In fact, looking at his small sample, he saw the denominator magnitudes very close to those specifically manipulated Karl, in particular a square root of chi-square law. Well, of course, remained the normality assumption, but at least we had some results for finite samples ! For the story, William Gosset suggested to use letter z for its statistics, the ratio between the mean and (empirical) standard deviation. But a few years later, statisticians became accustomed to use this letter for Gaussian distribution (i.e. when the variance is known), and it became the standard to use the letter t. Hence finally the present name of “Student-t distribution” and in regression outputs, we have the “t-test”. A legend (told by Harold Hotelling in his memoirs) claims that the Guinness family discovered this double life on the day of the death of William Gosset in 1937 when mathematicians requested financial assistance to print a volume of the works of their employee. But another legend claims that Mr Guinness himself would have suggested his nickname when he had expressed his intention to publish his research… So I guess we’ll never know. But at least, I’ll think about Bill when I’ll get my first Guiness tonight (but I will probably not be able to tell this story anymore when I’ll reach the fourth…) # Test, valeur critique et p-value Un petit complément suite au cours de mercredi dernier, pour insister sur l’importance de la $p$-value dans la lecture de la sortie d’un test. • Les erreurs dans un test statistique Mais avant, rappelons qu’un test est une prise de décision: accepter ou rejeter une hypothèse. Et qu’on peut commettre une erreur. Ou pour être plus précis, on peut commettre deux types d’erreur, • accepter l’hypothèse alors que cette dernière est fausse • rejeter l’hypothèse alors que cette dernière était vraie Pour reprendre une terminologie plus médicale, un test de grossesse peut dire à une femme qu’elle n’est pas enceinte, alors qu’elle l’est; ou dire qu’elle l’est, alors qu’elle ne l’est pas (voir tous les exemples dans les exercices de probabilités de l’examen P de la SOA, ou le cours ACT2121). Formellement, on a deux probabilités, • la probabilité d’accepter à tort notre hypothèse (on parlera d’erreur de second espèce), $\beta$ • la probabilité de rejeter à tort notre hypothèse (on parlera d’erreur de première espèce) $\alpha$ Dans un monde idéal on voudrait que les deux probabilités soient aussi petites que possibles… Mais c’est impossible, et le plus souvent, baisser une des probabilités se fait en augmentant l’autre. Les cas extrêmes étant • avoir un test de grossesse qui déclare tout le monde enceinte: on ne rejette alors jamais à tort (on ne rejette jamais tout court en fait), mais on a un fort taux d’acceptation à tort, • avoir un test de grossesse qui ne déclare personne enceinte: on n’accepte jamais à tort (car on n’accepte jamais) mais on a un fort taux de rejet à tort. Bref, on a un arbitrage à faire entre deux types d’erreurs. Souvent, en pratique on va demander à contrôler l’erreur de première espèce (i.e. $\alpha$ de l’ordre de 5%), et on chercher a un test qui, à $\alpha$ donné, possède la plus faible erreur de première espèce. Voilà en gros pour la théorie: on se donne un seuil de significativité $\alpha$, qui correspond à la probabilité d’erreur de premier type. Et on va chercher à tester si une hypothèse $H_0$ est vraie, l’alternative étant une hypothèse $H_1$.  $H_0$ vraie $H_1$ vraie accepter $H_0$ OK erreur type 2 rejeter $H_0$ erreur type 1 OK • La “valeur critique” La notion de valeur critique a été introduite dans Neyman & Pearson (1928). Cette valeur dépend de la forme de l’hypothèse alternative, en particulier savoir si le test est bilatéral, unilatéral à gauche, ou unilatéral à droite. Pour un test donné, la valeur critique peut-être vue comme la valeur limite a partir de laquelle on pourra rejeter $H_0$ avec un seuil de significativité donné. • La p-value La $p$-value a été introduite dans Gibbons & Pratt (1975), meme si on peut retrouve l’idée beaucoup plus tôt, comme Pearson (1900), qui propose de calculer “the probability that the observed value of the chi-square statistic would be exceeded under the null hypothesis“. La $p$-value est la probabilité, sous $H_0$, d’obtenir une statistique aussi extrême (pour ne pas dire aussi grande) que la valeur observée sur l’échantillon. Aussi, pour un seuil de significativité $\alpha$ donné, on compare $p$ et $\alpha$, afin d’accepter, ou de rejeter $H_0$, • si $p\leq\alpha$, on va rejeter l’hypothèse $H_0$ (en faveur de $H_1$) • si $p>\alpha$, on va rejeter $H_1$ (en faveur de $H_0$). On peut alors interpréter la $p$-value comme le plus petit seuil de significativité pour lequel l’hypothèse nulle est acceptée. Gibbons & Pratt (1975) reviennent longuement sur les interprétations, et surtout les mauvaises interprétations, de cette $p$-value. • Valeur critique versus p-value Si on formalise un peu, on peut vouloir tester $H_0:\theta=\theta_0$ contre $H_1:\theta>theta_0$ (par exemple). De manière très générale, on dispose d’une statistique de test $T$ qui a pour loi, sous $H_0$, $F_{\theta_0}(\cdot)$ (que l’on supposera continue). Notons qu’on peut considérer une hypothèse alternative de la forme $H_1:\theta\neq\theta_0$, c’est juste plus pénible parce qu’il faut travailler sur $\vert T\vert$, et calculer des probabilités à gauche, ou à droite. Donc pour notre exemple, on va prendre un test unilatéral. Dans l’approche classique (telle que présentée dans tous les cours de statistiques), on se donne un seul d’acceptation $\alpha$ petit (disons 5%), et on cherche une valeur critique $T_{1-alpha}$ telle que Pour ceux qui se souviennent de leur cours de stats, cela peut faire penser à la puissance du test, définie par $\pi(\theta\vert \alpha)=\mathbb{P}(T\geq T_{1-\alpha}\vert \theta)=1-F_{\theta}(T_{1-\alpha})$ Formellement, la $p$-value associée au test $T$ est la variable aléatoire $P$ définie par $P=1-F_{\theta_0}(T)$. Donc effectivement, la $p$-value et la puissance sont liées, puisque $\mathbb{P}(P\leq \alpha\vert \theta)=\pi(\theta\vert \alpha)$ autrement dit, la puissance peut-être vue comme la fonction de répartition de la $p$-value. • Intérêt computationnel de la p-value D’un point de vue computationnel, la $p$-value est l’outil le plus important pour interpréter la sortie d’un test. Commençons par un test simple, comme une comparaison de moyennes. On cherche ici à tester $H_0:\mu_X=\mu_Y$ contre $H_1:\mu_X>\mu_Y$ pour des moyennes calculées sur deux groupes. Pour reprendre l’exemple abordé dans un précédant billet, on a les notes obtenues en ACT6420 par deux groupes différents. Et on veut savoir s’ils sont vraiment différents (ci-dessous le nombre de bonnes réponses, sur 40 questions, on travaillera ensuite sur la note sur 100) image manquante La statistique de test est ici $T = \frac{\overline{X} - \overline{Y}}{\displaystyle{ \sqrt{ {s_X^2 \over n_X} + {s_Y^2 \over n_Y} }}}$ et sous $H_0$, $T$ va suivre une loi de Student à $\nu$ degrés de liberté, où $\nu$ est donné par la relation de Welch–Satterthwaite (d’après Satterwaite (1946) et Welch (1947)), $\nu = {{\left( {s_X^2 \over n_X} + {s_Y^2 \over n_Y}\right)^2 } \over {{s_X^4 \over n_X^2 \cdot \left({n_X-1}\right)}+{s_Y^4 \over n_Y^2 \cdot \left({n_Y-1}\right)}}}$ Numériquement, on a ici > Xbar=mean(X) > Ybar=mean(Y) > Sx2=var(X) > Sy2=var(Y) > nX=length(X) > nY=length(Y) > (T=(Xbar-Ybar)/sqrt(Sx2/nX+Sy2/nY))  -2.155754 et pour les degrés de liberté > (nu=(Sx2/nX+Sy2/nY)^2/(Sx2^2/nX^2/(nX-1)+ + Sy2^2/nY^2/(nY-1)))  36.35279 La valeur critique est obtenue en lisant dans les tables, (car ici on a des probabilité pour un test bilatéral dans la table) comme on apprenait dans les cours de statistique au siècle passé. D’un point de vue informatique, on cherche à savoir si on est à gauche, ou à droite de la valeur critique > qt(.05,df=nu)  -1.687865 image manquante On peut aussi calculer la $p$-value, > pt(T,df=nu)  0.01889768 Si on regarde, sous R, il existe des fonctions de tests, pour comparer des moyennes. Et dans ce cas, la sortie est > t.test(X,Y,alternative = "less") Welch Two Sample t-test data: X and Y t = -2.1558, df = 36.353, p-value = 0.0189 alternative hypothesis: true difference in means is less than 0 95 percent confidence interval: -Inf -1.772507 sample estimates: mean of x mean of y 48.75000 56.91667 Autrement dit, on a automatiquement la $p$-value, et qui permet rapidement d’interpréter le test. Moralité, si on sait interpréter une $p$-value (et que l’on vérifié au préalable les conditions d’application d’un test), on peut faire tous les tests que l’on veut ! Si on veut faire un peu plus compliqué, on peut regarder la distribution des notes, et se demander si une loi $\mathcal{N}(60,15^2)$ serait possible (par exemple, ça sera notre hypothèse $H_0$, l’hypothèse alternative étant que ce n’est pas cette loi). Pour faire ce test, il existe le test de Kolmogorov-Smirnov. La statistique de test est ici $T=\sup\{\vert \widehat{F}_n(x)-F_0(x)\vert ,x\in\mathbb{R}\}$ $F_0(\cdot)$ est la fonction de répartition de la loi $\mathcal{N}(60,15^2)$, et $\widehat{F}_n(\cdot)$ est la fonction de répartition empirique $\widehat{F}_n(x)=\frac{1}{n}\sum_{i=1}^n \mathbf{1}(x_i\leq x)$ La loi de $T$ n’est pas simple, ou moins simple qu’une loi de Student (cf Marsaglia, Tsang & Wang (2003) par exemple). En revanche, on a les $p$-values automatiquement, > ks.test(Y, "pnorm", 60, 15) One-sample Kolmogorov-Smirnov test data: Y D = 0.1421, p-value = 0.5796 alternative hypothesis: two-sided Aussi, on peut accepter ici l’hypothèse nulle. On peut d’ailleurs faire un petit dessin pour s’en convaincre, > Femp=function(x) mean(Y<=x) > plot(0:100,Vectorize(Femp)(0:100),type="s") > lines(0:100,pnorm(0:100,60,15),col="red") image manquante Et ça va nous servir dans ce cours ? A priori oui… parce qu’on parlera du test de Student (pour tester si une variable dans une régression est significative), du test de Fisher (pour tester si plusieurs variables dans une régression sont significatives, ou plus généralement si une contrainte – linéaire – sur les coefficients peut être acceptée), du test de Chow (pour tester des ruptures dans un modèle linéaire, mais c’est un test de Fisher un peu déguisé), du test d’Anderson-Darling (pour tester si des résidus sont Gaussiens), du test de Breuch-Pagan voire le test de White (pour tester si les résidus peuvent être considérés de variance constante), du test de Durbin-Watson (pour tester s’il n’y a pas d’auto-corrélation dans la série des résidus), du test de Dickey-Fuller (pour tester si une série temporelle est – ou n’est pas – stationnaire), des tests de Franses (pour tester si une série peut être considérée comme saisonnière, ou pas), du test de Ljung-Box (pour tester si un bruit est un bruit blanc)… Et j’en oublie un paquet. Donc quand il est dit (dans le plan de cours) que le cours de statistique est un prérequis, il ne s’agit pas de l’avoir suivi, mais bel et bien de l’avoir compris, car on passera notre temps à utiliser des notions entrevues dans ce cours. # (nonparametric) copula density estimation Today, we will go further on the inference of copula functions. Some codes (and references) can be found on a previous post, on nonparametric estimators of copula densities (among other related things). Consider (as before) the loss-ALAE dataset (since we’ve been working a lot on that dataset) > library(MASS) > library(evd) > X=lossalae > U=cbind(rank(X[,1])/(nrow(X)+1),rank(X[,2])/(nrow(X)+1)) The standard tool to plot nonparametric estimators of densities is to use multivariate kernels. We can look at the density using > mat1=kde2d(U[,1],U[,2],n=35) > persp(mat1$x,mat1$y,mat1$z,col="green",
+ xlab="",ylab="",zlab="",zlim=c(0,7))

or level curves (isodensity curves) with more detailed estimators (on grids with shorter steps)

> mat1=kde2d(U[,1],U[,2],n=101)
> image(mat1$x,mat1$y,mat1$z,col= + rev(heat.colors(100)),xlab="",ylab="") > contour(mat1$x,mat1$y,mat1$z,add=
+ TRUE,levels = pretty(c(0,4), 11)) Kernels are nice, but we clearly observe some border bias, extremely strong in corners (the estimator is 1/4th of what it should be, see another post for more details). Instead of working on sample  on the unit square, consider some transformed sample , where  is a given function. E.g. a quantile function of an unbounded distribution, for instance the quantile function of the  distribution. Then, we can estimate the density of the transformed sample, and using the inversion technique, derive an estimator of the density of the initial sample. Since the inverse of a (general) function is not that simple to compute, the code might be a bit slow. But it does work,

> gaussian.kernel.copula.surface <- function (u,v,n) {
+   s=seq(1/(n+1), length=n, by=1/(n+1))
+   mat=matrix(NA,nrow = n, ncol = n)
+ sur=kde2d(qnorm(u),qnorm(v),n=1000,
+ lims = c(-4, 4, -4, 4))
+ su<-sur$z + for (i in 1:n) { + for (j in 1:n) { + Xi<-round((qnorm(s[i])+4)*1000/8)+1; + Yj<-round((qnorm(s[j])+4)*1000/8)+1 + mat[i,j]<-su[Xi,Yj]/(dnorm(qnorm(s[i]))* + dnorm(qnorm(s[j]))) + } + } + return(list(x=s,y=s,z=data.matrix(mat))) + } Here, we get Note that it is possible to consider another transformation, e.g. the quantile function of a Student-t distribution. > student.kernel.copula.surface = + function (u,v,n,d=4) { + s <- seq(1/(n+1), length=n, by=1/(n+1)) + mat <- matrix(NA,nrow = n, ncol = n) + sur<-kde2d(qt(u,df=d),qt(v,df=d),n=5000, + lims = c(-8, 8, -8, 8)) + su<-sur$z
+ for (i in 1:n) {
+     for (j in 1:n) {
+ 	Xi<-round((qt(s[i],df=d)+8)*5000/16)+1;
+ 	Yj<-round((qt(s[j],df=d)+8)*5000/16)+1
+ 	mat[i,j]<-su[Xi,Yj]/(dt(qt(s[i],df=d),df=d)*
+ 	dt(qt(s[j],df=d),df=d))
+     }
+ }
+ return(list(x=s,y=s,z=data.matrix(mat)))
+ }

Another strategy is to consider kernel that have precisely the unit interval as support. The idea is here to consider the product of Beta kernels, where parameters depend on the location

> beta.kernel.copula.surface=
+  function (u,v,bx=.025,by=.025,n) {
+  s <- seq(1/(n+1), length=n, by=1/(n+1))
+  mat <- matrix(0,nrow = n, ncol = n)
+ for (i in 1:n) {
+     a <- s[i]
+     for (j in 1:n) {
+     b <- s[j]
+ 	mat[i,j] <- sum(dbeta(a,u/bx,(1-u)/bx) *
+     dbeta(b,v/by,(1-v)/by)) / length(u)
+     }
+ }
+ return(list(x=s,y=s,z=data.matrix(mat)))
+ } On those two graphs, we can clearly observe strong tail dependence in the upper (right) corner, that cannot be intuited using a standard kernel estimator…

# Copulas estimation and influence of margins

Just a short post to get back on results mentioned at the end of the course. Since copulas are obtained using (univariate) quantile functions in the joint cumulative distribution function, they are – somehow – related to the marginal distribution fitted. In order to illustrate this point, consider an i.i.d. sample from a Student-t distribution,

library(mnormt)
r=.5
n=200
X=rmt(n,mean=c(0,0),S=matrix(c(1,r,r,1),2,2),df=4)

Thus, the true copula is Student-t. Here, with 4 degrees of freedom. Note that we can easily get the (true) value of the copula, on the diagonal

dg=function(t) pmt(qt(t,df=4),mean=c(0,0),
S=matrix(c(1,r,r,1),2,2),df=4)
DG=Vectorize(dg)

Four strategies are considered here to define pseudo-copula base variates,

• misfit: consider an invalid marginal estimation: we have assumed that margins were Gaussian, i.e. • perfect fit: here, we know that margins were Student-t, with 4 degrees of freedom • standard fit: then, consider the case where we fit marginal distribution, but in the good family this time (e.g. among Student-t distributions), • ranks: finally, we consider nonparametric estimators for marginal distributions, Now that we have a sample with margins in the unit square, let us construct the empirical copula, Let us now compare those four approaches.

• The first one is to illustrate model error, i.e. what’s going on if we fit distributions, but not in the proper family of parametric distributions.
X0=cbind((X[,1]-mean(X[,1])/sd(X[,1])),
(X[,2]-mean(X[,2])/sd(X[,2])))
Y=pnorm(X0)

Then, the following code is used to compute the value of the empirical copula, on the diagonal,

diagonale=function(t,Z) mean((Z[,1]<=t)&(Z[,2]<=t))
diagY=function(t) diagonale(t,Y)
DiagY=Vectorize(diagY)
u=seq(0,1,by=.005)
dY=DiagY(u)

On the graph below, 1,000 samples of size 200 have been generated. All trajectories are the estimation of the copula on the diagonal. The black plain line is the true value of the copula Obviously, it is not good at all. Mainly because the distribution of can’t be a copula, since margins are not even uniform on the unit interval.

• a perfect fit. Here, we use the following code to generate our copula-type sample
U=pt(X,df=4)

This time, the fit is much better. • Using maximum likelihood estimators to fit the best distribution within the Student-t family
F1=fitdistr(X0[,1],dt,list(df=5),lower = 0.001)
F2=fitdistr(X0[,2],dt,list(df=5),lower = 0.001)
V=cbind(pt(X0[,1],df=F1$estimate),pt(X0[,2],df=F2$estimate))

Here, it is also very good. Even better than before, when the true distribution is considered. (it is like using Lillie test for goodness of fit, versus Kolmogorov-Smirnov, seehere for instance, in French).

• Finally, let us consider ranks, or nonparametric estimators for marginal distributions,
R=cbind(rank(X[,1])/(n+1),rank(X[,2])/(n+1))

Here it is even better then the previous one If we compare Box-plots of the value of the copula at point (.2,.2), we obtain the following, with on top ranks, then fitting with the good family, then using the true distribution, and finally, using a non-proper distribution. Just to illustrate one more time a result mentioned in a previous post, “in statistics, having too much information might not be a good thing“.

# Does the Student based confidence interval have any interest in practice ?

Friday in the course of statistics, we started the section on confidence interval, and like always, I got a bit confused with the degrees of freedom of the Student (should it be or ?) and which empirical variance (should we consider the one where we divide by or the one with ?).
And each time I start to get confused, the student obviously see it, and start to ask tricky questions… So let us make it clear now. The correct formula is the following: let then is a confidence interval for the mean of a Gaussian i.i.d. sample.
But the important thing is neither the n-1 that appear as degrees of freedom nor the that appear in the estimation of the standard error. Like always in mathematical result, the most important part of that result is not mentioned here: observations have to be i.i.d. and to be normally distributed. And not “almost” normally distributed….
Consider the following case: we have =20 observations that are almost normally distributed. Hence, I consider a student t distribution

n=20; X=rt(n,df=3)

An Anderson Darling normality test accepts a normal distribution in 2 cases out of 3.

for(s in 1:10000){
X=rt(n,df=3)
}
mean(pv>.05)
 0.6799

With a true normal distribution if would be 95% of the cases, so in some sense, I can pretend that I generate almost normal samples.
For those samples, we can look at bounds of the 90% confidence interval for the mean, with three different formulas, i.e. the correct one, or the one where I considered degrees of freedom instead of , and the one were we condired a Gaussian quantile instead of a Student t one, (and one might think to look at the non-unbiased estimator of the variance, also).
for(s in 1:10000){
X=rt(n,df=3)
m[s]=mean(X)
sd=sqrt(var(X))
IC1[s]=m[s]-qt(.95,df=n-1)*sd/sqrt(n)
IC2[s]=m[s]-qt(.95,df=n)*sd/sqrt(n)
IC3[s]=m[s]-qnorm(.95)*sd/sqrt(n)
}

One the graph below are plotted the distributions of the values obtained as lower bound of the 90% confidence interval, (the curves with and degrees of freedom in quantiles are the same, here).
The dotted vertical line is the true lower bound of the 90%-confidence interval, given the true distribution (which was not a Gaussian one).
If I get back to the standard procedure in any statistical textbook, since the sample is almost Gaussian, the lower bound of the confidence interval should be (since we have a Student t distribution)

mean(IC1)
 -0.605381

mean(IC3)
 -0.5759391

(obtained with a Gaussian distribution instead of a Student one). Actually, both of them are quite different from the correct one which was

quantile(m,.05)
5%
-0.623578

As I mentioned in a previous post (here), an important issue is that if we do not know a parameter and substitute an estimator, there is usually a cost (which means usually that the confidence interval should be larger). And this is what we observe here. From a teacher’s point of view, it is an important issue that should be mentioned in statistical courses….

But another important point is also that confidence interval is valid only if the underlying distribution is Gaussian. And not almost Gaussian, but really a Gaussian one.  So since with =20 observations everything might look Gaussian, I was wondering what should be done in practice… Because in some sense, using a Student quantile based confidence interval on some almost Gaussian sample is as wrong as using a Gaussian quantile based confidence interval on some Gaussian sample…

# In statistics, having too much information might not be a good thing

A common idea in statistics is that if we don’t know something, and we use anestimator of that something (instead of the true value) then there will be some additional uncertainty. For instance, consider a random sample, i.i.d., from a Gaussian distribution. Then, a confidence interval for the mean is where is the quantile of probability level of the standard normal distribution . But usually, standard deviation (the something is was talking about earlier) is usually unknown. So we substitute an estimation of the standard deviation, e.g. and the cost we have to pay is that the new confidence interval is where now is the quantile of the Student distribution, of probability level , with degrees of freedom.
We call it a cost since the new confidence interval is now larger (the Student distribution has higher upper-quantiles than the Gaussian distribution).
So usually, if we substitute an estimation to the true value, there is a price to pay.
A few years ago, with Jean David Fermanian and Olivier Scaillet, we were writing a survey on copula density estimation (using kernels,  here). At the end, we wanted to add a small paragraph on the fact that we assumed that we wanted to fit a copula on a sample i.i.d. with distribution , a copula, but in practice, we start from a sample with joint distribution (assumed to have continuous margins, and – unique – copula ). But since margins are usually unknown, there should be a price for not observing them.
To be more formal, in a perfect wold, we would consider but in the real world, we have to consider where it is standard to consider ranks, i.e. are empirical cumulative distribution functions.
My point is that when I ran simulations for the survey (the idea was more to give illustrations of several techniques of estimation, rather than proofs of technical theorems) we observed that the price to pay… was negative ! I.e. the variance of the estimator of the density (wherever on the unit square) was smaller on the pseudo sample than on perfect sample .
By that time, we could not understand why we got that counter-intuitive result: even if we do know the true distribution, it is better not to use it, and to use instead a nonparametric estimator. Our interpretation was based on the discrepancy concept and was related to the latin hypercube construction: With ranks, the data are more regular, and marginal distributions are exactlyuniform on the unit interval. So there is less variance.
This was our heuristic interpretation.
A couple of weeks ago, Christian Genest and Johan Segers proved that intuition in an article published in JMVA, Well, we observed something for finite , but Christian and Johan obtained an analytical result. Hence, if we denote the empirical copula in the perfect world (with known margins) and the one constructed from the pseudo sample, they obtained that, everywhere with nice graphs of ,  So I was very happy last week when Christian show me their results, to learn that our intuition was correct. Nevertheless, it is still a very counter-intuitive result…. If anyone has seen similar things, I’d be glad to hear about it !

# Happy birthday Arthur !

mais non, ce n’est pas mon anniversaire ! Mais l’occasion était trop belle… oui, c’est l’anniversaire d’Arthur Guinness (Arth Guinness pour ceux qui ont déjà eu une canette entre les mains…) Et comme je l’ai rappelé lors des rappels de statistiques, l’histoire de la statistique (j’entends par là statistique mathématique) est étroitement liée à Guiness.

Il y fort fort longtemps, il existait une Guiness Brewing Company of Dublin, qui comme son nom l’indique, était une brasserie irlandaise. Et au début du XXème siècle, le propriétaire, qui venait d’hériter de cette entreprise familiale décida d’attirer des jeunes étudiants, formés en chimie à Cambridge, ou Oxford.

En 1899, un  certain William Sealy Gosset, qui avait obtenu un double diplôme, en maths et en chimie, quitta Oxford pour Dublin. Et pour être tout à fait honnête, diplômé en maths signifiait alors qu’il avait étudié les équations différentielles ou l’astronomie. En gros, les maths, ça ne servait à rien chez Guiness, et il était arrivé là grâce à ses compétences en chimie. Et en fait, William s’avérait être en plus un très bon administrateur, mais on sort un peu de notre sujet.

Mais William avait gardé un bon souvenir de ses études en maths, et il se demandait s’il pouvait trouver un problème sur lequel se pencher…. il commença a travailler des études de qualité de fabrication, notant que les conditions varient tellement (température, provenance du houblon, du malt, conditions de fabrication…) que les données homogènes sont peu nombreuses. La “loi des erreurs” classique ne peut pas s’appliquer dans ces conditions, i.e. le théorème central limite.

Bref, Bill (maintenant qu’on se connaît un peu, on va l’appeler Bill) pris beaucoup de mesures, et observa une loi de Poisson, ce qui était pratique car cette loi n’avait qu’un paramètre. Bref, Bill arriva à contrôler la variance, et donc à faire moins de pertes lors de la production de bière…

Une jolie application comme celle là méritait publication dans une revue scientifique… Bon, bien sûr la loi de Poisson était connue depuis longtemps (on était en 1904 et quelques mois auparavant, Von Bortkiewicz avait trouvé d’élégantes applications à cette loi, en particulier sur les accidents de chevaux dans l’armée Prusse). Et à côté, il y avait un soucis de taille que connaissent bien des actuaires travaillant dans les compagnies d’assurance: le contrat de Bill lui interdisait de divulguer des secrets à la concurrence, qui plus est des secrets capables de leur faire gagner de l’argent !

Entre temps, Bill avait rencontré Karl Pearson, qui était alors éditeur de la revue Biometrika, et qui l’encouragea à publier ses résultats. En 1906, Bill qui avait fait gagné beaucoup d’argent à Guiness – et qui avait d’ailleurs montré l’intérêt de faire des maths appliqués – réussi à prendre une année sabbatique pour aller travailler avec Pearson au laboratoire de biométrie de Galton. Et Bill et Karl décidèrent de faire publier les travaux sous un pseudonyme, “Student”. La légende prétend qu’ils avaient aussi hésité à prendre “pupil“.

Et pendant presque 30 ans, “Mr Gosset“, honorable employé de Guinness menait une vie dissolue en publiant dans des revues de statistiques (après son travail dans la brasserie1), toujours sous le pseudonyme de “Student“.

En 1908, on lui doit un remarquable “the probable error of the mean“, qui remarque en tous les cas Ronald Fisher. En fait, Bill trouva qu’il existait une loi assez remarquable, mais qui – comme la loi normale – était délicate à manipuler pour obtenir des intervalles de confiance. Sans ordinateur, il eu malgré tout l’idée d’utiliser des méthodes de monte carlo pour tabuler les quantiles, et construire ses tables. Et il fut un des premiers à se pencher sur le problème des petits échantillons, contrairement à Karl Pearson, qui se plaçait toujours dans le cas asymptotique.

En fait, en regardant son petit échantillon, il vit apparaître au dénominateur des grandeurs très proches de celles que manipulait précisément Karl, en particulier une racine carré de loi du chi-deux. Bon, bien sûr, restait l’hypothèse (fondamentale) de normalité des observations, mais au moins, on avait des résultats à distance finie !

Pour la petite histoire, William Gosset proposa la lettre $z$ pour sa statistique, le ratio entre la moyenne et l’écart-type. Mais quelques années plus tard, les statisticiens prirent l’habitude d’utiliser cette lettre pour des lois normales (autrement dit quand la variance est connue), et la norme devint d’utiliser la lettre $t$. D’où finalement la dénomination actuelle de “Student-t distribution”.

Une légende (racontée par Harold Hotelling dans ses mémoires) veut que la famille Guinness découvrit cette double vie le jour de la mort de William Gosset, en 1937, lorsque des mathématiciens demandèrent une aide financière pour imprimer un volume des œuvres de leur employé. Mais une autre légende raconte que c’est Guinness qui lui aurait suggéré son pseudo lorsqu’il aurait fait part de ses intentions de publier ses travaux de recherche…

1 et si si, Bill avait une vie de famille… il avait d’ailleurs une femme qui devait le mener à la baguette, puisqu’elle était à l’époque capitaine de l’équipe nationale de Hockey ! il devait pas faire le malin avec ses calculs de loi quand c’était l’heure de faire la vaisselle….