Freakonometrics

To content | To menu | To search

Teaching › ACT6420-A2012

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.

Monday, December 3 2012

Méthodes de lissage en assurance

Dans le cours de modèles de prévision, j'avais abordé les méthodes de régression locale rapidement, en finissant la section sur les données individuelles. On verra plus d'outils la session prochaine dans le cours actuariat IARD, a.k.a. ACT2040 (que je donnerais cet hiver). En attendant, je mets en ligne les transparents du cours que vient de donner Julien Tomas à l'Institut de science financière et d'assurances (à Lyon, en France), dans un contexte d'assurance dépendance,

(et que Julien m'a autorisé à diffuser).

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)

Sunday, November 25 2012

My hobby

(via )

Wednesday, November 21 2012

Introduction aux modèles SARIMA

Quelques transparents en plus, qui devraient correspondre aux deux prochains cours de séries temporelles, sur les processus autorégressifs (AR) et moyennes mobiles (MA), les ARMA, les ARIMA (intégrés) et les SARIMA (saisonniers). J'ai mis des notes sur les tests de racine unités, je rajouterais quelques transparents la semaine prochaine sur les tests de saisonnalité, et quelques exemples pratiques de prévision. Les transparents sont en ligne ici,

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

Saturday, November 17 2012

Le mois des (boites à) moustaches

On célèbre cette année les 35 ans d'un livre qui a révolutionné la statistique, tant sur les aspects computationnels que graphiques: Exploratory Data Analysis par John Tukey. La légende prétend (on pourra relire a brief history of S par Richard Becker) que ce livre a inspiré ce qui allait devenir R, 
"Commercial software did not fit well into our research environment. It often used a ‘‘shotgun’’ approach — print out everything that might be relevant to the problem at hand because it could be several hours before the statistician could get another set of output. This was reasonable when computing was done in a batch mode. However, we wanted to be able to interact with our data, using Exploratory Data Analysis (Tukey, 1971) techniques. In addition, commercial statistical software usually didn’t compute what we wanted and was not set up to be modified."

En particulier, c'est dans ce livre qu'est introduit un outil graphique intéressant pour visualiser des données, le boxplot, appelé en français boite à moustaches (que tous les statisticiens se doivent d'utiliser en ce mois de Movembre). Dans la version originale, les boites sont verticales

mais mis horizontalement, on voit une boite, avec des moustaches. Cela dit, les moustaches ne sont pas toujours ce que l'on croit. Par exemple, on pense souvent (et moi le premier, il y a encore peu) que la boite à moustaches permettait de visualiser des valeurs extrêmes, voire des quantiles s'il y a trop de données (le maximum peut - en effet - être très très grand)

C'est malheureusement plus compliqué que ca... Afin de mieux comprendre ce qu'on trace, regardons le jeu de données suivant,
> set.seed(1)
> X=rlnorm(99)

(on est ici en échelle logarithmique). Pour commencer, on peut faire une boite à moustache avec une visualisation du minimum (à gauche) et du maximum (à droite)
> boxplot(X,horizontal=TRUE,log="x",range="0")

Et effectivement, on retrouve facilement tous les paramètres: ceux de la boite (en rouge)
> 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")
et ceux de la moustache (en bleu)
> abline(v=min(X),lty=2,col="blue")
> abline(v=max(X),lty=2,col="blue")
Mais ce n'est pas la version standard de la boite à moustaches, qui est - sous R - la suivante
> boxplot(X,horizontal=TRUE,log="x")

Si la borne inférieure de la moutache est toujours le minimum (en tous cas ici), pour la partie supérieure... c'est plus compliqué. Malheureusement ce n'est pas un quantile standard (ceux à 90% et 95% sont représentés ci-dessus en mauve). Pour comprendre ce que c'est, il faut regarder le code: ici, on rajoute au quartile supérieur (la partie de droite de la boite) un pourcentage de la distance inter-quartile (la longueur de la boite), et ce pourcentage est le suivant
> M=boxplot(X)$stats
> (M[5]-M[4])/(M[4]-M[2])
[1] 1.350628
> 2*qnorm(.75)
[1] 1.34898
Plusieurs sites ou ouvrages suggèrent d'utiliser 1.5 fois la longueur interquartile, mais sous R, on utilise . Il faut le savoir (je l'ai appris un peu par hasard en faisant une formation sur la régression quantile au R Montréal Group au printemps passé, cf. ici). Et effectivement, on retrouve les bornes de la moustache,
> abline(v=quantile(X,.75)+2*qnorm(3/4)*IQR(X),lty=2,col="blue")
> abline(v=min(X),lty=2,col="blue")
et je laisse les quantiles pour montrer que (malheureusement) ce n'est pas ce qui est utilisé ici,
> abline(v=quantile(X,.9),lty=2,col="purple")
> abline(v=quantile(X,.95),lty=2,col="purple")
Cela dit, comme je l'évoquais lors de la formation (ici), beaucoup de monde laisse planer une réelle ambiguité sur ce qui est représenté par ces boites à moustaches... Et le fait que ce ne soit pas un quantile me gene un peu... Par exemple, ici on avait la boite à moustaches suivante

Si je prends les 25 plus petites observations, et que je les divise par 4, on obtient

La partie de gauche bouge (c'est normal, c'est la partie qu'on a modifiée), et la médiane reste identique, ainsi que la partie droite de la boite (ce qui est normal, ce sont des quantiles au delà de 25%). Mais de manière un peu génante, la borne supérieure de la moustache s'est ralongée...  C'est normal, car on utilise une distance inter-quartile. Mais on n'a pas touché aux grandes observations, c'est donc génant de la voir bouger ainsi... non ?

Friday, November 16 2012

sur les transformations dans un modèle linéaire

Je voulais prendre 5 minutes pour reprendre une question posée par courriel, qui me permettra de poursuivre sur des choses évoquées mercredi dernier en cours. Je vais reformuler la question, mais en gros, cela disait: la méthode de Box-Cox (évoquée en début de semaine, ici) avait pour objet de choisir entre un modèle linéaire (que l'on étudié dans tous les sens depuis le début du cours) et un modèle log-linéaire (évoqué ici, par exemple). L'idée était d'associer le cas linéaire à (dans la transformée de Box-Cox) et à pour le cas multiplicatif (ou log-linéaire). Mais que se passe-t-il si la valeur optimale rejette ces deux cas, et est proche (disons) de ?
Considérons le cas suivant (exemple que je ne me lasse d'utiliser)

> reglm=lm(dist~speed,data=cars)
> library(MASS)
> boxcox(reglm)

La valeur optimale est effectivement proche de . En posant , on aurait envie de considérer un modèle de la forme
(comme le suggère la transformation de Box-Cox). On peut faire la régression,
> regsqrt=lm(sqrt(dist)~speed,data=cars)
+ summary(regsqrt)
 
Call:
lm(formula = sqrt(dist) ~ speed, data = cars)
 
Residuals:
Min      1Q  Median      3Q     Max
-2.0684 -0.6983 -0.1799  0.5909  3.1534
 
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  1.27705    0.48444   2.636   0.0113 *
speed        0.32241    0.02978  10.825 1.77e-14 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 1.102 on 48 degrees of freedom
Multiple R-squared: 0.7094,	Adjusted R-squared: 0.7034
F-statistic: 117.2 on 1 and 48 DF,  p-value: 1.773e-14
et effectivement, le résultat est concluant. Mais on ne peut pas en rester là, comme pour pour le passage au logarithme, notre variable d'intéret reste .
La prédiction sur le modèle transformé était
et
Comme avec le logarithme, on pourrait prendre comme prédiction pour  mais
Bref, on a un estimateur biaisé. Et là encore, comme on est sur des variables positives, l'inégalité de Jensen doit meme nous garantir que l'une des quantités domine toujours l'autre (je laisse les amateurs de convexité l'écrire dans le bon sens). Comment faire ?
La première solution est de noter que
Donc de manière très simple, on a notre prédiction (en prenant le carré, comme intuité, mais en rajoutant une terme - positif - lié à la variance).
Une seconde consiste à noter que Z (conditionnellement à la variable explicative) était gaussienne. Donc en prenant le carré (moyennant quelques changement d'échelle), on devrait tomber sur une loi du chi-deux qui est une loi que l'on connait bien (sinon je peux renvoyer ici).
Essayons de creuser un peu ces idées (en particulier, pour dériver des intervalles de confiance). Si on visualise la prédiction de ce modèle sur la racine carrée de la distance de freinage, on obtient
> plot(speed,sqrt(dist))
> x=seq(0,30,by=.2)
> distsqrtp=predict(regsqrt,newdata=
+ data.frame(speed=x),interval="prediction")
> polygon(c(x,rev(x)),c(distsqrtp[,2],
+ rev(distsqrtp[,3])),col="yellow",border=NA)
> lines(x,distsqrtp[,1],lwd=2,col="red")
Pour passer d'une loi normale à une loi du chi-deux, il faut retrancher la moyenne (pour centrer la variable) et diviser par l'écart-type (pour avoir une variance unitaire),
> s=summary(regsqrt)$sigma
> mu=predict(regsqrt)
> distsqrtp01=(sqrt(dist)-mu)/s
On a ainsi une variable qui, conditionnellement à la variable explicative est supposé suivre une loi ,
> plot(speed,distsqrtp01,ylim=c(-3,3))
> abline(h=qnorm(c(.025,.975),0,1),lty=2,col="red")
On prend ensuite le carré pour obtenir notre loi , on on trace la bande de confiance (il faut que la valeur à la deux soit faible),
> plot(speed,distsqrtp01^2)
> abline(h=qchisq(.95,df=1),lty=2,col="red")
On a maitenant nos intervalles de confiance, en utilisant cette loi du chi-deux, en inversant notre transformation,
> mu=predict(regsqrt,newdata=data.frame(speed=x))
> distsup=(mu+s*sqrt(qchisq(.95,df=1)))^2
> distinf=(mu-s*sqrt(qchisq(.95,df=1)))^2
que l'on peut visualiser simplement (sur les données de base, avec en ordonnée la distance)
> plot(cars)
> lines(x,distsup,lty=2,col="red")
> lines(x,distinf,lty=2,col="red")
Et la valeur prédite ? On va utiliser notre relation liant la variance, l'espérance du carré, et le carré de l'espérance,
> distesp=mu^2+s^2
> lines(x,distesp,lwd=2,col="red")
Nice, n'est-ce pas ?
Sinon, comme je le disais en cours, si la transformation optimale sur est de prendre , peut-etre pourrait-on envisager de prendre (de manière un peu duale) comme variable explicative. Encore une fois, c'est une idée, car rien ne le garantit. En effet
En particulier, le second modèle est un modèle linéaire classique, avec de l'homoscédasticité: on aura une dispersion uniforme autour de notre parabole. En revanche, avec le premier modèle, on a une espèce d'hétéroscédasticité (comme tenu du double produit), avec une variance du terme d'erreur qui croit avec (car les coefficients sont positifs). Bref, en terme d'intervalles de confiance, on devrait avoir des choses assez différentes.
Regardons la régression de la distance (cette fois) sur la vitesse, et le carré de la vitesse,
> reglm=lm(dist~speed+I(speed^2),data=cars)
> distp=predict(reglm,newdata=
+ data.frame(speed=x),interval="prediction")
Si on visualise la prédiction, avec un intervalle de confiance, on obtient
> plot(cars)
> polygon(c(x,rev(x)),c(distp[,2],
+ rev(distp[,3])),col="yellow",border=NA)
> lines(x,distp[,1],lwd=2,col="red")
Autrement dit, on suppose que notre modèle est homoscédastique. Les régions de confiance entre les deux approches sont clairement différentes... meme si les prédictions sont presque superposées (oui oui, il y a deux courbes, une rouge et une bleue).




Devoirs, ACT6420

Pour rappel, le premier devoir est à rendre pour dimanche soir (minuit - comme annoncé dans le plan de cours - date d'envoi du courriel faisant foi). Je veux un fichier pdf envoyé à mon adresse, avec sur la page de garde les codes étudiants. Pour les étudiants en groupes de deux, merci de mettre le coéquipier en copie du courriel. Je rappelle que des instructions sont toujours en ligne sur un vieux billet.

Histoire de clarifier les choses pour le second devoir, j'ai expliqué mercredi dernier les grandes lignes de ce qui serait demandé pour le devoir. Les instructions pour télécharger les données sont en ligne dans un précédant billet. Chaque groupe doit choisir un mot clé (de son choix, peu importe la langue et la localisation), extraire les données, et les transformer en données mensuelles. Les codes sont en ligne dans le  billet en question. Comme évoqué mercredi, compte tenu des configurations des comptes Google, et de R, il peut y avoir des soucis sur la manipulation des dates. Je serais disponible jusqu'au 1er décembre pour aider sur les aspects d'importation et de lecture des données, mais pas ensuite ! Pour le devoir, il faudra proposer une modélisation de la série, ainsi qu'une prédiction pour les 24 mois à venir, avec une intervalle de confiance. Pour rappel, ce second devoir sera à rendre pour le 23 décembre. Bon courage.

Wednesday, November 14 2012

Série temporelles

Pour la section sur les séries temporelles, je mets des notes de cours en ligne. Les slides seront en ligne d'ici la semaine prochaine. Sinon, quelques séries dans le fichier ci-dessous

> source(
+ "http://freakonometrics.blog.free.fr/public/data/sourcets.R")

Les données sont les suivantes, avec l'utilisation du mot clé épilation sur Google,

> ts.plot(epilation,col="red")

..

ou l'utilisation du mot clé gym,

> ts.plot(gym,col="blue")

On a aussi la série des décès suite à des accidents de la route en Ontario,

> ts.plot(ontario,col="purple")

ou la série du nombre de vente de voitures au Québec

> ts.plot(quebec,col="purple")

La série suivante est la température minimale journalière au parc Montsouris, à Paris,

> ts.plot(temperature,col="green")

Pour finir, deux séries de transport, avec le nombre de voitures qui ont emprunté l'autoroute A7 en France

> ts.plot(autoroute7,col="red")

ou le nombre de voyageur dans un aéroport,

> ts.plot(airline,col="blue")

- page 1 of 3