# On Cochran Theorem (and Orthogonal Projections)

Cochran Theorem – from The distribution of quadratic forms in a normal system, with applications to the analysis of covariance published in 1934 – is probably the most import one in a regression course. It is an application of a nice result on quadratic forms of Gaussian vectors. More precisely, we can prove that if $\boldsymbol{Y}\sim\mathcal{N}(\boldsymbol{0},\mathbb{I}_d)$ is a random vector with $d$ $\mathcal{N}(0,1)$ variable then (i) if $A$ is a (squared) idempotent matrix $\boldsymbol{Y}^\top A\boldsymbol{Y}\sim\chi^2_r$ where $r$ is the rank of matrix $A$, and (ii) conversely, if $\boldsymbol{Y}^\top A\boldsymbol{Y}\sim\chi^2_r$ then $A$ is an idempotent matrix of rank $r$. And just in case, $A$ is an idempotent matrix means that $A^2=A$, and a lot of results can be derived (for instance on the eigenvalues). The prof of that result (at least the (i) part) is nice: we diagonlize matrix $A$, so that $A=P\Delta P^\top$, with $P$ orthonormal. Since $A$ is an idempotent matrix observe that$$A^2=P\Delta P^\top=P\Delta P^\top=P\Delta^2 P^\top$$where $\Delta$ is some diagonal matrix such that $\Delta^2=\Delta$, so terms on the diagonal of $\Delta$ are either $0$ or $1$‘s. And because the rank of $A$ (and $\Delta$) is $r$ then there should be $r$ $1$‘s and $d-r$ $1$‘s. Now write$$\boldsymbol{Y}^\top A\boldsymbol{Y}=\boldsymbol{Y}^\top P\Delta P^\top\boldsymbol{Y}=\boldsymbol{Z}^\top \Delta\boldsymbol{Z}$$where $\boldsymbol{Z}=P^\top\boldsymbol{Y}$ that satisfies$\boldsymbol{Z}\sim\mathcal{N}(\boldsymbol{0},PP^\top)$ i.e. $\boldsymbol{Z}\sim\mathcal{N}(\boldsymbol{0},\mathbb{I}_d)$. Thus $$\boldsymbol{Z}^\top \Delta\boldsymbol{Z}=\sum_{i:\Delta_{i,i}-1}Z_i^2\sim\chi^2_r$$Nice, isn’t it. And there is more (that will be strongly connected actually to Cochran theorem). Let $A=A_1+\dots+A_k$, then the two following statements are equivalent (i) $A$ is idempotent and $\text{rank}(A)=\text{rank}(A_1)+\dots+\text{rank}(A_k)$ (ii) $A_i$‘s are idempotents, $A_iA_j=0$ for all $i\neq j$.

Now, let us talk about projections. Let $\boldsymbol{y}$ be a vector in $\mathbb{R}^n$. Its projection on the space $\mathcal V(\boldsymbol{v}_1,\dots,\boldsymbol{v}_p)$ (generated by those $p$ vectors) is the vector $\hat{\boldsymbol{y}}=\boldsymbol{V} \hat{\boldsymbol{a}}$ that minimizes $\|\boldsymbol{y} -\boldsymbol{V} \boldsymbol{a}\|$ (in $\boldsymbol{a}$). The solution is$$\hat{\boldsymbol{a}}=( \boldsymbol{V}^\top \boldsymbol{V})^{-1} \boldsymbol{V}^\top \boldsymbol{y} \text{ and } \hat{\boldsymbol{y}} = \boldsymbol{V} \hat{\boldsymbol{a}}$$
Matrix $P=\boldsymbol{V} ( \boldsymbol{V}^\top \boldsymbol{V})^{-1} \boldsymbol{V}^\top$ is the orthogonal projection on $\{\boldsymbol{v}_1,\dots,\boldsymbol{v}_p\}$ and $\hat{\boldsymbol{y}} = P\boldsymbol{y}$.

Now we can recall Cochran theorem. Let $\boldsymbol{Y}\sim\mathcal{N}(\boldsymbol{\mu},\sigma^2\mathbb{I}_d)$ for some $\sigma>0$ and $\boldsymbol{\mu}$. Consider sub-vector orthogonal spaces $F_1,\dots,F_m$, with dimension $d_i$. Let $P_{F_i}$ be the orthogonal projection matrix on $F_i$, then (i) vectors $P_{F_1}\boldsymbol{X},\dots,P_{F_m}\boldsymbol{X}$ are independent, with respective distribution $\mathcal{N}(P_{F_i}\boldsymbol{\mu},\sigma^2\mathbb{I}_{d_i})$ and (ii) random variables $\|P_{F_i}(\boldsymbol{X}-\boldsymbol{\mu})\|^2/\sigma^2$ are independent and $\chi^2_{d_i}$ distributed.

We can try to visualize those results. For instance, the orthogonal projection of a random vector has a Gaussian distribution. Consider a two-dimensional Gaussian vector

library(mnormt) r = .7 s1 = 1 s2 = 1 Sig = matrix(c(s1^2,r*s1*s2,r*s1*s2,s2^2),2,2) Sig Y = rmnorm(n = 1000,mean=c(0,0),varcov = Sig) plot(Y,cex=.6) vu = seq(-4,4,length=101) vz = outer(vu,vu,function (x,y) dmnorm(cbind(x,y), mean=c(0,0), varcov = Sig)) contour(vu,vu,vz,add=TRUE,col='blue') abline(a=0,b=2,col="red")

Consider now the projection of points $\boldsymbol{y}=(y_1,y_2)$ on the straight linear with directional vector $\overrightarrow{\boldsymbol{u}}$ with slope $a$ (say $a=2$). To get the projected point $\boldsymbol{x}=(x_1,x_2)$ recall that $x_2=ay_1$ and $\overrightarrow{\boldsymbol{x},\boldsymbol{y}}\perp\overrightarrow{\boldsymbol{u}}$. Hence, the following code will give us the orthogonal projections

p = function(a){ x0=(Y[,1]+a*Y[,2])/(1+a^2) y0=a*x0 cbind(x0,y0) }

with

P = p(2) for(i in 1:20) segments(Y[i,1],Y[i,2],P[i,1],P[i,2],lwd=4,col="red") points(P[,1],P[,2],col="red",cex=.7)

Now, if we look at the distribution of points on that line, we get… a Gaussian distribution, as expected,

z = sqrt(P[,1]^2+P[,2]^2)*c(-1,+1)[(P[,1]>0)*1+1] vu = seq(-6,6,length=601) vv = dnorm(vu,mean(z),sd(z)) hist(z,probability = TRUE,breaks = seq(-4,4,by=.25)) lines(vu,vv,col="red")

Or course, we can use the matrix representation to get the projection on $\overrightarrow{\boldsymbol{u}}$, or a normalized version of that vector actually

a=2 U = c(1,a)/sqrt(a^2+1) U [1] 0.4472136 0.8944272 matP = U %*% solve(t(U) %*% U) %*% t(U) matP %*% Y[1,] [,1] [1,] -0.1120555 [2,] -0.2241110 P[1,] x0 y0 -0.1120555 -0.2241110 

(which is consistent with our manual computation). Now, in Cochran theorem, we start with independent random variables,

Y = rmnorm(n = 1000,mean=c(0,0),varcov = diag(c(1,1)))

Then we consider the projection on $\overrightarrow{\boldsymbol{u}}$ and $\overrightarrow{\boldsymbol{v}}=\overrightarrow{\boldsymbol{u}}^\perp$

U = c(1,a)/sqrt(a^2+1) matP1 = U %*% solve(t(U) %*% U) %*% t(U) P1 = Y %*% matP1 z1 = sqrt(P1[,1]^2+P1[,2]^2)*c(-1,+1)[(P1[,1]>0)*1+1] V = c(a,-1)/sqrt(a^2+1) matP2 = V %*% solve(t(V) %*% V) %*% t(V) P2 = Y %*% matP2 z2 = sqrt(P2[,1]^2+P2[,2]^2)*c(-1,+1)[(P2[,1]>0)*1+1]

We can plot those two projections

plot(z1,z2)

and observe that the two are indeed, independent Gaussian variables. And (of course) there squared norms are $\chi^2_{1}$ distributed.

# Probabilistic Foundations of Econometrics, part 1

In a series of posts, I wanted to get into details of the history and foundations of econometric and machine learning models. It will be some sort of online version of our joint paper with Emmanuel Flachaire and Antoine Ly, Econometrics and Machine Learning (initially writen in French), that will actually appear soon in the journal Economics and Statistics. This is the first one…

The importance of probabilistic models in economics is rooted in Working’s (1927) questions and the attempts to answer them in Tinbergen’s two volumes (1939). The latter have subsequently generated a great deal of work, as recalled by Duo (1993) in his book on the foundations of econometrics, and more particularly in the first chapter “The Probability Foundations of Econometrics”. It should be recalled that Trygve Haavelmo was awarded the Nobel Prize in Economics in 1989 for his “clarification of the foundations of the probabilistic theory of econometrics”. Because as Haavelmo (1944) (initiating a profound change in econometric theory in the 1930s, as recalled in Morgan’s Chapter 8 (1990)) showed, econometrics is fundamentally based on a probabilistic model, for two main reasons. First, the use of statistical quantities (or “measures”) such as means, standard errors and correlation coefficients for inferential purposes can only be justified if the process generating the data can be expressed in terms of a probabilistic model. Second, the probability approach is relatively general, and is particularly well suited to the analysis of “dependent” and “non-homogeneous” observations, as they are often found on economic data.We will then assume that there is a probabilistic space $(\Omega,\mathcal{F},\mathbb{P})$ such that observations $(y_i,\mathbf{x}_i)$ are seen as realizations of random variables $(Y_i, \mathbf{X}_i)$. In practice, however, we are not very interested in the joint law of the couple $(Y, \mathbf{X})$ : the law of $\mathbf{X}$ is unknown, and it is the law of Y conditional on $\mathbf{X}$ that will be interested in. In the following, we will note $x$ a single observation, $\mathbf{x}$ a vector of observations, $X$ a random variable, and $\mathbf{X}$ a random vector. Abusively, $\mathbf{X}$ may also designate the matrix of individual observations (denoted $\mathbf{x}_i$), depending on the context.

## Foundations of mathematical statistics

As recalled in Vapnik’s (1998) introduction, inference in parametric statistics is based on the following belief: the statistician knows the problem to be analyzed well, in particular, he knows the physical law that generates the stochastic properties of the data, and the function to be found is written via a finite number of parameters[1]. To find these parameters, the maximum likelihood method is used. The purpose of the theory is to justify this approach (by discovering and describing its favorable properties). We will see that in learning, philosophy is very different, since we do not have a priori reliable information on the statistical law underlying the problem, nor even on the function we would like to approach (we will then propose methods to construct an approximation from the data at our disposal, as in (1998)). A “golden age” of parametric inference, from 1930 to 1960, laid the foundations for mathematical statistics, which can be found in all statistical textbooks, including today. As Vapnik (1998) states, the classical parametric paradigm is based on the following three beliefs:

1. To find a functional relationship from the data, the statistician is able to define a set of functions, linear in their parameters, that contain a good approximation of the desired function. The number of parameters describing this set is small.
2. The statistical law underlying the stochastic component of most real-life problems is the normal law. This belief has been supported by reference to the central limit theorem, which stipulates that under large conditions the sum of a large number of random variables is approximated by the normal law.
3. The maximum likelihood method is a good tool for estimating parameters.

In this section we will come back to the construction of the econometric paradigm, directly inspired by that of classical inferential statistics.

## Conditional laws and likelihood

Linear econometrics has been constructed under the assumption of individual data, which amounts to assuming independent variables $(Y_i, \mathbf{X}_i)$ (if it is possible to imagine temporal observations – then we would have a process $(Y_t, \mathbf{X}_t)$ – but we will not discuss time series here). More precisely, we will assume that, conditionally to the explanatory variables $\mathbf{X}_i$, the variables $Y_i$ are independent. We will also assume that these conditional laws remain in the same parametric family, but that the parameter is a function of $\mathbf{x}$. In the Gaussian linear model it is assumed that: $$(Y\vert \mathbf{X}=\mathbf{x})\overset{\mathcal{L}}{\sim}\mathcal{N}(\mu(\mathbf{x}),\sigma^2)~~~~ (1)$$where $\mu(\mathbf{x})=\beta_0+\mathbf{x}^T\mathbf{\beta}$ and $\mathbf{\beta}\in\mathbb{R}^{p}$.

It is usually called a ‘linear’ model since $\mathbb{E}[Y\vert \mathbf{X}=\mathbf{x}]=\beta_0+\mathbf{x}^T\mathbf{\beta}$ is a linear combination of covariates[2]. It is said to be a homoscedastic model if $Var[Y|\mathbf{X}=\mathbf{x}]=\sigma^2$, where $\sigma^2$ is a positive constant. To estimate the parameters, the traditional approach is to use the Maximum Likelihood estimator, as initially suggested by Ronald Fisher. In the case of the Gaussian linear model, log-likelihood is written:  $$\log\mathcal{L}(\beta_0, \mathbf{\beta},\sigma^2\vert \mathbf{y},\mathbf{x}) = -\frac{n}{2}\log[2\pi\sigma^2] - \frac{1}{2\sigma^2}\sum_{i=1}^n (y_i-\beta_0-\mathbf{x}_i^T\mathbf{\beta})^2$$Note that the term on the right, measuring a distance between the data and the model, will be interpreted as deviance in generalized linear models. Then we will set: $$(\widehat{\beta}_0,\widehat{\mathbf{\beta}},\widehat{\sigma}^2)=\text{argmax}\left\lbrace\log\mathcal{L}(\beta_0, \mathbf{\beta},\sigma^2\vert \mathbf{y},\mathbf{x})\right\rbrace$$The maximum likelihood estimator is obtained by minimizing the sum of the error squares (the so-called “least squares” estimator) that we will find in the “machine learning” approach.

The first order conditions allow to find the normal equations, whose matrix writing is $\mathbf{X}^T[\mathbf{y}-\mathbf{X}\mathbf{\beta}]=\mathbf{0}$, which can also be written $(\mathbf{X}^T \mathbf{X})\mathbf{\beta}=\mathbf{X}^T \mathbf{y}$. If $\mathbf{X}$ is a full (column) rank matrix, then we find the classical estimator:$$\widehat{\mathbf{\beta}}=(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{y}=\mathbf{\beta}+(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^{-1}\mathbf{\varepsilon}~~~(2)$$using residual-based writing (as often in econometrics), $y=\mathbf{x}^T\mathbf{\beta}+\varepsilon$. Gauss Markov’s theorem ensures that this estimator is the unbiased linear estimator with minimum variance. It can then be shown that $\widehat{\mathbf{\beta}}\sim\mathcal{N}(\mathbf{\beta},\sigma^2(\mathbf{X}^T\mathbf{X})^{-1})$, and in particular, if we simply need the first two moments : $$\mathbb{E}[\widehat{\mathbf{\beta}}]=\mathbf{\beta}~~~Var[\widehat{\mathbf{\beta}}]=\sigma^2 [\mathbf{X}^T\mathbf{X}]^{-1}$$In fact, the normality hypothesis makes it possible to make a link with mathematical statistics, but it is possible to construct this estimator given by equation (2) without that Gaussian assumption. Hence, if we assume that $Y|\mathbf{X}$ has the same distribution as $\mathbf{x}^T\mathbf{\beta}+\varepsilon$, where $\mathbb{E}[\varepsilon]=0$, $Var[\varepsilon]=\sigma^2$ and $Cov[X_j,\varepsilon]=0$ for all $j$, then $\widehat{\mathbf{\beta}}$ is an unbiased estimator of $\mathbf{\beta}$ with smallest variance[3] among unbiased linear estimators. Furthermore, if we cannot get normality at finite distance, asymptotically this estimator is Gaussian, with $$\sqrt{n}(\widehat{\mathbf{\beta}}-\mathbf{\beta})\overset{\mathcal{L}}{\rightarrow}\mathcal{N}(\mathbf{0},\mathbf{\Sigma})$$as $n\rightarrow\infty$, for some matrix $\mathbf{\Sigma}$.
The condition of having a full rank $\mathbf{X}$ matrix can be (numerically) strong in large dimensions. If it is not satisfied, $(\mathbf{X}^T \mathbf{X})^{-1}\mathbf{X}^T$ does not exist. If $\mathbb{I}$ denotes the identity matrix, however, it should be noted that $(\mathbf{X}^T \mathbf{X}+\lambda\mathbb{I})^{-1}\mathbf{X}^T$ still exists, whatever $\lambda>0$. This estimator is called the ridge estimator of level \lambda (introduced in the 1960s by Hoerl (1962), and associated with a regularization studied by Tikhonov (1963)). This estimator naturally appears in a Bayesian econometric context.

## Residuals

It is not uncommon to introduce the linear model from the distribution of the residuals, as we mentioned earlier. Also, equation (1) is written as often: $$y_i=\beta_0+\mathbf{x}_i^T\mathbf{\beta}+\varepsilon_i~~~~(3)$$where $\varepsilon_i$’s are realizations of independent and identically distributed random variables (i.i.d.) from some $\mathcal{N}(0,\sigma^2)$ distribution. With a vector notation, we will write $\mathbf{\varepsilon}\overset{\mathcal{L}}{\sim}\mathcal{N}(\mathbf{0},\sigma^2\mathbb{I})$. The estimated residuals are defined as: $$\widehat{\varepsilon}_i =y_i-[\widehat{\beta}_0+\mathbf{x}_i^T\widehat{\mathbf{\beta}}]$$ Those (estimated) residuals are basic tools for diagnosing the relevance of the model.

An extension of the model described by equation (1) has been proposed to take into account a possible heteroscedastic character: $$(Y\vert \mathbf{X}=\mathbf{x})\overset{\mathcal{L}}{\sim}\mathcal{N}(\mu(\mathbf{x}),\sigma^2(\mathbf{x}))$$where $\sigma^2(\mathbf{x})$ is a positive function of the explanatory variables. This model can be rewritten as: $$y_i=\beta_0+\mathbf{x}_i^T\mathbf{\beta}+\sigma^2(\mathbf{x}_i)\cdot\varepsilon_i$$where residuals are always i.i.d., with unit variance, $$\varepsilon_i=\frac{y_i-[\beta_0+\mathbf{x}_i^T\mathbf{\beta}]}{\sigma(\mathbf{x}_i)}$$ While residuals based equations are popular in linear econometrics (when the dependent variable is continuous), it is no longer popular in counting models, or logistic regression.

However, writing using an error term (as in equation (3)) raises many questions about the representation of an economic relationship between two quantities. For example, it can be assumed that there is a relationship (linear to begin with) between the quantities of a traded good, $q$ and its price $p$. This allows us to imagine a supply equation$$q_i=\beta_0+\beta_1 p_i+u_i$$($u_i$ being an error term) where the quantity sold depends on the price, but in an equally legitimate way, one can imagine that the price depends on the quantity produced (what one could call a demand equation), $$p_i=\alpha_0+\alpha_1 q_i+v_i$$($v_i$ denoting another error term). Historically, the error term in equation (3) could be interpreted as an idiosyncratic error on the variable $y$, the so-called explanatory variables being assumed to be fixed, but this interpretation often makes the link between an economic relationship and a complicated economic model difficult, the economic theory speaking abstractly about a relationship between a magnitude, the econometric model imposing a specific shape (what magnitude is $y$ and what magnitude is $x$) as shown in more detail in Morgan (1990) Chapter 7.

(references mentioned above are online here). To be continued…

[1] This approach can be compared to structural econometrics, as presented for example in Kean (2010).

[2] Here, we will try to distinguish $\beta_0$, the intercept, and the other parameters $\mathbf{\beta}$, since they are considered differently in many extensions (e.g. regularization). Nevertheless, in many expressions $\mathbf{\beta}$ will denote the joint vector $(\beta_0, \mathbf{\beta})$, for general formulas, to avoid too heavy notations.

[3] In the sense that the difference between variance matrices is a positive matrix.

# Simple Distributions for Mixtures?

The idea of GLMs is that given some covariates$X$$Y|X$ has a distribution in the exponential family (Gaussian, Poisson, Gamma, etc). But that does not mean that $Y$ has a similar distribution… so there is no reason to test for a Gamma model for $Y$ before running a Gamma regression, for instance. But are there cases where it might work? That the non-conditional distribution is the same (same family at least) than the conditional ones?

For instance, if $(X,Y)$ has a joint Gaussien distribution, then both marginals are Gaussian, but also $Y|X$. So, in that case, if the covariate is normally distributed, it is possible to have a Gaussian distribution also for $Y$. The econometric interpretation is that with a standard Gaussian linear model, if$X$ is normally distributed, not only the conditional distribution $Y|X$ is Gaussian but also the non-conditional distribution of $Y$.

> set.seed(1)
> n=1e3
> X=rnorm(n,10,2)
> Y=1+3*X+rnorm(n)
> plot(X,Y,xlim=c(4,20))

Indeed, here the distribution of $Y$ is also Gaussian

> library(nortest)

Anderson-Darling normality test

data:  Y
A = 0.23155, p-value = 0.802

> shapiro.test(Y)

Shapiro-Wilk normality test

data:  Y
W = 0.99892, p-value = 0.8293

(not only from a statistical point of view, the thoery of Gaussian random vectors confirms that the non-conditional distribution is Gaussian actually)

Here $X$ is continuous. What if we consider a finite mixture here, i.e.$X$ takes only a finite number of values? Actually, Teicher (1963) proved that it is not possible to have a non-conditional Gaussian distribution for $Y$. But in practice, would we really reject the Gaussian assumption, for $Y$? If the number of classes is to small, yes. But with a large number of classes (a sufficiently large number of mixture components), it is possible,

> pv=function(k=2){
+ n=1e4
+ X=rnorm(n,10,2)
+ Q=quantile(X,(0:k)/k)
+ Q[1]=0
+ Xc=cut(X,Q,labels=1:k)
+ XcN=tapply(X,Xc,mean)
+ Xn=XcN[as.numeric(Xc)]
+ Y=1+3*Xn+rnorm(n)
> b=summary(reg0)$dispersion > lines(x,dgamma(x,shape=1/b,scale=a*b),col="blue") Now, we need a covariate, to run some regressions. What I wanted is some variable slightly correlated with our previous variable. Slightly, just to make sure that our $p$-value in the regression will be close to 5% or 10%. So here, I did generate a variable so that the pair has Clayton copula, with coefficient 0.1 (which is small, extremely small) > a=.1 > set.seed(5) > n=200 > U=runif(n); > V=(U^(-a)*(runif(n)^(-a/(1+a))-1)+1)^(-1/a) > Y=-log(U) > X=qnorm(V) To visualize the copula of the variables, we can use > cop=function(u,v){ + (a+1)*(u*v)^(-(a+1))* + (u^(-a)+v^(-a)-1)^(-(2*a+1)/a) } > x=y=seq(.05,.95,by=.05) > z=outer(x,y,cop) > mat=persp(x,y,z,col="green",shade=TRUE,xlim=c(0,1),ylim=c(0,1),zlim=c(0,2),theta=-30, + ticktype ="detailed",zlab="") We should be not far away from the independence (actually, there is a negative – significant – correlation (Pearson’s correlation)). Now, consider two models, • a Gaussian model (here a standard linear model) • a gamma model, with a linear link function The outputs are the following (you will recognize the outputs given previously) > reg1=lm(Y~X) > reg2=glm(Y~X,family=Gamma(link="identity")) > summary(reg1) Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.92883 0.06391 14.534 <2e-16 *** X -0.12499 0.06108 -2.046 0.0421 * --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.9021 on 198 degrees of freedom Multiple R-squared: 0.02071, Adjusted R-squared: 0.01576 F-statistic: 4.187 on 1 and 198 DF, p-value: 0.04206 > summary(reg2) Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.92901 0.06270 14.817 <2e-16 *** X -0.09883 0.05816 -1.699 0.0909 . --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for Gamma family taken to be 0.9086447) Null deviance: 229.72 on 199 degrees of freedom Residual deviance: 226.58 on 198 degrees of freedom AIC: 379.22 Number of Fisher Scoring iterations: 10 And here are the two predictions, So, which model should we use? As usual, my answer will be “let’s have a look at the data” instead of looking only at tables of figures. Using some code posted a few days ago, let us visualize the two regressions. The Gaussian model is here (for the lower part, I do not go below 0 since we do have, here, a positive variable that we would like to model) while the gamma on is here And if we believe that the explanatory variable has no predictive power (since we can claim that the parameter is not significant in the regression), and we remove it from the regression, we get Here, I do believe that the gamma (not to say the exponential) model is better because it is clearly more coherent with properties of the variable of interest. I trust more the confidence interval obtained above on the gamma model, than the one obtained with a Gaussian distribution. Even if the parameter in the regression is “more significant”. # GLM, non-linearity and heteroscedasticity Last week in the non-life insurance course, we’ve seen the theory of the Generalized Linear Models, emphasizing the two important components • the link function (which is actually the key component in predictive modeling) • the distribution, or the variance function Just to illustrate, consider my favorite dataset ­lin.mod = lm(dist~speed,data=cars) A linear model means here $Y_i=\beta_0+\beta_1X_i+\varepsilon_i$ where the residuals are assumed to be centered, independent, and with identical variance. If we visualize that linear regression, we usually see something like that The idea here (in GLMs) is to assume $Y\vertX=x\sim\mathcal{N}(\beta_0+\beta_1x,\sigma^2)$ which will produce the same model as the one describe previously, based on some error term. That model can be visualized below, attach(cars) n=2 X= cars$speed
Y=cars$dist df=data.frame(X,Y) vX=seq(min(X)-2,max(X)+2,length=n) vY=seq(min(Y)-15,max(Y)+15,length=n) mat=persp(vX,vY,matrix(0,n,n),zlim=c(0,.1),theta=-30,ticktype ="detailed", box = FALSE) reggig=glm(Y~X,data=df,family=gaussian(link="identity")) x=seq(min(X),max(X),length=501) C=trans3d(x,predict(reggig,newdata=data.frame(X=x),type="response"),rep(0,length(x)),mat) lines(C,lwd=2) sdgig=sqrt(summary(reggig)$dispersion)
x=seq(min(X),max(X),length=501)
y1=qnorm(.95,predict(reggig,newdata=data.frame(X=x),type="response"), sdgig)
C=trans3d(x,y1,rep(0,length(x)),mat)
lines(C,lty=2)
y2=qnorm(.05,predict(reggig,newdata=data.frame(X=x),type="response"), sdgig)
C=trans3d(x,y2,rep(0,length(x)),mat)
lines(C,lty=2)
C=trans3d(c(x,rev(x)),c(y1,rev(y2)),rep(0,2*length(x)),mat)
polygon(C,border=NA,col="yellow")
C=trans3d(X,Y,rep(0,length(X)),mat)
points(C,pch=19,col="red")
n=8
vX=seq(min(X),max(X),length=n)
mgig=predict(reggig,newdata=data.frame(X=vX))
sdgig=sqrt(summary(reggig)$dispersion) for(j in n:1){ stp=251 x=rep(vX[j],stp) y=seq(min(min(Y)-15,qnorm(.05,predict(reggig,newdata=data.frame(X=vX[j]),type="response"), sdgig)),max(Y)+15,length=stp) z0=rep(0,stp) z=dnorm(y, mgig[j], sdgig) C=trans3d(c(x,x),c(y,rev(y)),c(z,z0),mat) polygon(C,border=NA,col="light blue",density=40) C=trans3d(x,y,z0,mat) lines(C,lty=2) C=trans3d(x,y,z,mat) lines(C,col="blue")} We do have two parts here: the linear increase of the average, $\mathbb{E}(Y\vert X=x)=\beta_0+\beta_1x$ and the constant variance of the normal distribution $\text{Var}(Y\vert X=x)=\sigma^2$. On the other hand, if we assume a Poisson regression, poisson.reg = glm(dist~speed,data=cars,family=poisson(link="log")) we have something like This time, two things have changed simultaneously: our model is no longer linear, it is an exponential one $\mathbb{E}(Y\vert X=x)=e^{\beta_0+\beta_1x}$, and the variance is also increasing with the explanatory variable $\text{Var}(Y\vert X=x)=e^{\beta_0+\beta_1x}$, since with a Poisson regression, $Y\vert X=x\sim\mathcal{P}(e^{\beta_0+\beta_1x})$ If we adapt the previous code, we get The problem is that we changed two things when we introduced the Poisson regression from the linear model. So let us look at what happens when we change the two components independently. First, we can change the link function, with a Gaussian model but this time a multiplicative model (with a logarithm link function) gaussian.reg = glm(dist~speed,data=cars,family=gaussian(link="log")) which is still, here, an homoscedasctic model, but this time non-linear. Or we can change the link function in the Poisson regression, to get a linear model, but heteroscedastic poisson.lin = glm(dist~speed,data=cars,family=poisson(link="identity")) So this is basically what GLMs are about…. # Bounding sums of random variables, part 2 It is possible to go further, much more actually, on bounding sums of random variables (mentioned in the previous post). For instance, if everything has been defined, in that previous post, on distributions on $\mathbb{R}^+$, it is possible to extend bounds of distributions on $\mathbb{R}$. Especially if we deal with quantiles. Everything we’ve seen remain valid. Consider for instance two $\mathcal{N}(0,1)$ distributions. Using the previous code, it is possible to compute bounds for the quantiles of the sum of two Gaussian variates. And one has to remember that those bounds are sharp. > Finv=function(u) qnorm(u,0,1) > Ginv=function(u) qnorm(u,0,1) > n=1000 > Qinf=Qsup=rep(NA,n-1) > for(i in 1:(n-1)){ + J=0:i + Qinf[i]=max(Finv(J/n)+Ginv((i-J)/n)) + J=(i-1):(n-1) + Qsup[i]=min(Finv((J+1)/n)+Ginv((i-1-J+n)/n)) + } Actually, it is possible to compare here with two simple cases: the independent case, where the sum has a $\mathcal{N}(0,2)$ distribution, and the comonotonic case where the sum has a $\mathcal{N}(0,2^2)$ distribution. > lines(x,qnorm(x,sd=sqrt(2)),col="blue",lty=2) > lines(x,qnorm(x,sd=2),col="blue",lwd=2) On the graph below, the comonotonic case (usually considered as the worst case scenario) is the plain blue line (with here an animation to illustrate the convergence of the numerical algorithm) Below that (strong) blue line, then risks are sub-additive for the Value-at-Risk, i.e. $H^{-1}(u)\leq F^{-1}(u)+G^{-1}(u)$ but above, risks are super-additive for the Value-at-RIsk. i.e. $H^{-1}(u)\geq F^{-1}(u)+G^{-1}(u)$ (since for comonotonic variates, the quantile of the sum is the sum of quantiles). It is possible to visualize those two cases above, in green the area where risks are super-additive, while the yellow area is where risks are sub-additive. Recall that with a Gaussian random vector, with correlation $r$ then the quantile is the quantile of a random variable centered, with variance $2(1+r)$. Thus, on the graph below, we can visualize case that can be obtained with this Gaussian copula. Here the yellow area can be obtained with a Gaussian copulas, the upper and the lower bounds being respectively the comonotonic and the countermononic cases. But the green area can also be obtained when we sum two Gaussian variables ! We just have to go outside the Gaussian world, and consider another copula. Another point is that, in the previous post, $C^-$ was the lower Fréchet-Hoeffding bound on the set of copulas. But all the previous results remain valid if $C^-$ is alower bound on the set of copulas of interest. Especially $\tau_{C^-,L}(F,G)\leq \sigma_{C,L}(F,G)\leq\rho_{C^-,L}(F,G)$ for all $C$ such that $C\geq C^-$. For instance, if we assume that the copula should have positive dependence, i.e. $C\geq C^\perp$, then $\tau_{C^\perp,L}(F,G)\leq \sigma_{C,L}(F,G)\leq\rho_{C^\perp,L}(F,G)$ Which means we should have sharper bounds. Numerically, it is possible to compute those sharper bounds for quantiles. The lower bound becomes $\sup_{u\in[x,1]}\left\{F^{-1}(u)+G^{-1}\left(\frac{x}{u}\right)\right\}$ while the upper bound is $\sup_{u\in[0,x]}\left\{F^{-1}(u)+G^{-1}\left(\frac{x-u}{1-u}\right)\right\}$ Again, one can easily compute those quantities on a grid of the unit interval, > Qinfind=Qsupind=rep(NA,n-1) > for(i in 1:(n-1)){ + J=1:(i) + Qinfind[i]=max(Finv(J/n)+Ginv((i-J)/n/(1-J/n))) + J=(i):(n-1) + Qsupind[i]=min(Finv(J/n)+Ginv(i/J)) + } We get the graph below (the blue area is here to illustrate how sharper those bounds get with the assumption that we do have positive dependence, this area been attained only with copulas exhibiting non-positive dependence) For high quantiles, the upper bound is rather close to the one we had before, since worst case are probably obtained when we do have positive correlation. But it will strongly impact the lower bound. For instance, it becomes now impossible to have a negative quantile, when the probability exceeds 75% if we do have positive dependence… > Qinfind[u==.75] [1] 0 # (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 $(U_i,V_i)$ on the unit square, consider some transformed sample $(Q(U_i),Q(V_i))$, where $Q:(0,1)\rightarrow\mathbb{R}$ is a given function. E.g. a quantile function of an unbounded distribution, for instance the quantile function of the $\mathcal{N}(0,1)$ 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 and tail dependence, part 1

As mentioned in the course last week Venter (2003) suggested nice functions to illustrate tail dependence (see also some slides used in Berlin a few years ago).

• Joe (1990)’s lambda

Joe (1990) suggested a (strong) tail dependence index. For lower tails, for instance, consider

i.e

• Upper and lower strong tail (empirical) dependence functions

The idea is to plot the function above, in order to visualize limiting behavior. Define

for the lower tail, and

for the upper tail, where is the survival copula associated with , in the sense that

while

Now, one can easily derive empirical conterparts of those function, i.e.

and

Thus, for upper tail, on the right, we have the following graph

and for the lower tail, on the left, we have

For the code, consider some real data, like the loss-ALAE dataset.

> library(evd)
> X=lossalae

The idea is to plot, on the left, the lower tail concentration function, and on the right, the upper tail function.

> U=rank(X[,1])/(nrow(X)+1)
> V=rank(X[,2])/(nrow(X)+1)
> Lemp=function(z) sum((U<=z)&(V<=z))/sum(U<=z)
> Remp=function(z) sum((U>=1-z)&(V>=1-z))/sum(U>=1-z)
> u=seq(.001,.5,by=.001)
> L=Vectorize(Lemp)(u)
> R=Vectorize(Remp)(rev(u))
> plot(c(u,u+.5-u[1]),c(L,R),type="l",ylim=0:1,
+ xlab="LOWER TAIL          UPPER TAIL")
> abline(v=.5,col="grey")

Now, we can compare this graph, with what should be obtained for some parametric copulas that have the same Kendall’s tau (e.g.). For instance, if we consider a Gaussian copula,

> tau=cor(lossalae,method="kendall")[1,2]
> library(copula)
> paramgauss=sin(tau*pi/2)
> copgauss=normalCopula(paramgauss)
> Lgaussian=function(z) pCopula(c(z,z),copgauss)/z
> Rgaussian=function(z) (1-2*z+pCopula(c(z,z),copgauss))/(1-z)
> u=seq(.001,.5,by=.001)
> Lgs=Vectorize(Lgaussian)(u)
> Rgs=Vectorize(Rgaussian)(1-rev(u))
> lines(c(u,u+.5-u[1]),c(Lgs,Rgs),col="red")

or Gumbel’s copula,

> paramgumbel=1/(1-tau)
> copgumbel=gumbelCopula(paramgumbel, dim = 2)
> Lgumbel=function(z) pCopula(c(z,z),copgumbel)/z
> Rgumbel=function(z) (1-2*z+pCopula(c(z,z),copgumbel))/(1-z)
> u=seq(.001,.5,by=.001)
> Lgl=Vectorize(Lgumbel)(u)
> Rgl=Vectorize(Rgumbel)(1-rev(u))
> lines(c(u,u+.5-u[1]),c(Lgl,Rgl),col="blue")

That’s nice (isn’t it?), but since we do not have any confidence interval, it is still hard to conclude (even if it looks like Gumbel copula has a much better fit than the Gaussian one). A strategy can be to generate samples from those copulas, and to visualize what we had. With a Gaussian copula, the graph looks like

> u=seq(.0025,.5,by=.0025); nu=length(u)
> nsimul=500
> MGS=matrix(NA,nsimul,2*nu)
> for(s in 1:nsimul){
+ Xs=rCopula(nrow(X),copgauss)
+ Us=rank(Xs[,1])/(nrow(Xs)+1)
+ Vs=rank(Xs[,2])/(nrow(Xs)+1)
+ Lemp=function(z) sum((Us<=z)&(Vs<=z))/sum(Us<=z)
+ Remp=function(z) sum((Us>=1-z)&(Vs>=1-z))/sum(Us>=1-z)
+ MGS[s,1:nu]=Vectorize(Lemp)(u)
+ MGS[s,(nu+1):(2*nu)]=Vectorize(Remp)(rev(u))
+ lines(c(u,u+.5-u[1]),MGS[s,],col="red")
+ }

(including – pointwise – 90% confidence bands)

> Q95=function(x) quantile(x,.95)
> V95=apply(MGS,2,Q95)
> lines(c(u,u+.5-u[1]),V95,col="red",lwd=2)
> Q05=function(x) quantile(x,.05)
> V05=apply(MGS,2,Q05)
> lines(c(u,u+.5-u[1]),V05,col="red",lwd=2)

while it is

with Gumbel copula. Isn’t it a nice (graphical) tool ?

But as mentioned in the course, the statistical convergence can be slow. Extremely slow. So assessing if the underlying copula has tail dependence, or not, it now that simple. Especially if the copula exhibits tail independence. Like the Gaussian copula. Consider a sample of size 1,000. This is what we obtain if we generate random scenarios,

or we look at the left tail (with a log-scale)

Now, consider a 10,000 sample,

or with a log-scale

We can even consider a 100,000 sample,

or with a log-scale

On those graphs, it is rather difficult to conclude if the limit is 0, or some strictly positive value (again, it is a classical statistical problem when the value of interest is at the border of the support of the parameter). So, a natural idea is to consider a weaker tail dependence index. Unless you have something like 100,000 observations…

# 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“.

# Border bias and weighted kernels

With Ewen (aka @3wen), not only we have been playing on Twitter this month, we have also been working on kernel estimation for densities of spatial processes. Actually, it is only a part of what he was working on, but that part on kernel estimation has been the opportunity to write a short paper, that can now be downloaded on hal.

The problem with kernels is that kernel density estimators suffer a strong bias on borders. And with geographic data, it is not uncommon to have observations very close to the border (frontier, or ocean). With standard kernels, some weight is allocated outside the area: the density does not sum to one. And we should not look for a global correction, but for a local one. So we should use weighted kernel estimators (see on hal for more details). The problem that weights can be difficult to derive, when the shape of the support is a strange polygon. The idea is to use a property of product Gaussian kernels (with identical bandwidth) i.e. with the interpretation of having noisy observation, we can use the property of circular isodensity curve. And this can be related to Ripley (1977) circumferential correction. And the good point is that, with R, it is extremely simple to get the area of the intersection of two polygons. But we need to upload some R packages first,

require(maps)
require(sp)
require(snow)
require(ellipse)
require(ks)
require(gpclib)
require(rgeos)
require(fields)

To be more clear, let us illustrate that technique on a nice example. For instance, consider some bodiliy injury car accidents in France, in 2008 (that I cannot upload but I can upload a random sample),

base_cara=read.table(
"http://freakonometrics.blog.free.fr/public/base_fin_morb.txt",
sep=";",header=TRUE)

The border of the support of our distribution of car accidents will be the contour of the Finistère departement, that can be found in standard packages

geoloc=read.csv(
"http://freakonometrics.free.fr/popfr19752010.csv",
colClasses=c(rep("character",5),rep("numeric",38)))
geoloc=geoloc[,c("dep","com","com_nom",
"long","lat","pop_2008")]
geoloc$id=paste(sprintf("%02s",geoloc$dep),
sprintf("%03s",geoloc$com),sep="") geoloc=geoloc[,c("com_nom","long","lat","pop_2008")] head(geoloc) france=map('france',namesonly=TRUE, plot=FALSE) francemap=map('france', fill=TRUE, col="transparent", plot=FALSE) detpartement_bzh=france[which(france%in% c("Finistere","Morbihan","Ille-et-Vilaine", "Cotes-Darmor"))] bretagne=map('france',regions=detpartement_bzh, fill=TRUE, col="transparent", plot=FALSE,exact=TRUE) finistere=cbind(bretagne$x[321:678],bretagne$y[321:678]) FINISTERE=map('france',regions="Finistere", fill=TRUE, col="transparent", plot=FALSE,exact=TRUE) monFINISTERE=cbind(FINISTERE$x[c(8:414)],FINISTERE$y[c(8:414)]) Now we need simple functions, cercle=function(n=200,centre=c(0,0),rayon) {theta=seq(0,2*pi,length=100) m=cbind(cos(theta),sin(theta))*rayon m[,1]=m[,1]+centre[1] m[,2]=m[,2]+centre[2] names(m)=c("x","y") return(m)} poids=function(x,h,POL) {leCercle=cercle(centre=x,rayon=5/pi*h) POLcercle=as(leCercle, "gpc.poly") return(area.poly(intersect(POL,POLcercle))/ area.poly(POLcercle))} lissage = function(U,polygone,optimal=TRUE,h=.1) {n=nrow(U) IND=which(is.na(U[,1])==FALSE) U=U[IND,] if(optimal==TRUE) {H=Hpi(U,binned=FALSE); H=matrix(c(sqrt(H[1,1]*H[2,2]),0,0, sqrt(H[1,1]*H[2,2])),2,2)} if(optimal==FALSE){H= matrix(c(h,0,0,h),2,2) before defining our weights. poidsU=function(i,U,h,POL) {x=U[i,] poids(x,h,POL)} OMEGA=parLapply(cl,1:n,poidsU,U=U,h=sqrt(H[1,1]), POL=as(polygone, "gpc.poly")) OMEGA=do.call("c",OMEGA) stopCluster(cl) }else {OMEGA=lapply(1:n,poidsU,U=U,h=sqrt(H[1,1]), POL=as(polygone, "gpc.poly")) OMEGA=do.call("c",OMEGA)} Note that it is possible to parallelize if there are a lot of observations, if(n>=500) {cl <- makeCluster(4,type="SOCK") worker.init <- function(packages) {for(p in packages){library(p, character.only=T)} NULL} clusterCall(cl, worker.init, c("gpclib","sp")) clusterExport(cl,c("cercle","poids")) Then, we can use standard bivariate kernel smoothing functions, but with the weights we just calculated, using a simple technique that can be related to one suggested in Ripley (1977), fhat=kde(U,H,w=1/OMEGA,xmin=c(min(polygone[,1]), min(polygone[,2])),xmax=c(max(polygone[,1]), max(polygone[,2]))) fhat$estimate=fhat$estimate*sum(1/OMEGA)/n vx=unlist(fhat$eval.points[1])
vy=unlist(fhat$eval.points[2]) VX = cbind(rep(vx,each=length(vy))) VY = cbind(rep(vy,length(vx))) VXY=cbind(VX,VY) Ind=matrix(point.in.polygon(VX,VY, polygone[,1], polygone[,2]),length(vy),length(vx)) f0=fhat f0$estimate[t(Ind)==0]=NA
return(list(
X=fhat$eval.points[[1]], Y=fhat$eval.points[[2]],
Z=fhat$estimate, ZNA=f0$estimate,
H=fhat$H, W=fhat$W))}
lissage_without_c = function(U,polygone,optimal=TRUE,h=.1)
{n=nrow(U)
IND=which(is.na(U[,1])==FALSE)
U=U[IND,]
if(optimal==TRUE) {H=Hpi(U,binned=FALSE);
H=matrix(c(sqrt(H[1,1]*H[2,2]),0,0,sqrt(H[1,1]*H[2,2])),2,2)}
if(optimal==FALSE){H= matrix(c(h,0,0,h),2,2)}
fhat=kde(U,H,xmin=c(min(polygone[,1]),
min(polygone[,2])),xmax=c(max(polygone[,1]),
max(polygone[,2])))
vx=unlist(fhat$eval.points[1]) vy=unlist(fhat$eval.points[2])
VX = cbind(rep(vx,each=length(vy)))
VY = cbind(rep(vy,length(vx)))
VXY=cbind(VX,VY)
Ind=matrix(point.in.polygon(VX,VY, polygone[,1],
polygone[,2]),length(vy),length(vx))
f0=fhat
f0$estimate[t(Ind)==0]=NA return(list( X=fhat$eval.points[[1]],
Y=fhat$eval.points[[2]], Z=fhat$estimate,
ZNA=f0$estimate, H=fhat$H,
W=fhat$W))} So, now we can play with those functions, base_cara_FINISTERE=base_cara[which(point.in.polygon( base_cara$long,base_cara$lat,monFINISTERE[,1], monFINISTERE[,2])==1),] coord=cbind(as.numeric(base_cara_FINISTERE$long),
as.numeric(base_cara_FINISTERE$lat)) nrow(coord) map(francemap) lissage_FIN_withoutc=lissage_without_c(coord, monFINISTERE,optimal=TRUE) lissage_FIN=lissage(coord,monFINISTERE, optimal=TRUE) lesBreaks_sans_pop=range(c( range(lissage_FIN_withoutc$Z),
range(lissage_FIN$Z))) lesBreaks_sans_pop=seq(min(lesBreaks_sans_pop)*.95, max(lesBreaks_sans_pop)*1.05,length=21) plot_article=function(lissage,breaks, polygone,coord){ par(mar=c(3,1,3,1)) image.plot(lissage$X,lissage$Y,(lissage$ZNA),
xlim=range(polygone[,1]),ylim=range(polygone[,2]),
breaks=breaks, col=rev(heat.colors(20)),xlab="",
ylab="",xaxt="n",yaxt="n",bty="n",zlim=range(breaks),
horizontal=TRUE)
contour(lissage$X,lissage$Y,lissage$ZNA,add=TRUE, col="grey") points(coord[,1],coord[,2],pch=19,cex=.1, col="dodger blue") polygon(polygone,lwd=2,)} plot_article(lissage_FIN_withoutc,breaks= lesBreaks_sans_pop,polygone=monFINISTERE, coord=coord) plot_article(lissage_FIN,breaks= lesBreaks_sans_pop,polygone=monFINISTERE, coord=coord) If we look at the graphs, we have the following densities of car accident, with a standard kernel on the left, and our proposal on the right (with local weight adjustment when the estimation is done next to the border of the region of interest), Similarly, in Morbihan, With those modified kernels, hot spots appear much more clearly. For more details, the paper is online on hal. # PhD defense on copulas This Wednesday I will be at Université Paris 1 Sorbonne as a member of the jury of the PhD thesis of Pierre-André Maugis, on conditional correlation and vine copula. Vine copulas were born in 2002 with thepaper of Tim Bedford and Roger M. CookeVines–a new graphical model for dependent random variables. The idea is to use the following decomposition for a multivariate density (from Bayes formula, with synthetic notations). Then using the relationship between a bivariate density and its copula (density) thus Using again Bayes formula, and we can write Since and , the previous expression becomes or to stress on the most important part (as I see it) It is common then to assume that this conditional copula does not depend on the conditioning parameter. The more detailed expression of that joint trivariate density is The (parametric) inference algorithm is defined in Cooke, Joe and Aas (2010) as follows The important assumption in vine copula models is that conditional copulas are constant. And this assumption might be relevant in some cases. For instance, in the Gaussian case (the observations have a Gaussian joint distribution – or at least copula – and we fit a vine model with Gaussian bivariate copulas). The code to fit a vine copula is the following, > library(CDVine) > library(mnormt) > SIGMA=matrix(c(1,.6,.7,.6,1,.8,.7,.8,1),3,3) > X=rmnorm(n=100000,varcov=SIGMA) > CDVineSeqEst(dat=X, family = c(1,1,1), + type = 1, method = "mle")$par
[1] 0.6001505 0.7023699 0.6698215

$par2 [1] 0 0 0 Note that it is consistent with the following algorithm where conditional copulas are fitted. In the following, for all values of the given component, we wit a Gaussian copula for the conditional remaining pair, > U=pnorm(X) > U1U2=U[,1:2] > U1U3=U[,c(1,3)] > GaussCop = normalCopula(param=.5, dim = 2) > U1U2=U[,1:2] > U1U3=U[,c(1,3)] > fit12.mpl = fitCopula(GaussCop, U1U2, method="mpl")@estimate > fit13.mpl = fitCopula(GaussCop, U1U3, method="mpl")@estimate > fit12.mpl [1] 0.5984932 > fit13.mpl [1] 0.7005185 > fit23a=fit23b=rep(NA,99) > for(i in 4:96){ + x=i/100 + C12=pcopula(normalCopula(param=fit12.mpl, dim = 2),U1U2) + C13=pcopula(normalCopula(param=fit13.mpl, dim = 2),U1U3) + U12=rank(C12)/(nrow(U)+1) + U13=rank(C13)/(nrow(U)+1) + U23=cbind(U12[abs(U[,1]-x)<.02],U13[abs(U[,1]-x)<.02]) + V23=cbind(rank(U23[,1])/(nrow(U23)+1), + rank(U23[,2])/(nrow(U23)+1)) + fit23.mpl = fitCopula(GaussCop, V23, method="mpl")@estimate + fit23a[i]=fit23.mpl + } > plot(X,fit23a,col="red") It looks like assuming the conditional copula as constant was a valid assumption here But note that if the true distribution is not Gaussian, then assuming the conditional copula as constant is not valid anymore (here a trivariate Clayton copula was generated) # 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) pv[s]=ad.test(X)$p.value
}
mean(pv>.05)
[1] 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)
[1] -0.605381

mean(IC3)
[1] -0.5759391
quantile(m,.05)
-0.623578