# What is a Linear Trend, by the way?

I had a very strange discussion on twitter (yes, another one), about regression curves. I think it started with a tweet based on some xkcd picture (just for fun, because it was New Year’s Day)

There were comments on that picture, by econometricians, mainly about ‘significant’ trends when datasets are very noisy. And I mentioned a graph that I saw earlier, a couple of days ago

Let us reproduce that graph (Roger kindly sent me the dataset)

db=data.frame(year=1990:2016, ratio=c(.23,.27,.32,.37,.22,.26,.29,.15,.40,.28,.14,.09,.24,.18,.29,.51,.13,.17,.25,.13,.21,.29,.25,.2,.15,.12,.12)) library(ggplot2)

The graph is here (with the same aesthetic conventions as Roger’s initial graph, i.e. using some sort of barplot)

ggplot(db, aes(year, ratio)) + geom_bar(stat="identity") + stat_smooth(method = "lm", se = FALSE) My point was that we miss the ‘confidence band’ of the regression

In R, at least, it is quite natural to get (and actually, it is the default version of the graph function)

ggplot(db, aes(year, ratio)) + geom_bar(stat="identity") + stat_smooth(method = "lm", se = TRUE) It is hard to claim that the ‘regression line’ is significant (in the sense “significantly non horizontal”). To be more specific, if we look at the output of the regression model, we get

summary(lm(ratio~year,data=db)) 

Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 9.158531 4.549672 2.013 0.055 . year -0.004457 0.002271 -1.962 0.061 . --- Signif. codes: 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(which is exactly what Roger used in his graph to plot his red straight line). The p-value of the estimator of the slope, in a linear regression model is here 6%. But I found Roger’s point puzzling

First of all, let us get back to a more standard graph, with a scatterplot, and not bars,

ggplot(db, aes(year, ratio)) + stat_smooth(method = "lm") + geom_point() Here, we observe points $\{y_{1990},y_{1991},\cdots,y_{2016}\}$. In order to draw that blue line, we assume (Econometrics 101, actually) that those observations are realizations of random variables $\{Y_{1990},Y_{1991},\cdots,Y_{2016}\}$. Randomness here does not come from a survey, or from ‘balls in an urn’. Randomness is because hurricanes and floods are themselves seen are realizations of random events. Yes, there might be measurement errors, but that’s not where randomness comes from (here). When we talk about ‘randomness’, it should be related to ‘model error’ i.e. the error we make if we consider a linear model (here), that is

$$?Y_t=\beta_0+\beta_1t+\varepsilon_t$$

Even if observations are not obtained from balls in an urn, there is some kind of randomness here. Randomness means that we might have errors (random errors) around the estimated value (that is on the blue curve), $y_t=\widehat{y}_t+\widehat{\varepsilon}_t$. One might consider a nonlinear model to reduce the error,

ggplot(db, aes(year, ratio)) + geom_point() + geom_smooth()

but in the case, the danger is to overfit So yes, when we fit a linear model, there is always some kind of randomness, and it is possible to get a ‘confidence band’, that will be very useful for predictions (e.g. for reinsurance purpose here).

# De la créativité des gangsters

Pendant mon séjour récent en Nouvelle Angleterre, j’ai survolé le livre de Leonard Mlodinow, the drunkard’s walk. Et au hasard de mes lectures, je suis tombé sur la petite histoire suivante Autrement dit, en utilisant les cinq derniers chiffres d’une quantité économique comme le us treasury balance, on aurait un générateur de nombre aléatoire…
Par contre la suite est un peu plus surprenante,  La loi de Benford (que j’avais pu évoquer ici ou ) parle des premiers chiffres, mais cette fois on parle deslast five digits. Donc visiblement l’évoquation n’est pas pertinente ici. Mais qui sait ? Ca reste malgré tout une histoire intéressante. Considérons – histoire de tester cette légende – les deux sources de données suivantes,http://treasurydirect.gov/ et http://economagic.com/. Sur ce dernier site, un petit travail de mise en forme des données est nécessaire.

b1=read.table("http://freakonometrics.free.fr/debtus1.txt",
X1=as.character(b1$Dollar.Amount) n1=nchar(X1) Y1=substr(X1,n1-8,n1-3) X1=as.numeric(substr(Y1,1,2))*1000+as.numeric(substr(Y1,4,6)) x=X1/100000 X2=b2$DEBT
Y2=trunc(as.numeric(X2))
X2=as.character(Y2)
n2=nchar(X2)
Y2=substr(X2,n2-4,n2)
y=as.numeric(Y2)/100000
y=y[y<1]

Pour rappel, un générateur aléatoire (standard) vérifie deux propriétés importantes

• les nombres doivent être tirés suivant une loi uniforme sur [0,1], i.e. ici, si on divise les nombres à 5 chiffres par 10000,
• les tirages doivent être indépendants entre eux.

La première propriété semble assez naturelle, et correspond à l’histoire racontée dans un commentaire posté ici (expliquant comment un casino avait été au bord de la faillite car une roulette faisait sortir certains chiffres trop souvent, et j’essayais de comprendre comment utiliser l’information qu’un chiffre sort plus souvent). La seconde est probablement encore plus importante.

• Visualisation de la distribution

La première idée est de visualiser la densité de nos séries de chiffres. Pour éviter les problèmes de bord (et comme c’est juste en introduction) on va utiliser un histogramme, et pas une estimation à noyau.

hist(x,col="red")
hist(y,col="blue")

On obtient pour la première série la courbe rouge, et pour la seconde la courbebleue, On note qu’a priori, pour la première série, l’hypothèse d’uniformité n’est peut être pas la plus réaliste…

• Test de Kolmogorov-Smirnov

On peut aussi mettre en œuvre un test de Kolmogorov-Smirnov afin de tester si la loi uniforme est adaptée.:

> ks.test(x,"punif")

One-sample Kolmogorov-Smirnov test

data:  x
D = 0.1047, p-value = 0.01645
alternative hypothesis: two-sided

> ks.test(y,"punif")

One-sample Kolmogorov-Smirnov test

data:  y
D = 0.0456, p-value = 0.3581
alternative hypothesis: two-sided

On retrouve ici confirmée l’intuition précédante: la loi uniforme est pertinente pour la seconde série, pas la première.

• Les autocorrélations de la série

Travaillons uniquement sur la seconde série. On peut étudier l’autocorrélation de notre série de nombres, ou peut-être un peu plus malin, sur les quantiles gaussiens associés (les autocorrélations étant intéressantes pour les séries gaussiennes),:

plot(acf(y))
plot(acf(qnorm(y)))

ie. pour la série brute et pour la série normalisée, Bref, on pourrait être tenté de valider l’hypothèse d’indépendance entre les tirages.
• Run test (de Wald–Wolfowitz)

L’idée est de comparer une série de chiffres à la médiane, s’ils sont plus grands, on note + (ou A) et sinon – (ou B). On crée alors une série du genre “+++−−++−−++++++−−−” et on compte les séries de + et les séries de -, les runs,

library(lawstat)
runs.test(y,plot=TRUE)

Runs Test - Two sided

data:  y
Standardized Runs Statistic = -0.2462, p-value = 0.8055 Bref, la légende me semble à prendre avec des pincettes (car fonction de la source considérée), même si l’idée est intéressante (si l’on met de côté les aspects d’aléa moral). Et l’analyse sur la loi de Benford ne semble pas valide: les derniers chiffres sur les grands nombres ne se comportent pas du tout comme les premiers.