Freakonometrics

To content | To menu | To search

Tag - ACT6420

Entries feed - Comments feed

Thursday, December 20 2012

Examen ACT6420 sur les séries temporelles, correction

L'examen de mercredi est en ligne ici (les données et les annexes étaient dans un précédant billet), et je poste des éléments de correction. Une correction avec les statistiques de réponses sera en ligne dans les jours qui viennent (le temps de corriger). Si quelqu'un voit des coquilles dans la correction, voire des réponses erronées, merci de me le dire avant que je ne valide les notes.

Sinon pour le second devoir, le nombre de pages attendu est de l'ordre d'une dizaine, maximum quinze, avec les noms et les codes étudiants sur la page de garde. Je veux - par courriel - un fichier pdf avant le 23 décembre, 23:59 (date limite, mais il est possible de l'envoyer plus tôt).

Friday, December 14 2012

ACT6420 examen final

Mercredi prochain, c'est l'examen final (qui compte pour 30%).

Au programme, 33 questions à choix multiples, en 3 heures, environ la moitié sur la base suivante


> base=read.table( "http://freakonometrics.blog.free.fr/public/data/TS-examen.txt", + sep=";",header=TRUE) > X=ts(base$X,start=c(base$A[1],base$M[1]),frequency=12) > plot(X)

Les annexes qu'il faudra commenter lors de l'examen sont en ligne ici. Inutile de préciser que je ne répondrais à aucune question sur les sorties.

Friday, December 7 2012

Modélisation et prévision, cas d'école

Quelques lignes de code que l'on reprendra au prochain cours, avec une transformation en log, et une tendance linéaire. Considérons la recherche du mot clé headphones, au Canada, la base est en ligne sur http://freakonometrics.blog.free.fr/...

> report=read.table(
+ "report-headphones.csv",
+ skip=4,header=TRUE,sep=",",nrows=464)
> source("http://freakonometrics.blog.free.fr/public/code/H2M.R")
> headphones=H2M(report,lang="FR",type="ts")
> plot(headphones)

Mais le modèle linéaire ne devrait pas convenir, car la série explose,

> n=length(headphones)
> X1=seq(12,n,by=12)
> Y1=headphones[X1]
> points(time(headphones)[X1],Y1,pch=19,col="red")
> X2=seq(6,n,by=12)
> Y2=headphones[X2]
> points(time(headphones)[X2],Y2,pch=19,col="blue")

Il est alors naturel de prendre le logarithme de la série,

> plot(headphones,log="y")

C'est cette série que l'on va modéliser (mais c'est bien entendu la première série, au final, qu'il faudra prévoir). On commence par ôter la tendance (ici linéaire)

> X=as.numeric(headphones)
> Y=log(X)
> n=length(Y)
> T=1:n
> B=data.frame(Y,T)
> reg=lm(Y~T,data=B)
> plot(T,Y,type="l")
> lines(T,predict(reg),col="purple",lwd=2)

On travaille alors sur la série résiduelle. 

> Z=Y-predict(reg)
> acf(Z,lag=36,lwd=6)
> pacf(Z,lag=36,lwd=6)

On peut tenter de différencier de manière saisonnière,

> DZ=diff(Z,12)
> acf(DZ,lag=36,lwd=6)
> pacf(DZ,lag=36,lwd=6)

On ajuste alors un processus ARIMA, sur la série différenciée,

> mod=arima(DZ,order=c(1,0,0),
+ seasonal=list(order=c(1,0,0),period=12))
> mod
 
Coefficients:
ar1     sar1  intercept
0.7937  -0.3696     0.0032
s.e.  0.0626   0.1072     0.0245
 
sigma^2 estimated as 0.0046:  log likelihood = 119.47

Mais comme c'est la série de base qui nous intéresse, on utilise une écriture SARIMA,

> mod=arima(Z,order=c(1,0,0),
+ seasonal=list(order=c(1,1,0),period=12))

On fait alors la prévision de cette série.

> modpred=predict(mod,24)
> Zm=modpred$pred
> Zse=modpred$se

On utilise aussi le prolongement de la tendance linéaire,

> tendance=predict(reg,newdata=data.frame(T=n+(1:24)))

Pour revenir enfin à notre série initiale, on utilise les propriétés de la loi lognormales, et plus particulièrement la forme de la moyenne, pour prédire la valeur de la série,

> Ym=exp(Zm+tendance+Zse^2/2)

Graphiquement, on a 

> plot(1:n,X,xlim=c(1,n+24),type="l",ylim=c(10,90))
> lines(n+(1:24),Ym,lwd=2,col="blue")

Pour les intervalles de confiance, on peut utiliser les quantiles de la loi lognormale,

> Ysup975=qlnorm(.975,meanlog=Zm+tendance,sdlog=Zse)
> Yinf025=qlnorm(.025,meanlog=Zm+tendance,sdlog=Zse)
> Ysup9=qlnorm(.9,meanlog=Zm+tendance,sdlog=Zse)
> Yinf1=qlnorm(.1,meanlog=Zm+tendance,sdlog=Zse)
> polygon(c(n+(1:24),rev(n+(1:24))),
+ c(Ysup975,rev(Yinf025)),col="orange",border=NA)
> polygon(c(n+(1:24),rev(n+(1:24))),
+ c(Ysup9,rev(Yinf1)),col="yellow",border=NA)

Wednesday, December 5 2012

Syntaxe pour les SARIMA

La forme générale (enfin, la plus simple) des modèles SARIMA est la suivante

http://latex.codecogs.com/gif.latex?(1-L)^d(1-L^s)^{d_s}\Phi(L)\Phi_s(L^s)%20(X_t-m)%20=%20\Theta(L)\Theta_s(L^s)\varepsilon_t

Afin de comprendre comment les écrire sous R, commençons par un processus simple, autorégressive, de la forme
http://latex.codecogs.com/gif.latex?(1-\phi_1%20L)X_t%20=%20\varepsilon_t

La syntaxe est ici

> arima(X, order = c(p=1, d=0, q=0))

Supposons que l'on souhaite rajouter une composante moyenne mobile, i.e.

http://latex.codecogs.com/gif.latex?(1-\phi_1%20L)X_t%20=%20(1-\theta_1%20L)\varepsilon_t

La syntaxe devient

> arima(X, order = c(p=1, d=0, q=1))

Si on suppose maintenant que pour la composante autorégressive, on a une racine unité et le processus s'écrit

http://latex.codecogs.com/gif.latex?(1-%20L)X_t%20=%20(1-\theta_1%20L)\varepsilon_t

L’estimation de ce modèle se fait avec la commande

> arima(X, order = c(p=0, d=1, q=1))

Si pour finir, on veut calibrer un modèle de la forme

http://latex.codecogs.com/gif.latex?(1-%20L)(1-\phi_1%20L-\phi_2%20L^2)X_t%20=%20(1-\theta_1%20L)\varepsilon_t

on utilise la commande

> arima(X, order = c(p=2, d=1, q=1))

Supposons maintenant que notre série http://latex.codecogs.com/gif.latex?(X_t) ait été obtenu par différenciation saisonnière d'une série http://latex.codecogs.com/gif.latex?(Y_t), au sens où http://latex.codecogs.com/gif.latex?X_t=(1-L^{12})Y_t=\Delta^{12}%20Y_t, alors on devrait écrire

http://latex.codecogs.com/gif.latex?(1-%20L)(1-\phi_1%20L-\phi_2%20L^2)(1-L^{12})Y_t=%20(1-\theta_1%20L)\varepsilon_t

Pour estimer un tel modèle, la syntaxe est alors

> arima(X, order = c(p=2, d=1, q=1),
+ seasonal = list(order = c(0, 1, 0), period = 12))

En particulier, on peut avoir deux modèles assez proches pour modéliser une série cyclique: un bruit blanc avec une intégration saisonnière, i.e.

http://latex.codecogs.com/gif.latex?(1-L^{12})%20X_t%20=%20\varepsilon_t

dont la syntaxe serait

> arima(X, order = c(p=0, d=0, q=0),
+ seasonal = list(order = c(0, 1, 0), period = 12))

et un simple processus autorégressif à l'ordre 12. Mais là encore, soit je veux imposer que seule la dernière composante soit non nulle i.e.

http://latex.codecogs.com/gif.latex?(1-\phi_{12}%20L^{12})%20X_t%20=%20\varepsilon_t

ce qui s'écrit

> arima(X, order = c(p=0, d=0, q=0),
+ seasonal = list(order = c(1, 0, 0), period = 12))

soit je prends un processus plus général, de la forme

http://latex.codecogs.com/gif.latex?(1-\phi_1%20L-\cdots%20\phi_{12}%20L^{12})%20X_t%20=%20\varepsilon_t

dont la syntaxe devient

> arima(X, order = c(p=12, d=0, q=0))

En fait, on peut introduire un polynôme autorégressive avec uniquement des retards à l'ordre 12, en plus de la différentiation saisonnière, i.e.

http://latex.codecogs.com/gif.latex?(1-L)(1-L^{12})(1-\phi_1%20L)(1-\phi_{12}%20L^{12})%20X_t%20=%20(1-\theta_1%20L)\varepsilon_t

Ce modèle s'écrire alors sous la forme

> arima(X, order = c(p=1, d=1, q=1),
+ seasonal = list(order = c(1, 1, 0), period = 12))

On peut bien sûr aller plus loin, en autorisant non seulement une composante autorégressive saisonnière, mais pourquoi pas, une composante moyenne mobile saisonnière, i.e.

http://latex.codecogs.com/gif.latex?(1-L)(1-L^{12})(1-\phi_1%20L)(1-\phi_{12}%20L^{12})%20X_t%20=%20(1-\theta_1%20L)\varepsilon_t(1-\theta_{12}L^{12})\varepsilon_t

qui s'écrit, sous R,

> arima(X, order = c(p=1, d=1, q=1),
+ seasonal = list(order = c(1, 1, 1), period = 12))

On a vu la forme la plus générale des SARIMA. Enfin, comme je le disais en préambule, c'est aussi la plus simple. Car si je suppose qu'une autre cycle se superpose au cycle annuel, je peux le faire. Mais on va peut être essayer d'éviter de trop compliquer les notations...

Examen ACT6420 sur les séries temporelles

L'examen final du cours de modèles de prévision aura lieu dans 15 jours. Comme annoncé ce matin, la forme sera proche de celle de l'examen intra,

  • quelques questions de compréhension générales sur la modélisation des séries temporelles,
  • quelques questions portant sur de l'analyse de sorties obtenues suite à une modélisation d'une série.

Cette session, la série à étudier sera celle obtenue sur la fréquentation d'un aéroport, sur une quinzaine d'années. Les données sont mensuelles, et sont en ligne via le code suivant

> base=read.table(
"http://freakonometrics.blog.free.fr/public/data/TS-examen.txt",
+ sep=";",header=TRUE)
> X=ts(base$X,start=c(base$A[1],base$M[1]),frequency=12)
> plot(X)

Je mettrais les sorties en ligne dans 10 jours.

Friday, November 30 2012

Estimation et prévision pour des séries temporelles

Pour la fin de cours , quelques transparents sur l'identification et l'estimation des modèles SARIMA, quelques compléments sur les tests (racines unités et non-stationnarité, ainsi que saisonnalité), et enfin, quelques pistes pour construire des prédictions (avec une quantification de l'incertitude). Les transparents sont en ligne ici. Maintenant que les transparents sont finis (et en ligne) les prochains billets seront orientés autour de la modélisation et des aspects computationnels.

Thursday, November 29 2012

Save R objects, and other stuff

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)

Tuesday, November 20 2012

Transparents, introduction aux séries temporelles

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, . 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).

Tuesday, November 13 2012

Non, la majorité n'a pas toujours raison

J'ai fini de corriger le premier devoir, la correction est en ligne ici, avec des statistiques par réponses. On y retrouve en particulier que la majorité n'a pas toujours raison. Une question a d'ailleurs été supprimée (sur la transformation logarithmique, nous en parlerons demain en cours car c'est un point important, cf un précédant billet): l'examen est donc noté sur 29 points. J'ai passé aussi pas mal de temps à essayer d'expliquer comment se fait la mise à jour en enlevant une observation. On en reparlera peut-être aussi en cours. Sinon j'encourage tout le monde à aller jeter un œil à cette correction.

Mercredi, on fini les données individuelles, et on commencera les séries temporelles.

Monday, November 12 2012

Qui sort en avance lors des examens?

Tous les profs qui surveillent leur examen se posent la même question: que signifie un départ au bout d'une heure et demi, si l'examen est supposé durer trois heures? Que l'examen était trop facile? trop dur? 

Loin de moi l'idée de blâmer les étudiants qui partent tôt ! Tout d'abord parce que j'ai toujours souvent été dans les premiers à quitter la salle lors d'examens, ou de concours. Et parce que je ne pense pas qu'il y ait un lien évident entre le fait de partir tôt, ou tard, et la note obtenue. 

Mais on peut regarder (et répondre au commentaire de @kl suite à l'examen de la session passée). Car l'examen du cours ACT6420 avait lieu la semaine passée, et il était court. Autrement dit, 95% des étudiants étaient sortis avant que j'annonce la fin de l'épreuve. Et j'ai noté sur toutes les copies l'heure de sortie (pour ceux qui veulent plus d'information sur le protocole de l'expérience, les étudiants ne savaient pas que je notais l'heure).

Il s'agissait d'un examen à choix multiples, ce qui garantit une forme d'impartialité dans la correction (sinon on pourrait objecter que si je note en prenant le paquet de copie tel qu'il est arrivé, les notes pour les élèves partis tôt - corrigés en dernier - et ceux partis à la fin - corrigés en premier - pourraient ne pas être vraiment comparables). Et comme à la dernière session, je demandais aux étudiants de prédire leur nombre de réponses. La question était facultative, mais les élèves pouvaient gagner un point en répondant (à condition de prédire exactement le bon nombre).

Les données sont (partiellement) en ligne ici

> act6420=read.table(
+ "http://freakonometrics.blog.free.fr/public/
data/baseact6420a1.txt",
+ sep=";",header=FALSE,skip=1)
> names(act6420)=c("n","heure","obs","pred")

On peut alors comparer les notes en fonction de l'heure de sortie, ainsi que la note prédites (les notes sont sur 30, ou plutôt les nombres de bonnes réponses), en faisant des régressions linéaires simples,

> summary(lm(obs~heure,data=act6420))
 
Call:
lm(formula = obs ~ heure, data = act6420)
 
Residuals:
Min      1Q  Median      3Q     Max
-6.1122 -1.6326 -0.1876  1.9351  8.9224
 
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)   3.1519     6.0546   0.521   0.6039
heure         1.3209     0.5391   2.450   0.0162 *
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 2.819 on 93 degrees of freedom
Multiple R-squared: 0.06063,	Adjusted R-squared: 0.05053
F-statistic: 6.002 on 1 and 93 DF,  p-value: 0.01616
 
> summary(lm(pred~heure,data=act6420))
 
Call:
lm(formula = pred ~ heure, data = act6420)
 
Residuals:
Min       1Q   Median       3Q      Max
-11.9427  -2.5208   0.0052   3.2404   8.3733
 
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)   1.1178     8.8373   0.126   0.8996
heure         1.7234     0.7889   2.185   0.0317 *
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 3.999 on 85 degrees of freedom
(8 observations deleted due to missingness)
Multiple R-squared: 0.05316,	Adjusted R-squared: 0.04202
F-statistic: 4.773 on 1 and 85 DF,  p-value: 0.03167 

On observe qu'à la fois les bons élèves (ceux avec une bonne note) et les élèves qui se pensent bons (ceux qui ont prédit une note élevée) sortent plutôt tard. Avec dans les deux cas,  une significativité de la variable heure de sortie sur la prédiction de la note. Mais peut-être est-ce plus complexe que ça, et qu'il se cache un effet non-linéaire. On essayer de lisser un peu, avec des splines

> D=data.frame(Y=note,X=heure)
> library(splines)
> rg=lm(Y~bs(X,df=5),data=D)
> newheure=seq(10.2,12.1,by=.01)
> notep=predict(rg,newdata=data.frame(X=
+ newheure),interval="confidence")

On obtient les régressions (non-linéaires) suivantes

ou, si on souhaite rajouter des intervalles de confiance autour de la prédiction faite par nos modèles,

(avec en abscisse l'heure de sortie, l'examen finissant un peu après midi). L'analyse est un peu plus complexe qu'il n'y paraissait, avec quelques bons élèves qui finissent très tôt. Si on souhaite d'ailleurs mieux comprendre le lien entre la note obtenue (ou plutôt, encore une fois le nombre de bonnes réponses données, car aucun bonus n'est pris en compte ici), on peut faire un graphique proche de celui fait à la session passée

Le trait bleu est ici la régression obtenue à la session passée. On retrouve exactement la même chose sur les deux sessions, étonnant, non?

Mais j’arrête là pour l'analyse de ces notes. Pour ceux qui veulent en savoir plus, je mettrais bientôt en ligne une correction, avec des statistiques par réponse. On en profitera pour faire une discussion sur l'utilisation de règles démocratiques, à savoir, la majorité a-t-elle toujours raison? Et pour ceux qui n'auraient pas encore compris, on a fait ici une analyse de corrélations entre les heures de sortie, et les notes, et non pas une analyse de causalité. Donc inutile de se forcer à attendre la fin de l'examen pour sortir la prochaine fois, ça ne fera a pas (a priori) monter la note.

Mise en bouche pour l'analyse des séries temporelles

Mise en bouche n'est peut-être pas le terme adapté. Disons qu'avant d'attaquer la modélisation des séries temporelles, je vais prendre un peu de temps mercredi en fin de cours, pour rappeler quelques résultats sur les suites définies par récurrence. Considérons une suite réelle,  avec . On suppose qu'elle suit une relation de récurrence de la forme

(commencons par le degré 2, histoire de poser les choses). La stratégie pour trouver la forme (générale) d'un des termes de la suite est de passer par une forme vectorielle. On va alors considérer  définie par

alors

i.e. . Avec cette écriture, on note que l'on peut itérer .

Classiquement, pour résoudre cette équation, et donner la forme de la solution générale, on va passer par les valeurs propres de la matrice . Or en algèbre linéaire, on apprend qu'il suffit de passer par le polynôme caractéristique (cf http://fr.wikipedia.org/),

i.e. ici . Pour résoudre, ou utilise les formules classiques sur les polynômes de degré 2, en posant , et on en déduit la forme générales des racines. Supposons, pour commencer qu'il existe deux solutions, distinctes (dans  ou dans , peu importe) notées  et . On notera (par exemple) que

et donc, en multipliant par  on note que

Aussi, la suite   vérifie la relation de récurrence. Et pareil pour . On peut alors noter que toute combinaison linéaire vérifie la relation, i.e.

avec . Bref, on a trouvé un espace de solutions de dimension 2. Or quand on construit une suite par récurrence, cette dernière est entièrement caractérisée par  et . Autrement dit, l'espace des suites définies par la relation considérée est également de dimension 2. On a donc trouvé la forme de toutes les solutions de la relation de récurrence: toute suite vérifiant

est de la forme 

avec . On notera que l'on peut trouver les deux paramètres en fonction des deux premières valeurs de la suite.

Regardons maintenant un peu le cas où les racines sont dans . Comme on a un polynôme avec des coefficients dans , les solutions sont des nombres complexes conjugués, i.e. de la forme . La solution générale est alors maintenant de la forme

avec . On a ainsi des sinusoïdes amorties.

Pour ceux qui veulent jouer un peu, on peut tout programmer en R. On commence par se donner un polynôme (ou plutot les coefficients)

> PM=c(1,-2.5,4)
> polyroot(PM)
[1] 0.3125+0.3903124i 0.3125-0.3903124i

On peut ensuite construire une suite par récurrence (on va tirer au hasard les premières valeurs),

> set.seed(1); PM0=PM[2:length(PM)]
> x=runif(length(PM0)) > for(s in (length(PM0)+1):30){ + y= sum(-PM0*x[(s-1):(s-length(PM0))]) + x=c(x,y) +}

On peut faire un petit dessin,

> plot(x,cex=.2,type="b",xlab="",ylab="")
> points(1:length(PM0),x[1:length(PM0)],
+ pch=19,col="red",cex=.6)
> points((length(PM0)+1):30,x[(length(PM0)+1):30],
+ pch=19,col="blue",cex=.6)

On note que les racines sont à l'intérieur du disque unité, dans , la suite explose. On peut visualiser les racines du polynomes à l'aide des fonctions suivantes,

> plot(Re(polyroot(PM)),Im(polyroot(PM)))
> u=seq(-1,1,by=.01)
> polygon(c(u,rev(u)),c(sqrt(1-u^2),rev(-sqrt(1-u^2))),
+ col="yellow",border=NA)
> lines(u,sqrt(1-u^2),col="red",lwd=1.5)
> lines(u,-sqrt(1-u^2),col="red",lwd=1.5)
> points(Re(polyroot(PM)),Im(polyroot(PM)),
+ pch=19,col="blue")

On note que si des racines sont à l'extérieur du disque unité, dans , la suite explose. 

En revanche, si les racines sont dans le disque unité de , alors on observe une sinusoïde amortie (la longueur du cycle étant liée à l'argument des racines complexes conjuguées).

On va retrouver ces suites à plusieurs reprises lorsque l'on modélisera des séries temporelles. De manière très explicites avec les relations de Yule-Walker sur les autocorrélations. Mais aussi, en construisant les séries autorégressives comme des suites définies par récurrence, mais avec un bruit additionnel.

Wednesday, November 7 2012

Examen ACT6420 sur les modèles de régression

L'examen intra avait lieu ce matin. L'examen est en ligne ici, avec les annexes, en ligne . Des éléments de correction seront mis en ligne très bientôt.

Tuesday, October 30 2012

Créer sa propre base de données

Pour répondre à plusieurs questions, pour constituer sa propre base, il est important d'avoir différentes variables pour un même individu. Un individu peut être une entité plus large, comme une région, ou un quartier. Supposons que l'on s'intéresse aux votes lors d’élections. Désolé de reprendre des votes en France, mais les données sont faciles d'accès. Par exemple, pour Paris, on a les données de l'élection présidentielle de 2012, sur http://opendata.paris.fr/

> baseE=read.csv(
+ "http://freakonometrics.blog.free.fr/public/
data/election-paris-quartiers.csv",
+ header=TRUE,sep=";")

J'ai mes les données sur mon blog car c'est plus simple pour les lire. Les données ressemblent à ça,

> baseE[1:6,1:5]
ARR              QUARTIERS GRD_QUART INSCRITS EMARGEMENTS
1   1 SAINTGERMAINLAUXERROIS   7510101     1075         854
2   1                 HALLES   7510102     5667        4529
3   1            PALAISROYAL   7510103     1986        1610
4   1           PLACEVENDOME   7510104     1762        1427
5   2                GAILLON   7510205      799         613
6   2               VIVIENNE   7510206     1929        1485

Les données sont "par arrondissement et par quartier" comme on le lit dans le descriptif. Une ligne est un quartier. Maintenant, si on veut expliquer le taux de vote, ou le taux obtenu par tel ou tel candidat, il faut des variables explicatives. 

Pour cela, on peut aller sur http://www.recensement.insee.fr/ qui met en ligne beaucoup de données. En particulier, avec le même niveau de granularité, par quartier !

> base1=read.table(
+ "http://freakonometrics.blog.free.fr/public/
data/IRIS-PARIS-ACTIV.csv",
+ header=TRUE,sep=";")
Si on regarde dans la base,
> base1[1:15,c(1,2,3,4,5,7,8)]
IRIS REG DEP UU2010   COM TRIRIS GRD_QUART
1  751010101  11  75    851 75101 750011   7510101
2  751010102  11  75    851 75101 750011   7510101
3  751010103  11  75    851 75101 750011   7510101
4  751010104  11  75    851 75101 750011   7510101
5  751010105  11  75    851 75101 750011   7510101
6  751010199  11  75    851 75101 750011   7510101
7  751010201  11  75    851 75101 750021   7510102
8  751010202  11  75    851 75101 750021   7510102
9  751010203  11  75    851 75101 750021   7510102
10 751010204  11  75    851 75101 750021   7510102
11 751010205  11  75    851 75101 750021   7510102
12 751010206  11  75    851 75101 750021   7510102
13 751010301  11  75    851 75101 750011   7510103
14 751010302  11  75    851 75101 750011   7510103
15 751010303  11  75    851 75101 750011   7510103
En fait, le niveau d'agrégation est encore plus fin que les quartiers utilisés pour les élections. Mais on peut faire des sous-totaux, par quartier,
> base1t=aggregate(base1[,-(1:11)],
+ by=list(GRD_QUART=base1$GRD_QUART), FUN=sum)
> base1t[1:5,c(1,4,5,7,8)]
GRD_QUART P08_POP1564 P08_POP1524 P08_POP5564 P08_H1564
1   7510101        1278         222         250       662
2   7510102        7420        1187        1061      3987
3   7510103        2305         475         392      1171
4   7510104        2111         376         308      1010
5   7510205         875         111         136       426
On a enfin des bases construites sur les mêmes individus. Reste à fusionner,
> baseT=merge(baseE,base1t,by="GRD_QUART")
On a bien d'autres bases que cette base de population,
> base2=read.table(
+ "http://freakonometrics.blog.free.fr/public/
data/IRIS-PARIS-POPUL.csv",
+ header=TRUE,sep=";")
> base3=read.table(
+ "http://freakonometrics.blog.free.fr/public/
data/IRIS-PARIS-LOGEMENT.csv",
+ header=TRUE,sep=";")
> base4=read.table(
+ "http://freakonometrics.blog.free.fr/public/
data/IRIS-PARIS-FORMATION.csv",
+ header=TRUE,sep=";")
qu'il faut encore agréger par quartier,
> base2t=aggregate(base2[,-(1:11)],
+ by=list(GRD_QUART=base2$GRD_QUART), FUN=sum)
(etc) que l'on fusionne avec les autres bases
> baseT=merge(baseT,base2t,by="GRD_QUART")
On a alors enfin une base avec des données individuelles pour modéliser le taux de vote, par exemple,
> plot((baseT$P08_FNSCOL15P_SUP/baseT$P08_POP),
+ (baseT$VOTANTS/baseT$INSCRITS))
On notera qu'il est possible d'extraire une sous-base, avec les variables d’intérêt,
> base=data.frame(
+ Thollande=(baseT$HOLLANDE/baseT$INSCRITS),
+ Tsarkozy=(baseT$SARKOZY/baseT$INSCRITS),
+ T5pieces=(baseT$P08_RP_5PP/baseT$P08_RP),
+ Tfemmes=(baseT$P08_POPF/baseT$P08_POP),
+ Tpop2554=(baseT$P08_POP2554/baseT$P08_POP1564),
+ Tenfants018=((baseT$P08_POP0002+baseT$P08_POP0305
+ +baseT$P08_POP1117+baseT$P08_POP0610)/baseT$P08_POP),
+ Tpop5564=(baseT$P08_POP5564/baseT$P08_POP1564),
+ Tactifs=(baseT$P08_ACT5564/baseT$P08_POP1564),
+ Tchomeur=(baseT$P08_CHOM1564/baseT$P08_POP1564),
+ Tetudiant=(baseT$P08_ETUD1564/baseT$P08_POP1564),
+ T65ansplus=(baseT$P08_POP65P/baseT$P08_POP),
+ Npop=baseT$INSCRITS)
et c'est parti pour faire un modèle de régression...

Le passage au log dans les modèles linéaires

Un billet rapide pour compléter et illustrer le passage au log dans un modèle linéaire (que l'on abordera cette semaine en cours). Le point de départ est le modèle linéaire, où on suppose que, conditionnellement à  suit une loi normale. Pour rappel, si on a une loi normale, , alors  et . Les intervalles de confiance à 90% et 95% sont symétriques par rapport à la moyenne (qui est aussi la médiane, soit dit en passant),

Dans un modèle Gaussien avec homoscédasiticité,  i.e.  alors que . On a alors les bandes de confiance suivantes, pour un modèle de régression linéaire,

Bon, maintenant, que se passe-t-il si on prend l'exponentiel ? Pour la loi normale, rappelons que l'on obtient une loi lognormale, i.e. , les deux paramètres étant liés à la loi normale sous jacente, car désormais

alors que

Graphiquement, on a la loi suivante, avec les intervalles de confiance à 90% et 95% représentés ci-dessous. Le point noir est   alors que le point bleu est l'espérance de la loi lognormale.

On notera que le quantile de la loi log-normale est l'exponentiel du quantile de la loi normale. En effet, si http://latex.codecogs.com/gif.latex?\mathbb{P}(Y\leq%20q)=\alpha alors http://latex.codecogs.com/gif.latex?\mathbb{P}(\exp(Y)\leq%20\exp(q))=\alpha. En particulier,  n'est pas la moyenne de , mais la médiane (puisque  était la médiane de ).

Mais il n'est pas rare de voir utilisé un intervalle de confiance de la forme

qui est la forme classique de l'intervalle de confiance Gaussien (symétrique autour de la moyenne). Ici, on aurait les niveaux suivants

Notons qu'il n'y a aucune raison ici d'avoir une probabilité http://latex.codecogs.com/gif.latex?1-\alpha d'être dans l'intervalle de confiance obtenu avec les quantiles http://latex.codecogs.com/gif.latex?q_{1-\alpha/2} de la loi normale.

Maintenant, si on prend l'exponentiel d'un modèle linéaire (i.e. le logarithme de la variable d’intérêt est modélisé par un modèle linéaire) on a

avec une variance (conditionnelle) qui dépend de la variable explicative

Là encore, le plus naturel est d'utiliser comme bornes de l'intervalle de confiance des quantiles associés à la loi lognormale,

mais il n'est pas rare de voir utilisé des intervalles de type Gaussiens,

On perd là encore en interprétation car les bornes n'ont plus rien à voir avec les quantiles.

Wednesday, October 17 2012

ACT6420 examen (et mise en abyme)

L'examen intra du cours ACT6420 méthodes de prévisions de cet automne aura lieu (comme promis) le premier mercredi de novembre. La forme sera proche de l'examen donné cet hiver, avec 40 questions: une partie de questions "théoriques", portant sur la compréhension du cours, et une seconde partie portant sur l'analyse de sortie et de modélisation de données (avec uniquement des questions sur la régression sur données individuelles cette fois). Pour l'examen, la partie de modélisation portera sur la base suivante.

> base=read.table(
+ "http://freakonometrics.free.fr/examen-act6420-A.txt",
+ sep=";")

Il s'agit de données collectée lors de l'examen de cet hiver.

> tail(base)
NOTE PRED SEXE  AGE MAT3080
36 57.5 65.0    F 25.0    61.5
37 62.5   NA    F 21.4    95.0
38 50.0   NA    H 24.6    52.5
39 65.0 75.0    H 27.2    74.5
40 57.5 62.5    H 28.8    57.5
41 47.5 70.0    F 23.4    61.5

On dispose de 41 observations, avec la note obtenue à l'examen final (choix multiple, note sur 100), du sexe de l'étudiant ou de l'étudiante, de son age au moment de passer l'examen, de sa note (ou un proxy de sa note, la lettre étant ici convertie arbitrairement en chiffres) à l'examen de statistique MAT3080 (qui est un prérequis à ce cours), et une prédiction de la note, faite par chaque étudiant au moment de remettre sa copie (cette question était libre, et sans conséquence sur la note... même si un bonus était accordé à toute personne qui prédirait parfaitement son nombre de bonnes réponses). L'examen intra, des questions porteront sur la modélisation de la note des étudiants, voire de l'écart de perception entre la note obtenue, et la note espérée.

- page 1 of 4