Freakonometrics

Yesterday, Christopher asked me how to store an R object, to save some time, when working on the project. 

First, download the csv file for searches related to some keyword, via http://www.google.com/trends/, for instance "sunglasses". Recall that csv files store tabular data (numbers and text) in plain-text form, with comma-separated values (where csv term comes from). Even if it is not a single and well-defined format, it is a text file, not an excel file!

Recherche sur le Web : intérêt pour sunglasses
Dans tous les pays; De 2004 à ce jour
 
Intérêt dans le temps
Semaine,sunglasses
2004-01-04 - 2004-01-10,48
2004-01-11 - 2004-01-17,47
2004-01-18 - 2004-01-24,51
2004-01-25 - 2004-01-31,50
2004-02-01 - 2004-02-07,52

(etc) The file can be downloaded from the blog,

> report=read.table(
+ "http://freakonometrics.blog.free.fr/public/data/glasses.csv",
+ skip=4,header=TRUE,sep=",",nrows=464)

Then, we just run a function to transform that data frame. It can be found from a source file

> source("http://freakonometrics.blog.free.fr/public/code/H2M.R")

In this source file, there is function that transforms a weekly series into a monthly one. The output is either a time series, or a numeric vector,

> sunglasses=H2M(report,lang="FR",type="ts")

Here, we asked for a time series,

> sunglasses
Jan      Feb      Mar      Apr
2004 49.00000 54.27586 66.38710 80.10000
2005 48.45161 58.25000 69.93548 80.06667
2006 49.70968 57.21429 67.41935 82.10000
2007 47.32258 55.92857 70.87097 84.36667
2008 47.19355 54.20690 64.03226 79.36667
2009 45.16129 50.75000 63.58065 76.90000
2010 32.67742 44.35714 58.19355 70.00000
2011 44.38710 49.75000 59.16129 71.60000
2012 43.64516 48.75862 64.06452 70.13333
May      Jun      Jul      Aug
2004 83.77419 89.10000 84.67742 73.51613
2005 83.06452 91.36667 89.16129 76.32258
2006 86.00000 92.90000 93.00000 72.29032
2007 86.83871 88.63333 84.61290 72.93548
2008 80.70968 80.30000 78.29032 64.58065
2009 77.93548 70.40000 62.22581 51.58065
2010 71.06452 73.66667 76.90323 61.77419
2011 74.00000 79.66667 79.12903 66.29032
2012 79.74194 82.90000 79.96774 69.80645
Sep      Oct      Nov      Dec
2004 56.20000 46.25806 44.63333 53.96774
2005 56.53333 47.54839 47.60000 54.38710
2006 51.23333 46.70968 45.23333 54.22581
2007 56.33333 46.38710 44.40000 51.12903
2008 51.50000 44.61290 40.93333 47.74194
2009 37.90000 30.38710 28.43333 31.67742
2010 50.16667 46.54839 42.36667 45.90323
2011 52.23333 45.32258 42.60000 47.35484
2012 54.03333 46.09677 43.45833  

that we can plot using

> plot(sunglasses)

Now we would like to store this time series. This is done easily using

> save(sunglasses,file="sunglasses.RData")

Next time we open R, we just have to use

> load("/Users/UQAM/sunglasses.Rdata")

to load the time series in R memory. So saving objects is not difficult.

Last be not least, for the part on seasonal models, we will be using some functions from an old package. Unfortunately, on the CRAN website, we see that

but nicely, files can still be found on some archive page. On Linux, one can easily install the package using (in R)

> install.packages(
+ "/Users/UQAM/uroot_1.4-1.tar.gz",
+ type="source")

With a Mac, it is slightly more complicated (see e.g. Jon's blog): one has to open a Terminal and to type

R CMD INSTALL /Users/UQAM/uroot_1.4-1.tar.gz

On Windows, it is possible to install a package from a zipped file: one has to download the file from archive page, and then to spot it from R.

The package is now installed, we just have to load it to play with it, and use functions it contains to tests for cycles and seasonal behavior,

> library(uroot)
> CH.test
function (wts, frec = NULL, f0 = 1, DetTr = FALSE, ltrunc = NULL)
{
s <- frequency(wts)
t0 <- start(wts)
N <- length(wts)
if (class(frec) == "NULL")
frec <- rep(1, s/2)
if (class(ltrunc) == "NULL")
ltrunc <- round(s * (N/100)^0.25)
R1 <- SeasDummy(wts, "trg")
VFEalg <- SeasDummy(wts, "alg")

(etc)

Last week, I went to the Econometrics seminar of Montréal, at UdM, where Alfred Galichon was giving a great talk on marriage market. Alfred is a former colleague (from France), a co-author, an amazing researcher, and above all, a friend of mine. And he has always be supportive about my blogging activities. So while we were having lunch, after the seminar, Alfred mentioned my blogging activity to the other researchers. I should say researchers in Econometrics (yes, with a capital E, since it is a Science, as mention in an old paper by David Hendry by the end of the 70's). Usually, when I am involved in this kind of meeting, I start with some apologies, explaining that I do like theoretical econometrics (if not, I would not come to the seminar), but I do like my freakonometrics activity. I do like to use econometrics (or statistical techniques) to figure out (at least to try) why some things works the way they do. I try to find data, and then try to briefly analyze them to answer some simple questions. Or sometime, I just run simulations to answer more theoretical questions (or at least to give clues).

But above all, I like the fact that blogging gives me the opportunity to interact with people I would never meet without this activity. For instance, last May, I was discussing (on Twitter) with @coulmont, @joelgombin and @imparibus about elections in France. Then @coulmont asked me "yes, everyone knows that there should be some ecological fallacies behind my interpretation, but I am not so sure since I have data with a small granularity. As an econometrician, what do you think ?" Usually, I hate having a label, like "... I ask you since you're a mathematician", or "as an economist, what do you think of...". Usually, when people ask me economic questions, I just claim being a mathematician, and vice-versa. But here, I even put on the front of my blog the word "econometrics" (more or less). So here, I could not escape... And the truth is, that while I was a student, I never heard anything about this "ecological fallacy". Neither did I as a researcher (even if I have been reading hundreds of econometric articles, theoretical and applied). Neither did I as a professor (even if I have been teaching econometrics for almost ten years, and I have read dozens textbooks to write notes and handouts). How comes ? How come researchers in sociology and in political sciences know things in econometrics that I have never heard about ?

The main reason - from my understanding - is the following: if everyone talks about "interdisciplinarity" no one (perhaps a few) is really willing to pay the price of working on different (not to say many) areas. I tried, and trust me, I found it difficult. It is difficult to publish a paper in a climate journal when you're not specialist in climate (and you just want to give your opinion as a statistician). It is difficult to assume that you might waste weeks (not to say months) reading articles in geophysics if you want to know more about earthquakes risks, going to seminars, etc. Research is clearly a club ("club" as defined in Buchanan (1965)) story.

This week, I planned to go to some journal club in biology and physics, at McGill (kindly, a colleague there invited me, but we got a time misunderstanding)... this has nothing to do with my teaching, nor with my research activities. But I might learn something ! Yes, I do claim that I am paid just to have fun, to read stuff that I do find interesting, trying to understand the details of a proof, trying to understand how data were obtained. In most cases, it might (and should) be a complete waste of time, since I will not publish anything (anything serious, published in some peer reviewed journal) on that topic... but should I really care ? As I explained earlier (in French), I do also claim that I have a moral obligation to return everything I have seen, heard, read. And since I am not a big fan of lectures (and that I do not think I have skills for that) I cannot give my opinion, neither on economics facts (as @adelaigue or @obouba might do on their blogs) or on science results (as @tomroud does). But I think I can help on modeling and computational issues. My point being: never trust what you read (even on my blog) but please, try to do it yourself! You read that "90% of French executive think about expatriation" (as mentioned here)? Then try to find some data that should confront that statement. And see if you come up with the same conclusion... And since it might be a bit technical sometimes, here are some lines of code, to do it on your own... Academics have a legitimacy when they give their opinions on technical issues. At least they can provide with a list of reference everyone should read to get an overview of the topic. But academics can also help people read graphs, or data. To give them "numeracy" (or a culture in numbers) necessary to understand complex issues.

To conclude, I should mention that I understood what this "ecological fallacy" was from  Thomsen (1987) and many more documents could be found on Soren Thomsen's page http://www.mit.ps.au.dk/srt/. But I got most of the information I was looking for from a great statistician, who happens to be also an amazing blogger: Andrew Gelman (see http://andrewgelman.com/). I will probably write a post someday about this, since I found the question extremely interesting, and important.

(via http:/xkcd.com/605/ ) )))

Depuis que le blog (ou le précédent) existe, je sais (par les messages que je reçois par courriel) qu’il est beaucoup lu en Afrique (francophone). J’ai reçu avant hier un livre publié par Aymric Kamega et Frédéric Planchet, tiré de la thèse de doctorat d’Aymric (en ligne sur http://halshs.archives-ouvertes.fr/) défendue en décembre dernier. J’étais alors rapporteur extérieur, et j’avais souligné l’intérêt des travaux d’Aymric dans le cadre de la volonté de la CIMA (autorité de contrôle régionale des marchés d’assurance pour l’Afrique subsaharienne francophone) de fournir aux assureurs de la région des outils adaptés. En particulier des tables de mortalité d’expérience, propres à la région. En effet, (comme le rappelle Aymric), suite aux états généraux (de l’assurance vie)  en 2007, le principe de la construction de tables de mortalité d’expérience a été adopté en remplacement des tables de mortalité de la population générale française entre 1960 et 1964 (tables dites PM 60-64 et PF 60-64) jusqu’alors imposées. Aymric avait travaillé sur ce sujet, et a soutenu une thèse de doctorat qui avait fait l'unanimité. J’avais pris beaucoup de plaisir à lire la thèse, j’en pense que j’en aurais à lire le livre (disponible sur le site de l’éditeur http://seddita.com/).

Nous allons commencer demain la modélisation des séries temporelles. Les transparents de la séance sont en ligne ici. Je rappelle que les notes de cours (complètes) sont également en ligne, là. Je continuerai à poster régulièrement des billets contenant des commandes R.

Au programme cette semaine, l'utilisation des méthodes de régression pour extraire tendance et cycle, le lissage exponentielle (simple, double et saisonnier), et une présentation des notions importantes dans le cours (stationnarité, autocorrélations, bruit blanc, etc).

A couple of days ago, Jeffrey sent me an e-mail, encouraging me to write about the international year of statistics, hoping that I would participate to the event. The goals of Statistics2013 are (as far as I understood) to increase public awareness of the power and impact of statistics on all aspects of society, nuture Statistics as a profession and promote creativity and development in Statistics.

So I will try to write more posts (yes, I surely can try) on statistical methodology, with sexy applications... So see you in a few days on this blog, to celebrate statistics ! I can even invite guests to write here... all contributors are welcome !

> set.seed(1)
> X=rlnorm(99)

> boxplot(X,horizontal=TRUE,log="x",range="0")

> abline(v=median(X),lty=2,col="red")
> abline(v=quantile(X,.25),lty=2,col="red")
> abline(v=quantile(X,.75),lty=2,col="red")

> abline(v=min(X),lty=2,col="blue")
> abline(v=max(X),lty=2,col="blue")

> boxplot(X,horizontal=TRUE,log="x")

> M=boxplot(X)$stats
> (M[5]-M[4])/(M[4]-M[2])
[1] 1.350628
> 2*qnorm(.75)
[1] 1.34898

> abline(v=quantile(X,.75)+2*qnorm(3/4)*IQR(X),lty=2,col="blue")
> abline(v=min(X),lty=2,col="blue")

> abline(v=quantile(X,.9),lty=2,col="purple")
> abline(v=quantile(X,.95),lty=2,col="purple")