# Convergence and Asymptotic Results

Last week, in our mathematical statistics course, we’ve seen the law of large numbers (that was proven in the probability course), claiming that

$\overline{X}_n\ \xrightarrow{\text{a.s.}}\ \mathbb{E}(X)$

given a collection $\{X_1,\cdots,X_n\}$ of i.i.d. random variables, with

$\overline{X}_n=\frac1n(X_1+\cdots+X_n)$

To visualize that convergence, we can use

> m=100
> mean_samples=function(n=10){
+   X=matrix(rnorm(n*m),nrow=m,ncol=n)
+   return(apply(X,1,mean))
+ }
> B=matrix(NA,100,20)
> for(i in 1:20){
+   B[,i]=mean_samples(i*10)
+ }
> colnames(B)=as.character(seq(10,200,by=10))
> boxplot(B)

It is possible to visualize also the $\sqrt{n}$ bounds (used in the central limit theorem to get a limiting non degenerated distribution)

> u=seq(0,21,by=.2)
> v=sqrt(u*10)
> lines(u,1.96/v,col="red")
> lines(u,-1.96/v,col="red")

Yesterday, we’ve been discussing properties of the empirical cumulative distribution function,

$\hat F_n(x)=\frac{1}{n}\sum_{i=1}^n \boldsymbol{1}(X_i\in(-\infty,x])$

We’ve seen Glivenko-Cantelli theorem, which states that (under mild assumptions)

$\|\hat F_n-F\|_\infty \equiv \sup_{t\in\mathbb{R}} \big|\hat F_n(t)-F(t)\big|\ \xrightarrow{a.s.}\ 0.$

To visualize that convergence use the following code. Here I use the trick

$\max\{a,b\}=\frac{a+b}{2}+\frac{\vert b-a\vert}{2}$

to get the maximum (componentwise) between two matrices

> m=100
> inf_sample=function(n=10){
+ X=matrix(rnorm(n*m),nrow=m,ncol=n)
+ Xs=t(apply(X,1,sort))
+ Pe_inf=matrix(rep((0:(n-1))/n,
+ each=m),nrow=m,ncol=n)
+ Pe_sup=matrix(rep((0:n)/n,each=m),
+ nrow=m,ncol=n)
+ Pt=pnorm(Xs)
+ D1=abs(Pe_inf-Pt)
+ D2=abs(Pe_sup-Pt)
+ Df=(D1+D2)/2+abs(D2-D1)/2
+ return(apply(Df,1,max))
+ }
> B=matrix(NA,100,20)
> for(i in 1:20){
+   B[,i]=inf_sample(i*10)
+ }
> colnames(B)=as.character(seq(10,200,by=10))
> boxplot(B)

We have also discussed the pointwise asymptotic normality of the empirical cumulative distribution function

$\sqrt{n}\big(\hat F_n(t) - F(t)\big)\ \ \xrightarrow{\mathcal{L}}\ \ \mathcal{N}\Big( 0, F(t)\big(1-F(t)\big) \Big).$

Here again, it is possible to visualize it. The first step is to compute several trajectories for empirical cumulative distribution function

> u=seq(-3,3,by=.1)
> plot(u,u,ylim=c(0,1),col="white")
> M=matrix(NA,length(u),1000)
> for(m in 1:1000){
+ n=100
+ x=rnorm(n)
+ Femp=Vectorize(function(t) mean(x<=t))
+ v=Femp(u)
+ M[,m]=v
+ lines(u,v,col='light blue',type="s")
+ }

Note that we can compute (pointwise) confidence bands

> lines(u,apply(M,1,mean),col="red",type="l")
> lines(u,apply(M,1,function(x) quantile(x,.05)),
+ col="red",type="s")
> lines(u,apply(M,1,function(x) quantile(x,.95)),
+ col="red",type="s")

Now, if we focus on one specific point, we can visualize the asmptotic normality (i.e. the almost normality when we have a sample of size 100)

> x0=-1
> y=M[which(u==x0),]
> hist(y,probability=TRUE,
+ breaks=seq(.015,0.55,by=.01))
> vu=seq(0,1,by=.001)
> lines(vu,dnorm(vu,pnorm(x0),
+ sqrt((pnorm(x0)*(1-pnorm(x0)))/100)),
+ col="red")

# Defining Properly MA(∞) Time Series

In order to properly define $MA(\infty)$ series, we need to get back on some properties of infinite sequences, as briefly mentioned yesterday in the MAT8181 course. Consider some sequence $(a_i)_{i\in\mathbb{N}}$. The sequence is said to be summable if

$S_n=\sum_{i=0}^n a_i$

is convergent, i.e. if the limit of $S_n$ exists when $n\rightarrow\infty$.

From Cauchy criterion$\sum a_i$ converges if and only if for each $\eta>0$, there is $m\in\mathbb{N}$ for which

$\vert a_i+a_{i+1}+\cdots+a_{j-1}+a_j\vert<\eta$

when $i,j>m$. The sequence $(a_i)_{i\in\mathbb{N}}$ is said to be absolutely summable if

$\sum_{i=0}^\infty \vert a_i\vert <\infty$

and square-summable if

$\sum_{i=0}^\infty a_i^2 <\infty$

Observe that absolute summability will imply square summability (since for $j$‘s large enough $\vert a_j\vert \leq1$, and then $a_j^2\leq\vert a_j\vert$)

Consider now some $MA(\infty)$ time series

$X_t=\sum_{h=0}^\infty \theta_h \varepsilon_{t-h}$

If the sequence of coefficients $(\theta_i)$ is square-summable, then

$S_T = \sum_{h=0}^T \theta_h \varepsilon_{t-h}$

converges in $L_2$  to some random varible as $T\rightarrow\infty$. This can be proved easily using Cauchy criteria, in the sense that for any $\eta>0$, there is a $T$ large enough such that, for any $h$,

$\underbrace{\mathbb{E}\left(\sum_{i=T}^h\theta_i \varepsilon_{t-i}\right)^2}_{\sigma^2 \cdot \sum_{i=T}^h \theta_i^2} \leq \eta$

In that case, if the sequence of coefficients $(\theta_i)$ is square-summable, then $(X_t)$ is stationary (in the $L_2$ sense) since the process is centered, and

$\gamma(h)=\sigma^2 \cdot \sum_{i=0}^\infty \theta_i \theta_{i+h}$

for all $h\in\mathbb{N}$.

Further, ergodicity of the time series, define as the absolute summability of the autocovariance sequence, is obtained when the sequence of coefficients $(\theta_i)$ is absolutely summable.

# Simple and heuristic optimization

This week, at the Rmetrics conference, there has been an interesting discussion about heuristic optimization. The starting point was simple: in complex optimization problems (here we mean with a lot of local maxima, for instance), we do not necessarily need extremely advanced algorithms that do converge extremly fast, if we cannot ensure that they reach the optimum. Converging extremely fast, with a great numerical precision to some point (that is not the point we’re looking for) is useless. And some algorithms might be much slower, but at least, it is much more likely to converge to the optimum. Wherever we start from.
We have experienced that with Mathieu, while we were looking for maximum likelihood of our MINAR process: genetic algorithm have performed extremely well. The idea is extremly simple, and natural. Let us consider as a starting point the following algorithm,

1. Start from some
2. At step , draw a point  in a neighborhood of
• either  then
• or  then

This is simple (if you do not enter into details about what such a neighborhood should be). But using that kind of algorithm, you might get trapped and attracted to some local optima if the neighborhood is not large enough. An alternative to this technique is the following: it might be interesting to change a bit more, and instead of changing when we have a maximum, we change if we have almost a maximum. Namely at step ,

• either then
• or  then

for some . To illustrate the idea, consider the following function

> f=function(x,y) { r <- sqrt(x^2+y^2);
+ 1.1^(x+y)*10 * sin(r)/r }
(on some bounded support). Here, by picking noise and  values arbitrary, we have obtained the following scenarios
> x0=15
> MX=matrix(NA,501,2)
> MX[1,]=runif(2,-x0,x0)
> k=.5
> for(s in 2:501){
+  bruit=rnorm(2)
+  X=MX[s-1,]+bruit*3
+  if(X[1]>x0){X[1]=x0}
+  if(X[1]<(-x0)){X[1]=-x0}
+  if(X[2]>x0){X[2]=x0}
+  if(X[2]<(-x0)){X[2]=-x0}
+  if(f(X[1],X[2])+k>f(MX[s-1,1],
+    MX[s-1,2])){MX[s,]=X}
+  if(f(X[1],X[2])+k<=f(MX[s-1,1],
+    MX[s-1,2])){MX[s,]=MX[s-1,]}
+}

It does not always converge towards the optimum,

and sometimes, we just missed it after being extremely unlucky

Note that if we run 10,000 scenarios (with different random noises and starting point), in 50% scenarios, we reach the maxima. Or at least, we are next to it, on top.

What if we compare with a standard optimization routine, like Nelder-Mead, or quasi gradient ?Since we look for the maxima on a restricted domain, we can use the following function,

> g=function(x) f(x[1],x[2])
> optim(X0, g,method="L-BFGS-B",
[1] 42.98
> X=rep(1,nrow(cars)); Y=cars$dist > (solve(t(X)%*%X))%*%(t(X)%*%Y) [,1] [1,] 42.98 • sur les propriétés des estimateurs Une question était “doit-on vérifier les propriétés de l’estimateur (sans biais, efficace etc…) du summary de la Reg ? pour chaque régression ? Bon, la réponse est non car c’est impossible ! Pour répéter tranquillement ce que j’avais dit plusieurs fois oralement: quand on fait un modèle (théorique), on part d’hypothèses (par exemple ici les résidus sont centrées, de variance constante). Sous ces hypothèses on a des propriétés: ce sont des théorèmes, autrement dit de la théorie. Le plus connu étant le théorème de Gauss-Markov sur le modèle linéaire, qui garantie que l’estimateur par moindre carré est BLUE (best linear unbiased estimator). Autrement dit (si je retraduit ce que dit ce théorème), si les hypothèses sont valides, alors la théorie nous garantie que les estimateurs vérifient des propriétés, dont celle d’être sans biais, par exemple. On ne peut pas vérifier qu’un estimateur est sans biais, on peut juste vérifier ex post que les hypothèses sont valides (ou non), ce qui garantit (ou pas) l’absence de biais. Mais comme on va le voir par la suite, le biais ce n’est pas forcément gênant, ce qui est important, c’est surtout la convergence. Comme le notait Clive Granger “if you can’t even get a consistant estimator, you shouldn’t be in this business“. • et sur des points plus théoriques sinon la personne qui avait des soucis numériques à calculer les coefficients me demandait si cela avait un “rapport avec les plim X’E/n et plim (X’X)-1/ n“. Pour revenir sur ce point, une hypothèse forte dans le modèle linéaire (de base) est que le bruit soit soit un vrai bruit, c’est à dire non corrélé avec la variable explicative, i.e. En effet, sinon autrement dit l’estimateur par moindre carré est biasé. La solution la plus naturelle pour éviter ce “problème” est d’utiliser un instrument (ou une variable instrumentale). Mais surtout, sur cet exemple, on voit qu’on ne pourra avoir la convergence de notre estimateur par moindres carrés que si la corrélation entre le bruit et la variable explicative tend vers 0 asymptotiquement. Il existe en effet une différence fondamentale entre le biais et la consistance (ou la convergence vers la vraie valeur). • un estimateur sans biais est généralement convergent (et de manière consistante) • un estimateur biaisé peut ne pas converger Prenons un petit exemple pour illustrer ces points. On supposera n1 < n2 < n3par la suite. Comme je l’ai expliqué auparavant, la distribution de l’estimateur est une propriété théorique qui ne peut se voir sur un jeu de données (à moins de faire des simulations mais ça sort du cadre de mon billet d’aujourd’hui). La figure ci-dessous correspond au cas sans biais, et convergent La figure ci-dessous correspond au cas biaisé, mais convergent, et consistant, au sens où asymptotiquement, l’estimateur sera sans biais. Enfin, un dernier cas correspondant au cas biaisé, convergent, mais non consistant. Asymptotiquement l’estimateur sera toujours biaisé, Mais revenons un peu sur la formaliation sous-jacente au modèle linéaire. Les propriétés sur les résidus sont conditionnelles aux variables explicatives, en particulier. (c’est l’hypothèse d’homoscédasticité et d’indépendance des résidus). Bon, il faut aussi l’absence de relation linéaire entre les variables explicatives, ce qui se traduit parfois comme une absence de collinéarité. On peut s’intéresser aux propriétés asymptotiques de l’estimateur. Pour cela, il faut peut être rappeler ce que signifie la convergence pour une variable aléatoire (en statistique mathématique, les estimateurs sont vus comme des variables aléatoires). On dira que converge en probabilité vers , parfois noté si Par exemple la loi des grands nombres garantie que, pour un échantillon i.i.d. si l’espérance des variables est finie, alors Alors sous les hypothèses mentionnées auparavant, on peut montrer que, dans le modèle linéaire, On parlera alors d’estimateur convergent. On peut aussi montrer que Mais formellement, si Voilà en gros quelques éléments pour répondre aux questions, mais tout ça sera repris en détails dans le cours d’économétrie 1 qui commencera très bientôt… et surtout d’économétrie 2 qui abordera les aspects dynamiques. # Delta method avec R, ou transformation des paramètres Un petit billet aujourd’hui sur une question qui m’a été posée (et qui m’a permis de découvrir des nouvelles fonctions sous R). On suppose que l’on a un modèle économétrique de la forme que l’on peut encore écrire Supposons que ce soit la formulation en terme de ces derniers coefficients qui nous intéresse, et que l’on cherche un intervalle de confiance. La Delta Method nous dit que si on dispose d’un estimateur asymptotiquement gaussien (mais c’est le cas pour la plupart des estimateurs que l’on construit en pratique), tout transformation bijective de cet estimateur reste asymptotiquement Gaussien. Bref, on va pouvoir utiliser ces résultats pour en déduire un intervalle de confiance (asymptotique) pour le nouveau paramètre. En effet, l’estimateur par moindre carrés vérifie la normalité asymptotique (sous quelques hypothèses) ce qui donne avec les notations précédentes Autrement dit, on peut facilement en déduire un intervalle de confiance (asymptotique). Pour cela, on a trois méthodes (au moins) avec R. Pour la suite, je vais considerer la régression sur la base de données cars (comme toujours), et on s’intéresse à l’estimation de (ne me demandez pas pourquoi, il me fallait un exemple numérique, j’ai pris un truc nonlinéaire pour qu’au moins ça puisse avoir un vague intérêt). • En faisant un peu de calcul formel Pour les calculs symboliques, il y a la fonction D dans le package de base, > f <- expression(sin(cos(x + y^2))) > ( Dfx <- D(f, "x") ) -(cos(cos(x + y^2)) * sin(x + y^2)) ce qui pourra alors être utilisé pour faire les calculs de la variance asymptotique. > x=c(0,1,2); y=1 > eval(D.sc) [1] -0.7216061 -0.8316919 -0.0774320 On peut également utiliser cette fonction pour faire plein de jolis dessins, > x <-seq(-5, 5, length=100) > fx=x^2 > Df <-D(expression(x^2), "x") > Dfx<-eval(Df) > plot(x, fx, type = "l", las = 1,ylim=c(-2,2)) > for (i in 1:length(x)) { + abline(a = fx[i] - Dfx[i] * x[i], b = Dfx[i], + col = rainbow(100)[i]) + } Mais bon, c’était juste pour rigoler un peu… On peut utiliser cette fonction pour calculer l’écart type de l’estimateur de la grandeur qui nous intéresse. > reg1 = lm(dist~speed, data = cars) > estmean = coef(reg1) > estvar = summary(reg1)$cov.unscaled * summary(reg1)$sigma^2 > f = function(x,y){ + -y/(x+y) + } > X=f(estmean[1],estmean[2]) > g=expression(-y/(x+y)) > x<-as.numeric(estmean[1]); y<-as.numeric(estmean[2]) > V=eval(c(eval(D(g, "x")),eval(D(g, "y")))) > sqrt(t(V) %*% estvar %*% V) [,1] [1,] 0.1063255 Voilà pour la première méthode. Petite remarque: pour faire le calcul symbolique de la dérivée, j’utilise expression et non pas function. • En utilisant la fonction deltamethod de la librarie msm Ici on commence comme avant, sauf que la phase de calcul est immédiat… > reg2 <- lm(dist~speed, data = cars) > estmean = coef(reg2) > estvar = summary(reg2)$cov.unscaled * summary(reg2)$sigma^2 > deltamethod (~ -x2 / (x1 + x2), estmean, estvar) [1] 0.1063255 et d’ailleurs, ça peut s’appliquer pour un paramètre en dimension 2 > deltamethod (list(~ -x2 / (x1 + x2), ~ x1+x2), estmean, estvar) [1] 0.1063255 6.3664368 En revanche, on obtient des estimateur des écart-types, et pas de la matrice de variance-covariance des estimateurs. Autrement dit on ne peut pas en déduire l’ellipse de confiance, seulement des intervalles de confiance marginaux (en utilisant la normalité asymptotique). • En utilisant la fonction delta.method de librarie alr3 Pour finir, une dernière fonction. La seule chose qui change c’est la synthaxe. Mais la sortie est plus élégante, > reg3 = lm(dist~speed, data = cars) > delta.method(reg3, "-b1/(b1+b0)") Functions of parameters: expression(-b1/(b1+b0)) Estimate = 0.2881585 with se = 0.1063255 En fait, l’écart-type indiqué ici est précisément celui qui est utilisé pour écrire l’intervalle de confiance. Dans cet exemple, l’intervalle de confiance approché (car on suppose que la normalité – en théorie seulement asymptotique – est vérifiée) et estimé (car les paramètres sont inconnus et estimés) est de la forme suivant, pour un seuil de confiance à 95% > DM=delta.method(reg3, "-b1/(b1+b0)") Functions of parameters: expression(-b1/(b1+b0)) Estimate = 0.2881585 with se = 0.1063255 > DM$estimate+qt(.975,df=nrow(cars)-1)*DM$se [1] 0.5018277 > DM$estimate+qt(.025,df=nrow(cars)-1)*DM\$se
[1] 0.07448936

C’est à dire qu’on a un intervalle de confiance de la forme [0.074;0.502].