Somewhere else, part 112

Some writings worth reading, found somewhere, on the internet

et un peu de lecture en français

Did I miss something interesting?

Inference for MA(q) Time Series

Yesterday, we’ve seen how inference for https://latex.codecogs.com/gif.latex?AR(p) time series was possible.  I started  with that one because it is actually the simple case. For instance, we can use ordinary least squares. There might be some possible bias (see e.g. White (1961)), but asymptotically, estimators are fine (consistent, with asymptotic normality). But when the noise is (auto)correlated, then it is more complex. So, consider here some  time series

for some white noise https://latex.codecogs.com/gif.latex?(\varepsilon_t).

> theta1=.25
> theta2=.7
> n=1000
> set.seed(1)
> e=rnorm(n)
> Z=rep(0,n)
> for(t in 3:n) Z[t]=e[t]+theta1*e[t-1]+theta2*e[t-2]
> Z=Z[800:1000]
> plot(Z,type="l")

  • Using the empirical autocorrelations

The first idea might be to use the first two (empirical) autocorrelations (the two that are supposed to be – theoretically – non null).

with  when . We also have the following relationship on the variance of the process

With those three equations, for three unknown parameters, https://latex.codecogs.com/gif.latex?\theta_1https://latex.codecogs.com/gif.latex?\theta_2 and https://latex.codecogs.com/gif.latex?\sigma, we simply have to solve (numerically) that system of equations,

> v=c(as.numeric(acf(Z)$acf[2:3]),var(Z))
> v
[1] 0.1658760 0.3823053 1.6379498
> library(rootSolve)
> seteq=function(x){
+ F1=v[1]-(x[1]+x[1]*x[2])/(1+x[1]^2+x[2]^2)
+ F2=v[2]-(x[2])/(1+x[1]^2+x[2]^2)
+ F3=v[3]-(1+x[1]^2+x[2]^2)*x[3]^2
+ return(c(F1,F2,F3))}
> multiroot(f=seteq,start=c(.1,.1,1))
$root
[1] 0.1400579 0.4766699 1.1461636

$f.root
[1]  7.876355e-10  4.188458e-09 -2.839977e-09

$iter
[1] 5

$estim.precis
[1] 2.605357e-09

We are a bit far away from the true values, used to generate our sample. And if we consider 1,000 sample (instead of only one), we still have the bias, and a large variance for our three estimators,

http://freakonometrics.hypotheses.org/files/2014/01/Capture-d%E2%80%99e%CC%81cran-2014-01-29-a%CC%80-11.34.46.png

  • Using least square techniques

We can try something quite different here. The problem we have is that we do not observe the noise https://latex.codecogs.com/gif.latex?(\varepsilon_t), we only observe our series https://latex.codecogs.com/gif.latex?(X_t). But we can try to rebuild that series (call it  since we’re not sure it will be a reconstruction of the noise). As suggested in Box & Jenkins (1967), assume that the first two values are null. And then, use

and then, we can use least square techniques

The code will be

> V=function(p){
+ theta1=p[1]
+ theta2=p[2]
+ u=rep(0,length(Z))
+ for(t in 3:length(Z)) u[t]=Z[t]-theta1*u[t-1]-theta2*u[t-2]
+ return(sum(u^2))
+ }

If we try to minimize the sum of the squares of the residuals, we get

> optim(par=c(.1,.1),V)
$par
[1] 0.2751667 0.6723909

$value
[1] 225.8104

$counts
function gradient 
      77       NA 

$convergence
[1] 0

$message
NULL

which is close to the true value. Another good thing is that, if we compare that rebuilt noise with the true one (since we actually have it), then we have the same vector,

> plot(e[800:1000],col="blue",type="l")
> theta1=0.2751667
> theta2=0.6723909
> u=rep(0,length(Z))
> for(t in 3:length(Z)) u[t]=Z[t]-theta1*u[t-1]-theta2*u[t-2]
> lines(1:201,u,col="red")

So far, so good. And if we look at 1,000 samples, we get

It looks like we have some bias here. And since the two estimators should be negatively correlated, one over-estimates, while the other one under-estimates.

  • Using the (global) maximum likelihood technique

And a final method might be to use the maximum likelihood technique (globally). Again, if we assume that we have a Gaussian i.i.d noise, then the vector https://latex.codecogs.com/gif.latex?\boldsymbol{Y}=(Y_1,\cdots,Y_t) is Gaussian, with a simple variance matrix (since a lot of elements will be null),

> library(mnormt)
> GlobalLogLik=function(A,TS){
+ n=length(TS)
+ theta1=A[1];  theta2=A[2]
+ sigma=A[3]
+ SIG=matrix(0,n,n)
+ rho=rep(0,n)
+ rho[1]=1
+ rho[2]=(theta1+theta1*theta2)/(1+theta1^2+theta2^2)
+ rho[3]=(theta2)/(1+theta1^2+theta2^2)
+ for(i in 1:n){for(j in 1:n){
+ SIG[i,j]=rho[abs(i-j)+1]}}
+ gamma0=(1+theta1^2+theta2^2)*sigma^2
+ SIG=gamma0*SIG
+ return(dmnorm(TS,rep(0,n),SIG,log=TRUE))}
> LogL=function(A) -GlobalLogLik(A,TS=Z)
> optim(c(.1,.1,1),LogL)
$par
[1] 0.2584144 0.6826530 1.0669820

$value
[1] 298.8699

$counts
function gradient 
      86       NA 

$convergence
[1] 0

$message
NULL

Here, the values that minimize the likelihood are rather close to the ones used to generate our sample. And if we run this algorithm on 1,000 samples, we can see that those estimates are fine,

I could not find other ideas, to estimate those parameters. I guess we can use the partial autocorrelation function, since we have relationships that can be related to Yule-Walker equations for https://latex.codecogs.com/gif.latex?AR(p) time series.

Journalisme Scientifique

Il y a quelques jours, Passeur de Sciences a annoncé que le site du Monde.fr allait abriter plusieurs nouveaux blogs de science. Mon voisin @tomroud a parfaitement résumé la situation en 140 caractères, “plus de blogs sur les sciences sur le monde.fr http://passeurdesciences.blog.lemonde.fr/…. Mais toujours pas (vraiment) de blogs de scientifiques” (il ajoutait quelques minutes plus tard “il semble donc que l’avenir du blog sur les sciences soit le journalisme scientifique. Pas sûr que ce soit un progrès…“), et @marc_rr a poursuivi la discussion sur son blog – tout se passe comme si – avec un long débat passionnant. N’ayant pas l’habitude de poster des commentaires, je vais plutôt mettre un billet sur mon blog, car je plussoie à la vision de @tomroud !

Je suis en train de finaliser depuis quelques jours un article sur les blogs académiques, et à plusieurs reprises, je me suis interrogé sur la distinction entre l’activité de bloggeur académique, et l’activité de journaliste scientifique. Et le débat des derniers jours (je devrais dire des derniers heures) m’a beaucoup éclairé. Je pourrais faire beaucoup de reproches à @passeursciences (et je l’avais déjà fait une fois, dans Martingale et Journalisme Scientifique, suite à la publication d’un article qui, en quelques lignes, avait balayé le travail de centaines de scientifiques, dont deux chercheurs qui auraient un prix Nobel quelques mois plus tard pour ces travaux, justement). Et beaucoup de compliments (j’adore la chronique “Improbablologie“, même si j’ai pris l’habitude de la lire en la version originale sur http://www.improbable.com/, car j’adore l’humour de @marcabrahams). Et dans un commentaire sur le blog de @marc_rr, j’ai eu l’impression que le débat se simplifiait, et que la distinction entre scientifiques et journaliste devenait plus claire,

“Les chercheurs-blogueurs : vous allez sans doute me trouver sévère mais si je dis que j’ai du mal à trouver des chercheurs qui bloguent a/ très régulièrement (au moins une ou deux fois par semaine) b/ en français et sans fautes c/ avec un excellent niveau de vulgarisation (donc sans s’adresser en jargon aux happy few) et d/ sur l’actualité, c’est parce que cela correspond à la réalité. Mais c’est aussi sans doute parce que mes critères sont trop élevés ou trop journalistiques” (a continuer…)

Le problème du jargon scientifique est un vieux problème, qui semble atrocement gêner les journalistes, car ils semblent ne pas le maîtriser. Mais je me suis toujours énervé contre cette tentative de nivellement par le bas (ce qui fait que je prend le clavier aujourd’hui). Cette distanciation forcée enlève beaucoup, et fait rapidement perdre toute crédibilité. Dans son Histoire de la Lecture, Alberto Manguel citait un poète gaulois, Ausone, qui écrivait

Tu as acheté des livres et rempli des rayons, ô amoureux des Muses.
Cela signifie-t-il que tu es désormais savant ?
Si tu achètes aujourd’hui des instruments à cordes, plectre et lyre :
Crois-tu que demain le royaume de la musique t’appartiendra ?

Effectivement, je pense qu’un passionné d’astronomie ne se contente pas de lire Ciel et Espace, il va s’acheter un télescope, et essayer de voir les étoiles par lui même. Et celui qui s’intéresse aux mathématiques ira lire Choux Romanesco etc, avec tout le jargon qui va avec, et voudra aller plus loin. Lorsqu’on a vu les figures d’Esscher, avec les enfants, on a essayé de comprendre pourquoi ça marche ! Alors on a pris nos crayons, nos compas, et on s’est lancé. La science, ça se vit. Je pense que c’est la grande distinction entre les journalistes scientifiques et les bloggers académiques. Dans un vieux billet, sur l’importance du do-it-yourself, j’avais justement expliqué qu’il était important que les gens mettent les mains à la pâte. Je peux citer, par exemple, un autre billet (né d’une discussion avec @imparibus) à l’époque des élections présidentielles. J’avais mis en ligne des graphiques et des codes informatiques pour faire un peu de prévision. Au lieu d’attendre que quelqu’un d’autre le fasse. Dans quelques jours, je dois participer a une discussion (dite grand public) sur le big data. Et un de mes points est qu’il est important de comprendre ce qui se passe, comment fonctionnent les algorithmes qui, soit disant, nous gouvernent. L’affiche de la conférence (j’en reparlerais bientôt sur le blog) associe big data et big brother, et je pense que cette association est (partiellement) fausse : on pense à Big Brother si on se sent dépossédé, et je ne pense pas que ce soit le cas. A condition de comprendre ce qui peut se faire avec du big data, et ce qui reste de l’utopie… Il faut se lancer, et s’approprier le jargon, au contraire ! Je trouve cette démission face à la difficulté incroyable !

Quant à parler de happy-few, c’est tellement méprisant ! Au risque de choquer, il y a du monde qui veut des choses techniques, avec du jargon dedans ! Dans l’article mentionné au début, @passeursciences se félicitait d’avoir 20 millions de pages vues en un peu plus de décembre 2011 (et pas lues comme le dit l’article, malheureusement, on ne saura jamais trop ce que les gens font quand ils arrivent sur nos pages). Je l’en félicite. Personnellement, je suis loin derrière, car depuis décembre 2012 (j’ai fait ma migration assez tardivement) j’en totalise à peine 2 millions. Étant statisticien, je serais le premier à éviter de donner trop de poids à ces chiffres (surtout que je publie, à l’occasion, en anglais, ce qui peut faire augmenter le nombre de lecteurs, et beaucoup de mes articles sont resyndiqués sur des blogs comme , r-bloggers, architects.dzone ou statsblogs – parmi ceux que je connais – ce qui fait baisser le nombre de lecteurs), mais quand on regarde rapidement, ça veut dire que j’ai juste 5 fois moins de pages vues qu’un blog hébergé par le plus grand quotidien français !? En mettant des maths et des formules de codes dans presque tous mes billets ?! J’ai du mal a comprendre l’idée de happy-few, je suis désolé… Cela dit, j’espère que ceux qui viennent sont effectivement happy

PS: sur le français sans faute, je suis malheureusement tout seul sur mon blog, et j’ai beaucoup de mal à me relire moi-même. Mais promis, un jour je me payerais le luxe d’avoir quelqu’un pour relire mes articles avant que je ne les mette en ligne ! Ah ah ah, on peut toujours rêver…