# Forecasts with ARIMA Models

In our time series class this morning, I was discussing forecasts with ARIMA Models. Consider some simple stationnary AR(1) simulated time series

> n=95
> set.seed(1)
> E=rnorm(n)
> X=rep(0,n)
> phi=.85
> for(t in 2:n) X[t]=phi*X[t-1]+E[t]
> plot(X,type="l")

If we fit an AR(1) model,

> model=arima(X,order=c(1,0,0),
+             include.mean = FALSE)
> plot(P$pred) > lines(P$pred+2*P$se,col="red") > lines(P$pred-2*P$se,col="red") > abline(h=0,lty=2) > abline(h=2*P$se[20],lty=2,col="red")
> abline(h=-2*P$se[20],lty=2,col="red") we observe the exponential decay of the forecast towards 0, and the increasing confidence interval (where the variance increases, from the variance of the white noise to the variance of the stationnary time series). Plain lines are conditional forecast (given our latest observation, since the AR(1) is a first order Markov process), and dotted lines are unconditional. Let us store some values, to use them as benchmark > s=P$se[20]
> y=P$pred If we fit a MA(1) model > model=arima(X,order=c(0,0,1), + include.mean = FALSE) > P=predict(model,n.ahead=20) > plot(P$pred)
> lines(P$pred+2*P$se,col="red")
> lines(P$pred-2*P$se,col="red")
> abline(h=0,lty=2)
> abline(h=2*s,lty=2,col="red")
> abline(h=-2*s,lty=2,col="red")
> lines(y,col="grey")

after two lags, the forecast is null, and the (conditional) variance remains constant. But if we consider a moving average process with a longer order,

> model=arima(X,order=c(0,0,14),
+             include.mean = FALSE)
> plot(P$pred) > lines(P$pred+2*P$se,col="red") > lines(P$pred-2*P$se,col="red") > abline(h=0,lty=2) > abline(h=2*s,lty=2,col="red") > abline(h=-2*s,lty=2,col="red") > lines(y,col="grey") we get an output that can be compared with the AR(1) processes. Which makes sense since our AR(1) process can also be seen as a MA(∞), with infinite order. But if we think that our time series is not stationary, an we fit an integrated model > model=arima(X,order=c(0,1,0), + include.mean = FALSE) > P=predict(model,n.ahead=20) > plot(P$pred)
> lines(P$pred+2*P$se,col="red")
> lines(P$pred-2*P$se,col="red")
> abline(h=0,lty=2)
> abline(h=2*s,lty=2,col="red")
> abline(h=-2*s,lty=2,col="red")
> lines(y,col="grey")

we observe the (standard) martingale property: the forecast is flat, and the confidence interval keeps increasing, and actually, the variance increases towards infinity (at a linear rate). So one should be very careful when differentiation a time series… it will have a huge impact on the forecasts….

# Temperatures Series as Random Walks

Last year, I did mention in a post that unit-root tests are dangerous, because they might lead us to strange models. For instance, in a post, I did obtain that the temperature observed in January 2013, in Montréal, might be considered as a random walk process (or at leat an integrated process). The code to extract the data has changed (since the website has been updated), so here, we use

library(RCurl)
library(XML)
options(RCurlOptions = list(useragent = "R"))
HEURE=0:23
extracttemp=function(Y,M,D){
url=paste(
"http://climate.weather.gc.ca/climateData/hourlydata_e.html?timeframe=1&Prov=QC&StationID=5415&Year=",Y,"&Month=",
M,"&Day=",D,sep="")
wp <- getURLContent(url)
doc <- htmlParse(wp, asText = TRUE)
docName(doc) <- url
basejour=data.frame(Year=Y,Month=M,Day=D,
Hour=HEURE,Temp=as.numeric(as.character(data.frame(tmp[2])[,2]))[2:25])
return(basejour)}
B=NULL
for(y in 1955:2013){
for(d in 1:31){
B=rbind(B,extracttemp(y,1,d))}}

Here are all the temperatures observed, and 2013,

plot(B$X,B$Temp,cex=.5,col="light blue",xlab="January, in Montreal",ylab="Temperature (Celsius)")
I=which(B$Year==2013) lines(B$X[I],B$Temp[I],col="red") In the previous post, one test only was used, and one year was considered. I was wondering if this behavior was observed only with temperature of 2013 (or not), and how the other tests (mentioned in a previous post too) were performing. I might need a function, because those tests cannot be used if there is a missing value, even only one. So I did use the value observed one hour before (just to make sure that the tests can be done) correcty=function(Y){ I=which(is.na(Y)) if(length(I)==0){Yc=Y} if(length(I)>0){Yc=Y;for(i in I) Yc[i]=Yc[i-1]} return(Yc) } Now, we can compute the p-values, for all the years, and the three different three (keeping in mind that two test if the series is non-stationary, and one if the series is stationary) DF=matrix(NA,2013-1954,3) library(urca) for(y in 1955:2013){ Z=B$Temp[which(B$Year==y)] Zc=correcty(Z) DF[y-1954,2]=as.numeric(pp.test(Zc)$p.value)
DF[y-1954,1]=as.numeric(kpss.test(Zc)$p.value) DF[y-1954,3]=as.numeric(adf.test(Zc)$p.value)
}

Visually, if red means stationary, and blue means non-stationary, we get

DFP=DF
DFP[,1]=DF[,1]<.05
DFP[,2:3]=DF[,2:3]>.05
library(RColorBrewer)
CL=brewer.pal(6, "RdBu")
plot(0:1,0:1,xlim=c(1950,2015),ylim=c(0,3),axes=FALSE,xlab="",ylab="")
axis(1)
text(1952,.5,"KPSS")
text(1952,1.5,"PP")
for(y in 1955:2013){
for(i in 1:3){
polygon(y+c(-1,-1,1,1)/2.2,i-.5+c(-1,1,1,-1)/2.2,col=CL[1+(DFP[y-1954,i]==1)*5],border=NA)}}

Quite frequently, we conclude that the temperature is a random walk. Which does not make sense (from a physical point of view). But again, it might come from the fact that temperature are stationary, but with some fractional behavior (as suggested in the previous post).

# Central Limit Theorem

This week, in the MAT8595 course, before proving Fisher-Tippett theorem, we will get back on the proof of the Central Limit Theorem, and the class of stable distribution (in Lévy’s sense). In order to illustrate the problem of heavy tails, on the behavior of the mean, consider a sequence of i.i.d. Gaussian random variables $X_i$‘s. On top, we visualize the sequence, and below, we visualize the associate random walk

$S_n=\sum_{i=1}^n X_i$

(the central limit theorem will give a limiting distribution for $n^{-1}S_n$ in the case where the variance of the $X_i$‘s is finite)

If we consider a sequence of i.i.d. random variables $X_i$‘s whith heavier tails (possibly with infinite variance), we can still define $S_n$, but as we can see below, $S_n$ can be quite heratic.

As we will see this Thursday, the key to derive stable distribution for the central limit theorem, or possible limiting distributions for the maximum is Cauchy’s function equation. I strongly recommand to look at the proof.

# Martingale et journalisme scientifique

Être journaliste scientifique ne doit pas être facile. J’imagine qu’il faut être à l’écoute des nouvelles scientifiques, et d’informer, aussi justement que possible. C’est ce que fait avec brio Pierre Barthélémy (aka @PasseurSciences sur Twitter) dans sa chronique hebdomadaire sélection scientifique de la semaine. Il essaye ainsi de parler de sciences dans un journal qui a (trop) souvent confondu science et technologie (technologie étant aussi souvent un mot savant utilisé pour masquer de la publicité pour des appareils technologiques avancés). Et comme l’espace qui lui est imparti est restreint (c’est le moins qu’on puisse dire), la plupart des informations passent par son blog http://passeurdesciences.blog.lemonde.fr/, qui est riche ! Sylvestre Huet fait aussi ce travail pour Libération, avec son blog http://sciences.blogs.liberation.fr/ (là encore, on peut lire bien plus de choses en ligne que que dans la version papier, ce qui est d’autant plus intéressant pour les lecteurs de l’autre coté de l’océan atlantique). Sylvestre Huet arrive à parler géologie, puis démographie, avant d’évoquer dans un article truffé de références les réformes des instances qui régentent la recherche en France. Oui, les deux blogs sont malgré tout très français. Pour une vision plus internationale, on pourra citer l’admirable travail de Pascal Lapointe, par exemple (aka @paslap) sur http://sciencepresse.qc.ca/.

Mais si je commence à mettre les pieds dans le plat, je poserais la question de la légitimité d’un journaliste scientifique, à l’heure où la blogosphère scientifique explose. Par exemple, si je veux apprendre des choses en physique, ou en biologie, je sais que le blog de Tom Roud (que l’on peut suivre sur Twitter @tomroud) http://tomroud.cafe-sciences.org est une source infinie d’information (ok, peut-être pas “infinie” car comme nous tous, il a une vie en dehors de son blog, mais disons qu’il y a de quoi lire). En sciences humaines, il existe des centaines de blogs hébergés sur http://hypotheses.org/ tenus par des universitaires, ou http://www.cafe-sciences.org/ pour des blogs de sciences, en français. Les bloggers scientifiques s’expriment – le plus souvent – en restant dans leur champ d’expertise. Et c’est tant mieux. Le journaliste scientifique, lui, se doit de parler de tous les sujets scientifiques. Et c’est là que l’exercice devient délicat, car le journaliste se doit d’être critique. Compte tenu de la concurrence féroce qui existe dans le monde académique, on nous pousse à faire croire qu’on vient d’écrire l’article qui va révolutionner la science, sinon le monde. Que notre approche est novatrice, et qu’en plus, on vient de montrer ce qui pourrait ce qui pourrait être le sain Graal dans notre communauté. C’est le travail du journaliste de se demander si c’est vrai….

J’en arrive au point de mon billet. Le dernier article  le hasard, martingale boursière ? (en ligne sur http://lemonde.fr/sciences/…) de Pierre Barthélémy m’a un peu agacé. Pas sur le thème abordé, mais l’impression générale qui m’est restée après avoir lu l’article (samedi matin). Pour l’histoire complète, c’est ma femme qui est tombée la première sur l’article, et qui a été étonnée de lire un article pareil dans Le Monde, en tant que mathématicienne (et probabiliste). Et quand je l’ai lu, c’est en tant qu’économiste que je suis resté sans voix. En résumé (trop succinct, j’en convient), on nous résume l’article de physiciens italiens qui “expliquent qu’au grand dam des traders et autres analystes financiers, en quête perpétuelle de justifications rationnelles aux fluctuations des cours boursiers, les marchés demeurent obstinément imprévisibles“. Pour parcourir l’article mentionné, mais non cité – probablement car il s’agit de la version papier ? – l’article est en ligne sur http://arxiv.org/1303.4351  (trouvé via http://improbable.com/…). On retrouve cette idée dans le texte sous la forme suivante “our main result, which is independent of the market considered, is that standard trading strategies and their algorithms, based on the past history of the time series, although have occasionally the chance to be successful inside small temporal windows, on a large temporal scale, perform on average not better than the purely random strategy, which, on the other hand, is also much less volatile. In this respect, for the individual trader, a purely random strategy represents a costless alternative to expensive professional financial consulting, being at the same time also much less risky, if compared to the other trading strategies.” On apprend dans l’article paru dans Le Monde que “pour Pluchino et compagnie, le CAC 40 est une loterie, ce que les spécialistes ne veulent pas admettre“. Damned, rien que ça ? Pour les amateurs de lyrisme (il doit y en avoir dans les lecteurs du monde, moins parmi les lecteurs plus scientifiques qui aiment les phrases moins alambiquées), je passe les moments où Pierre Barthélémy est clairement plus journaliste que scientifique, comme lorsqu’il nous explique que ces quatre chercheurs “viennent d ‘attaquer à l’acide de l’aléatoire” les modèles financiers. Cela dit, cela a le mérite de mettre les points sur les i: ce sont des physiciens. Pas des économistes.

Je pense qu’il y a des milliers de raisons d’attaquer les traders (si ce sont bien les “spécialistes” visés). Mais qu’on ne les prenne pas pour plus bêtes qu’ils ne le sont (je le dis d’autant plus volontiers que j’en ai formé un paquet, à l’ENSAE ou à Polytechnique, au cours des 10 dernières années).

Avant de revenir sur la petite histoire de martingale financière (tel que le résume le titre de l’article) faisons un court détour par celle brillamment racontée dans Mansuy (en ligne sur http://math.harvard.edu/~ctm/… et traduit en anglais dans http://jehps.net/Mansuy.pdf) sur l’origine de la notion de martingale. D’un point de vue mathématique (disons comme propriété des processus stochastiques), il faut remonter aux travaux de Joseph Bernstein, Paul Lévy, Émile Borel, et surtout Joseph Leo Doob, au milieu au XXième siècle. Cela dit, le mot est plus ancien: il entre dans le dictionnaire de l’Académie Française en 1762. Jouer à la martingale, c’est jouer toujours tout ce que l’on a perdu.  Mais quelques années plus tot, on pouvait déjà trouver le mot sous la plume de l’Abbé Prévost (le jeu qu’il décrit en 1750 comme variante du jeu du pharaon est aussi appelé martingale d’Alembert) Cela dit, l’origine que je préfère est de relier le mot martingale à l’expression provençale “jouga a la martegalo” qui signifierait jouer de manière incompréhensible, absurde, comme l’évoque Frédéric Mistral dans Lou Tresor dòu Felibrige ou dictionnaire provençal-français. On retrouve une origine proche dans le dictionarie of the French and English tongues de Randle Cotgrave, datant de 1611, qui mentionne l’expression “à la martingale” avec le sens absurdly, foolishly, untowardly, grossely, rudely, in the homeliest manner. Il cite même l’usage de l’expression philosopher à la martingale (sans citer Bernard Henri Levy, mais il semble que ce soit l’idée). Cela dit, même dans Lapinot on parle de martingale,

Bref, à partir de ce comportement incompréhensible, voire absurde, des économistes vont définir une notion importante en finance, que l’on appellera efficience (qui est un mot dangereux car il évoque aux oreilles de tous les économiste une notion d’optimalité au sens de Pareto, on pourra relire Malkiel (2003) The Efficient Market Hypothesis and Its Critics sur ce point). L’idée d’utiliser une marche aléatoire pour modéliser le cours d’un actif est ancienne. Par exemple (je ne remontrais pas au XVIème siècle, mais http://e-m-h.org/history.html le suggère) on pourra relire les travaux de Bachelier datant du début du XXième siècle, mais comme le dit la légende, Bachelier a été peu lu, à l’époque (en tous les cas par des économistes). Publié quelques années plus tard, on pourra relire Cowles & Jones (1937), Some A Posteriori Probabilities in Stock Market ou surtout Samuelson (1965) Proof That Properly Anticipated Prices Fluctuate Randomly. Cet article est passionnant, si on prend le temps de le lire “There is no way of making an expected profit by extrapolating past changes in the future price, by chart or any other esoteric devices of magic or mathematics.” Paul Samuelson est d’une modestie remarquable pour un chercheur, “I have not here discussed what the basic probability distributions are supposed to come from. In whose mind are they ex ante? Is there any ex post validation of them? Are they supposed to belong to the market as a whole? And what does that mean? Are they supposed to belong to the “representative individual”, and who is he? Are there some defensible or necessitous compromises of divergent anticipations patterns? Do price quotations somehow produce a Pareto-optimal configuration of ex ante subjective probabilities? This paper has not attempted to pronounce on these interesting questions“. Il peut y avoir des bulles, des processus non Gaussien, un peu tout ce qu’on peut imaginer, le point important étant qu’il est impossible d’utiliser le passé pour prédire le futur, en finance.

Oui, depuis 50 ans, les économistes savent que la marche aléatoire peut être un bon modèle pour décrire le prix des actifs. Ou pour être plus précis (et rendre à Paul Samuelson ce que Paul Samuelson a dit le premier) les martingales. Mais il faudra surtout attendre les travaux d’Eugène Fama, en particulier Fama (1965) The Behavior of Stock Market Prices et Fama (1970) Efficient Capital Markets: A Review of Theory and Empirical Work, qui vont de manière définitive marquer le début de l’utilisation des martingales en finance. Eugène Fama y retient trois notions d’efficience (pour aller plus loin que le fameux “A market in which prices at any time “fully reflect” available information is called “efficient”“)

1. Expected Returns or “Fair Game” Models
2. Submartingale Models
3. Random Walk Model

(en 1965, l’efficience devait être reliée, pour Eugène Fama, à la notion de marche aléatoire, alors que Paul Samuelson utilisait déjà la notion de martingale). Ces modèles disent tous (on est en 1970) qu’il est inutile d’utiliser l’information passée ou présente “to predict the future in a way which makes expected profits greater than they would be under a naive buy-and-hold model“. N’en déplaise à Pierre Barthélémy , depuis presque 50 ans ce principe est énoncé dans les cours d’asset pricing. Bachelier disait la meme chose en 1900, “l’espérance mathématique du spéculateur est nulle” dans Théorie de la spéculation. Et cela sera confirmé quelques années plus tard par les résultats empiriques de Cowles et Jones. Un an plus tard sera publié Black (1971) Implications of the random walk hypothesis for portfolio management, et trois ans plus tard, Hagerman & Richmond (1973) Random Walks, Martingales and the OTC. La littérature économique va se multiplier pendant les années 70, à tel point qu’en 1978, Michael Jensen, alors professeur à Harvard écrivait “I believe there is no other proposition in economics which has more solid empirical evidence supporting it than the Efficient Market Hypothesis” (cité dans http://economist.com/14030296). Plus récemment, à la fin des années 80 était publié le remarquable LeRoy (1989) Efficient Capital Markets and Martingales, qui étoffait la critique initiée dans LeRoy (1976). Efficient capital markets: Comment. LeRoy et Samuelson ont été les premiers à parler de martingale pour modéliser les prix des actifs. On continue (avec les articles qui font référence, et que tous les étudiants qui ont suivi un cours d’asset pricing ont lu) ? L’année suivante était publié Lehman (1990) Fads, Martingales, and Market Efficiency. Les martingales sont des outils incroyablement riches. Ils ne sont pas équivalent à une marche aléatoire: on peut avoir une martingale, sans avoir de marche aléatoire (par exemple avec un processus ARCH, on pourra relire sur le sujet les notes de cours Predictability of Asset Returns, en particulier le premier paragraphe, ou en français, avec des cours de licence, par Francis Diener,  http://math.unice.fr/~diener/…). Pour ceux qui veulent moins de formalisme, une martingale (en simplifant outrageusement) c’est, comme l’explique Nicolas Poupon,

Mais ces articles posent essentiellement les bases de la théorie financière. Qu’en est-il des aspects empiriques ? Andrew Lo et Craig MacKinlay ont publié un livre, en 2001, qui recense plusieurs articles sur le sujet, A Non-Random Walk Down Wall Street. Qui remet en cause l’idée de marche aléatoire, moins celle de marginale. On pourra aussi penser à Beechey Gruen & Vickrey (2000) et leur étude The Efficient Markets Hypothesis: A Survey, en ligne sur la Reserve Bank of Australia. Enfin, Pour une méthodologie plus proche de celle évoquée dans l’article en ligne sur arxiv, on pourra lire Martingales, the Efficient Market Hypothesis, and Spurious Stylized Facts de Joseph McCauley, Kevin Bassler et Gemunu Gunaratne (si on souhaite aussi utiliser des propriétés de mémoire longue). Moralité ? Non, ce n’est pas nouveau que l’on teste l’hypothèse d’efficience et de martingale. Et si je cite des articles académiques, il faut  ajouter que les journalistes économiques connaissent également cette littérature. Il y a 10 ans, Justin Fox posait la même question que nos chercheurs dans un article dans Fortune intitulé Is The Market Rational ? No, say the experts. But neither are you–so don’t go thinking you can outsmart it. On pourra aussi relire Efficiency and beyond (mentionné auparavant) qui revient sur la notion d’efficience des marchés. On y retrouve que les spécialistes savent tout ca: “on such ideas, and on the complex mathematics that described them [i.e martingales], was founded the Wall Street profession of financial engineering“. Enfin, en 2009, Richard Thaler faisait quelques rappels sur la difficulté d’interpréter correctement l’hypothèse d’efficience des marchés dans Markets can be wrong and the price is not always right. Pour conclure sur l’histoire des martingales en finance, en 2010, Fama répondait d’ailleurs de manière très claire à la question dans une interview – interview with Eugene Fama – lorsque John Cassidy lui demandait “the fundamental insight of the efficient market hypothesis [is] that you can’t beat the market ?“, et qu’il répondait “Right—that’s the practical insight. No matter what research gets done, that one always looks good“.

Pour revenir sur l’article de Pierre Barthélémy, je trouve qu’écrire “le CAC 40 est une loterie, ce que les spécialistes ne veulent pas admettre” est incroyablement méprisant envers des centaines, pour ne pas dire des milliers de chercheurs en mathématiques financières, et en économétrie de la finance. Mais rassure toi Pierre, je vais continuer à lire tes chroniques. C’est juste que ma grand mère est encore abonnée au Monde, et j’aimerais bien qu’elle continue de croire que je peux sauver l’humanité (et non pas que je refuse d’admettre ce que des chercheurs italiens viennent d’établir, même si cela confirmerait – à ses yeux – que je peux etre têtu comme une mule quand je veux).

Maintenant, pour conclure sur l’article mentionné, et sur l’analyse faite, on ne peut pas dire que les conclusions soient révolutionnaires. Disons, pour quiconque ayant lu un peu de littérature économique dans les 50 dernières années. En fait, pour être honnête, je pense qu’on peut dire que Alessandro Pluchino et Andrea Rapisarda se sont amusé (on inclura Biondo et Helbing). Sans plus de prétention. Pour ceux qui ne se souviennent pas, ils avaient eu un IGNobel en 2011, pour leur étude sur organizations would become more efficient if they promoted people at random. J’avais aussi beaucoup aimé leur article sur les aspects computationnels du principe de Peter, en ligne sur http://arxiv.org/0907.0455. Quand je lis leurs études, j’ai l’impression qu’ils s’amusent. Et le terme n’est absolument pas péjoratif dans ma bouche, loin de là. C’est un peu ce que disait Lucile Quillet la semaine passée dans http://etudiant.lefigaro.fr/…  lorsqu’en mentionnant une étude de Baptiste Coulmont, elle écrivait “c’est le jeu auquel c’est amusé Baptiste Coulmont sur son blog“. Je ne sais pas si le terme était voulu, si c’était basé sur le fait qu’un blog, c’est fait pour jouer… mais j’ai toujours revendiqué que je faisais de la recherche pour mon plaisir, pour m’amuser. Et comme me le rappelait Baptiste, le mot schola signifie école, en latin, et son étymologie grecque est le mot σχολή, signifiant loisir.

# A random walk ? What else ?

Consider the following time series,

What does it look like ? I know, this is a stupid game, but I keep using it in my time series courses. It does look like a random walk, doesn’t it ? If we use Philipps-Perron test, yes, it does,

> PP.test(x)

Phillips-Perron Unit Root Test

data:  x
Dickey-Fuller = -2.2421, Truncation lag parameter = 6, p-value = 0.4758

If we look at the autocorrelation function, we do observe some persistence,

> acf(x,100)

Perhaps this persistence can be related to long range dependence, or to some fractional random walk. A natural idea could be estimate Hurst parameter, using for instance Beran (1992) estimator – based on Whittle (1956) – where we assume that the autocorrelation function satisfies

$\rho(h)=\text{cor}(X_t,X_{t-h})\sim \alpha |h|^{2H-2}$

as $h\rightarrow\infty$ for some $H\in(1/2,1)$ (the so called Hurst index). But here, we start to observe unexpected ouputs,

> library(longmemo)
> (d  <- WhittleEst(x))
'WhittleEst' Whittle estimator for  fractional Gaussian noise ('fGn');	 call:
WhittleEst(x = x)
time series of length  n = 759.

H = 0.9899335
coefficients 'eta' =
Estimate Std. Error z value   Pr(>|z|)
H 0.98993350 0.02468323 40.1055 < 2.22e-16
<==> d := H - 1/2 = 0.49 (0.025)

$vcov : num [1, 1] 0.000609 ..- attr(*, "dimnames")=List of 2 .. ..$ : chr "H"
.. ..$: chr "H"$ periodogr.x: num [1:379] 1479.3 1077.3 371.7 287.2 51.2 ...
\$ spec       : num [1:379] 62.5 31.7 21.3 16.1 12.9 ...

or more precisely some non-expected values for Hurst parameter, which should be in $(1/2,1)$

> confint(d)
2.5 %   97.5 %
H 0.9415553 1.038312

Oops, perhaps, we did miss something, because it looks like there is extremely strong persistence on our time series,

> plot(d)

It is probablty time to ask where I found that series… To be honest, I did borrow  it from a great canadian website http://climate.weatheroffice.gc.ca/climateData/. For instance, it you want the temperature we did experience a few days ago, you can use

> Y=2013
> M=1
> D=25
> url=paste(
"http://climate.weatheroffice.gc.ca/climateData/hourlydata_e.html?
timeframe=1&Prov=QC&StationID=5415&hlyRange=1953-01-01|2013-02-
01&Year=",Y,"&Month=",M,"&Day=",D,sep="")
> page=scan(url,what="character")

Yes, that series is the temperature we did experience in Montréal last month (hourly time seies). On the graph below, you can actually compare it with temperature experienced in Januarys over the past 60 years,

So it is not that surprising to see long range dependence models appearing (I did write a paper on that topic precisely a few years ago). What I found puzzeling is that persistence is large, extremely large. And the problem is that I do not see how we can explain ‘jumps’ that we do observe on that series. For instance the behavior of the series while I was in Europe, before January 20th: within 3 days, the temperature went down, from 0°C to -20°C, and up from -20°C to 0°C, and then down again, from 0°C to -20°C (a nice И if we use cyrillic letters). Or how can we explain the oscillating behavior observed the week after, where the temperature went up, from -25°C to (almost) +10°C in a few days. Within 10 days, we did observe also two ‘jumps’ (or ‘crashes‘ if we want to use the terminology of financial time series) with a decrease of 25 degrees in less than 24 hours ! Obviously, we need to find other classes of model to replicate that kind of behavior we observe on temperatures…

# Ruin probability and infinite time

A couple of weeks ago, I had a discussion with a practitioner, working in some financial company, about ruin, and infinite time. And it reminded me a weird result. Well, not a weird result, but a result I found disturbing, at first, when I was a student (that I rediscovered with the eyes of someone dealing with computational issues, seeing here a difficult theoretical question). Consider a simple ruin problem. A player has wealth . Then he flips a coin: tails he has a gain of 1, heads he experiences a loss of 1. At time , his wealth is where  is associated to the th coin:  is equal to 1 with probability (tails), and -1 with probability  (heads). It is also possible to write

where  can be interpreted as the net gain of the player. In order to get a good understanding of results that can be obtained. Assume  to be given. Let denote the number of heads and  the number of tails. Then , while . Let  denote the number of paths to go from point A (wealth  at time ) to point B (wealth  at time ). Note that this is a Markovian problem, that can be modeled using Markov chains

But here, we will focus on combinatorial results. Hence,

In order to derive probabilities to reach , let  denote the number of paths going from  to . And let denote the number of paths going from  to  that do reach  at some point between  and . Using a simple reflexion property, then if  and  are positive,

Based on those reflexions, two results can be derived (focusing on probability, instead of counting paths). First, we can obtain that

(given that n and x have the same parity). The second result we can obtain is that

Based on those two expressions, if  denotes the first time  become null, given ,

then

This can be computed easily,

> x=10
> p=.55
> ProbN=function(n){
+ pb=0
+ if(abs(n-x) %% 2 == 0)
+ pb=x/n*choose(n,(n+x)/2)*(1-p)^((n+x)/2)*(p)^((n-x)/2)
+ return(pb)}
> plot(Vectorize(ProbN)(1:1000),type="s")

That looks nice… But if we look closer, we can wonder what

would be ? Since we have the distribution of a probabilty measure, we might expect one. But here

> sum(Vectorize(ProbN)(1:1000))
[1] 0.134385

And this is not due to calculation mistakes that we do not get 1 here. Actually, we should write

which might be interpreted as the probability of ruin, starting from , that we denote  from now on. The term on the left can be approximated using monte-carlo simulations

> p=.55
> x=10
> m=1000
> simul=10000
> S=sample(c(-1,1),size=m*simul,replace=TRUE,prob=c(1-p,p))
> MS=matrix(S,simul,m)
> for(k in 2:m) MS[,k]=MS[,k]+MS[,k-1]
> T0=function(vm) which(vm<=(-x))[1]
> MTmin=apply(MS,1,T0)
> mean(is.na(MTmin)==FALSE)
[1] 0.1328

To check the validity of the relationship above, a simple (theoretical) recursive formula can be derived for the term on the right (ruin probability), namely

with a boundary conditions , and . Then is comes that

Note that it might be tricky to check using monte carlo simulation… since we cannot have an infinite number of runs. And we’re dealing precisely with things that do occur when time is infinite. Actually, we can still check convergence, considering an upper limit  for the number of runs, and then letting  go to infinity. Note that an explicit formula can then be derived (using additional border condition )

Using the following code, it is possible to calculate ruin probability, in order to estimate .

> MSmin=apply(MS,1,min)
> mean(MSmin<=(-x))
[1] 0.1328
> (((1-p)/p)^x-((1-p)/p)^m)/(1-((1-p)/p)^m)
[1] 0.1344306

The following graph shows the evolution of ruin probability as a function of initial wealth (with monte carlo simulation, with a fixed horizon – including a confidence interval – versus the analytical expression)

Hence, with stopping times, one should remember that

and that those two terms can be approximated simply using simulations or standard approximations.

# Halloween: zombies walk versus random walk

I don’t know if it’s because of Halloween (or the new season of walking dead), but I had weird thoughts recently. For instance, I wondered if zombies walk was a random walk… Note that I am not the only one: there are papers on that topic on Wired or even Nature (at least on one of their blogs). Before modeling that walk, it might be better to describe the process, and to generate zombies walk through simulations.

# Marches aléatoires en dimension 1, 2 ou 3, ou histoire éthylique

Revenons à nos histoires d’alcoolique un instant (les précédentes sont en ligne iciet ). L’avant dernière fois (ici), nous avions parlé de la marche aléatoire, et du temps mis avant de trouver la sortie de la rue de la soif.

Et j’avais évoqué la probabilité la probabilité qu’il ne récupère pas les clés (i.e. qu’il ne repasse pas par l’origine) s’il peut sortir un jour de la rue. En probabilité, on appelle cela la probabilité de repasser par l’origine et cela vérifie une propriété étrange dès lors que l’on change de dimension… On va suivre notre alcoolique et noter  sa position à la date , puis poser

autrement dit la probabilité de se retrouver à sa position initiale à la date n (et de retrouver ainsi ses clés), le  pouvait être en n’importe quelle dimension, on posera

correspondant à la probabilité de premier retour à la position initiale, et enfin

qui est la probabilité de revenir un jour au point de départ, c’est à dire … C’est cette dernière grandeur qui nous intéresse….

• De la probabilité de récupérer ses clés dans la rue

La rue est ici infinie (contrairement à l’autre fois, où on pouvait en sortir, soit vers la place de Lices, soit vers la place Sainte Anne). Autrement dit, on s’intéresse à la probabilité qu’une marche aléatoire sur  issue de  repasse un jour par . S’il va à droite ou à gauche avec probabilités  et  respectivement, notons que , car il faut un nombre de pas pairs pour revenir à l’origine, et

Il est possible de montrer que

et en particulier

On peut relier les  et les  en notant que

Or  constitue une distribution de probabilité, dont on peut définir la fonction génératrice

On notera que que l’on peut écrire

soit, en utilisant l’expression donnée auparavant

Soit finalement

ce qui se traduit par  si la série  converge. Aussi,

donc si , notre alcoolique repassera presque sûrement par le point d’origine, alors que si , il ne repasse pas par le bar où il a oublié ses clés avec probabilité |.

• De la probabilité de récupérer ses clés sur une place

La place est ici infinie de part et d’autre. Autrement dit, on s’intéresse à la probabilité qu’une marche aléatoire sur  issue de  repasse un jour par . Il a alors 1 chance sur 4 d’aller devant lui, autant d’aller en arrière, d’aller à gauche ou d’aller à droite. On peut noter que là encore
, car il faut un nombre de pas pairs pour revenir à l’origine, et

Or on peut réécrire cette somme en notant que

et que

aussi,

que l’on peut approcher (via la formule de Stirling) sous le forme

La série  est alors divergente, et donc . Aussi, même sur une place, notre alcoolique reviendra presque sûrement à sa place initiale (et ce une infinité de fois, car une fois qu’il y est revenu et qu’il en repart, il repassera presque sûrement à cette position, etc).

• De la probabilité de récupérer ses clés dans l’espace

Bon, cette fois notre alcoolique passe en dimension 3…. j’admets, il faut de l’imagination. On s’intéresse à la probabilité qu’une marche aléatoire sur  issue de  repasse un jour par . Il a alors 1 chance sur 6 d’aller devant lui, autant d’aller en arrière, d’aller à gauche ou d’aller à droite, d’aller en haut ou d’aller en bas.

Il y a désormais 6 voisins, et en raisonnant comme avant (le maillage est un peu plus complexe), on obtient que

Je laisse les amateurs de combinatoires faire des calculs ou des approximations, mais on peut montrer que

La série \sum u_n est alors convergente, et

qui vaut un peu plus de 1/3. Autrement dit, on a seulement 1 chance sur 3 de revenir à l’origine et de retrouver ses clés…
Bref, en dimension 1 ou 2, on repasse presque sûrement une infinité de fois à l’origine, mais en dimension 3, on a 2 chances sur 3 de ne jamais y repasser…. J’ai toujours trouvé ça étonnant….