Freakonometrics

To content | To menu | To search

Tag - ARIMA

Entries feed - Comments feed

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

Friday, September 21 2012

Transformation logarithmique de séries temporelles

Pour poursuivre une discussion amorcée en fin de cours, dans certains cas, on peut avoir l'impression que modéliser une série pourrait être compliqué,

plot(X,xlim=c(1,length(X)+20))

Mais on peut avoir l'intuition que modéliser le logarithme de la série pourrait être plus simple,

> X=log(Y)
> plot(X,xlim=c(1,length(X)+20))

On va alors tenter une modélisation par un processus ARMA de cette dernière série,

> md=arima(X,c(12,0,1))
> P=predict(md,24)
> E=P$pred
> V=P$se^2

On peut alors faire une prévision sur cette série plus simple à modéliser, et visualiser cette prévision.

> temps=length(X)+1:24
> ciu=(E+2*sqrt(V))
> cil=(E-2*sqrt(V))
> polygon(c(temps,rev(temps)),c(ciu,
+ rev(cil)),col="yellow",border=NA)
> lines(temps,E,col="red",lwd=2)
> lines(temps,ciu,col="red",lty=2)
> lines(temps,cil,col="red",lty=2)

Maintenant, on va devoir remonter. On va utiliser un résultat que l'on a vu sur la transformation logarithmique dans une régression: si après avoir pris le logarithme, on a un modèle simple, Gaussien, c'est que le modèle initial était log-normal. On peut alors utiliser les propriétés de la loi lognormale, dont on connait les moments à partir de ceux de la loi Gaussienne sous-jacente. Pour la prévision, on n'a pas trop le choix,

> mu=exp(E+.5*V)

Par contre, pour construire un intervalle de confiance, soit on utilise la variance de notre loi lognormale pour avoir la variance de notre processus, et on oublie cette histoire de loi lognormale pour construire un intervalle Gaussien,

> sig2=(exp(V)-1)*exp(2*E+V)
> ci1u=mu+2*sqrt(sig2)
> ci1l=mu-2*sqrt(sig2)

ou alors on utilise le fait que comme la transformation est monotone, l'intervalle de confiance peut etre vu comme une transformation du précédant intervalle de confiance,

> ci2u=exp(E+2*sqrt(V))
> ci2l=exp(E-2*sqrt(V))

Si on compare visuellement les deux, on a dans le premier cas,

> plot(Y,xlim=c(1,length(X)+20))
> temps=length(X)+1:24
> polygon(c(temps,rev(temps)),c(ci1u,
> rev(ci1l)),col="yellow",border=NA)
> lines(temps,mu,col="red",lwd=2)
> lines(temps,ci1u,col="red",lty=2)
> lines(temps,ci1l,col="red",lty=2)

(qui est symétrique et centré sur notre prévision) et dans le second

> plot(Y,xlim=c(1,length(X)+20))
> temps=length(X)+1:24
> ci1u=mu+2*sqrt(sig2)
> ci1l=mu-2*sqrt(sig2)
> polygon(c(temps,rev(temps)),c(ci2u,
> rev(ci2l)),col="yellow",border=NA)
> lines(temps,mu,col="red",lwd=2)
> lines(temps,ci2u,col="red",lty=2)
> lines(temps,ci2l,col="red",lty=2)

modélisation des données du Nil, et prévision

Un billet pour remettre le code de ce soir, sur la modélisation du Nil (sans le détail de tous les graphiques). Tout d'abord, dans un premier temps, on notait qu'il devait y avoir une rupture dans la série, et on optait pour un modèle basé sur les dernières observations (on enlevait les 30 premières),

seuil=30
library(datasets)
X=as.numeric(Nile)
Xc=X[-(1:seuil)]
plot(1:length(X)-seuil,X,xlim=c(-seuil,90),type="l")
lines(Xc,lwd=3)

On va alors commencer à la date 0 à partir de la trentième observation. On va translater la série pour avoir une série centrée sur 0 (dans un premier temps).

PP.test(Xc)
Y=Xc-mean(Xc)
plot(acf(Y))	

On regardant un peu la série, on rejette l'hypothèse de bruit blanc, et de cycle. Donc on se lance dans un modélisation par un processus http://latex.codecogs.com/gif.latex?ARMA(p,q).

pacf(Y,30)
model1=arima(Y,order=c(4,0,0))
model1
acf(residuals(model1))
Box.test(residuals(model1),lag=12,type='Box-Pierce')

On peut commencer par valider une modélisation http://latex.codecogs.com/gif.latex?AR(4). Ensuite,

acf(Y,30)
model2=arima(Y,order=c(0,0,4))
model2
acf(residuals(model2))
Box.test(residuals(model2),lag=12,type='Box-Pierce')

on va valider une modélisation http://latex.codecogs.com/gif.latex?MA(4). Enfin, on peut voir parmi la sélection automatique de modèle,

library(caschrono)
armaselect(Y,nbmod=5)

Le première modèle que l'on va tester, un processus http://latex.codecogs.com/gif.latex?MA(1) est acceptée.

model3=arima(Y,order=c(0,0,1))
model3
acf(residuals(model3))
Box.test(residuals(model3),lag=6,type='Box-Pierce')
 
BP=function(h) Box.test(residuals(model3),lag=h,
type='Box-Pierce')$p.value
plot(1:24,Vectorize(BP)(1:24),type='b',
ylim=c(0,1),col="red")
abline(h=.05,lty=2,col="blue")

En revanche les autres modèles testés sont rejetés.

model4=arima(Y,order=c(1,0,1))
model4
model5=arima(Y,order=c(2,0,1))
model5

On a trois modèles, et c'est déjà pas mal. Faisons un peu de backtesting pour comparer.

AIC(model1)
AIC(model2)
AIC(model3)
b=10
T=length(Y)
sY=Y[-((T-b):T)]
smodel1=arima(sY,order=c(4,0,0))
smodel2=arima(sY,order=c(0,0,4))
smodel3=arima(sY,order=c(0,0,1))
Yp1=forecast(model1,b+1)
Yp2=forecast(model2,b+1)
Yp3=forecast(model3,b+1)
M=cbind(Y[((T-b):T)],Yp1$mean,Yp2$mean,Yp3$mean)
sum((M[,1]-M[,2])^2)
sum((M[,1]-M[,3])^2)
sum((M[,1]-M[,4])^2)

Le premier modèle semble plus intéressant. On va donc le retenir pour faire de la prévision.

plot(1:length(X)-seuil,X,xlim=c(-seuil,90),type="l")
abline(h=mean(Xc),col="green")
lines(Xc,lwd=3)
polygon(c(temps,rev(temps)),
c(Yhat$lower[,2],
rev(Yhat$upper[,2]))+mean(Xc),col="yellow",border=NA)
polygon(c(temps,rev(temps)),c(Yhat$lower[,1],
rev(Yhat$upper[,1]))+mean(Xc),col="orange",border=NA)
lines(temps,Yhat$mean+mean(Xc),col="red")

Maintenant, on peut se demander s'il n'y avait pas une tendance linéaire (on va faire comme si, pour l'exercice)

plot(1:length(X)-seuil,X,xlim=c(-seuil,90),type="l")
lines(Xc,lwd=3)
temps=1:length(Xc)-seuil
abline(lm(Xc~temps),lwd=2,col="red")

Dans ce cas, on modèlise la série à laquelle on a retranché cette tendance linéaire. On reprend tous les étapes précédentes et cette fois (pour changer un peu) on va faire comme si on retenait un modèle moyenne-mobile.

Y=Xc-predict(lm(Xc~temps))
plot(Y,type="l")
pacf(Y,30)
acf(Y,30)
model2=arima(Y,order=c(0,0,4))
model2
acf(residuals(model2))
Box.test(residuals(model2),lag=12,type='Box-Pierce')

Pour la prévision, on va projeter la série stationnaire et tout remonter, en tenant compte de la tendance,

Xc=X[-(1:seuil)]
tps=1:length(Xc)
base=data.frame(Xc,tps)
tendance=lm(Xc~tps,data=base)
T=length(Y)
temps=T+1:20
pred=predict(tendance,newdata=
data.frame(tps=c(temps,rev(temps))))
plot(1:length(X)-seuil,X,xlim=c(-seuil,90),type="l")
lines(Xc,lwd=3)
abline(tendance,col="green")
polygon(c(temps,rev(temps)),c(Yhat$lower[,2],
rev(Yhat$upper[,2]))+pred,col="yellow",border=NA)
polygon(c(temps,rev(temps)),c(Yhat$lower[,1],
rev(Yhat$upper[,1]))+pred,col="orange",border=NA)
lines(temps,Yhat$mean+pred[1:length(temps)],col="red")

L'étape suivante sera du faire du backtesting avec ces deux types de modèles, avec constante versus avec tendance... Mais je ne vais pas tout faire non plus...

Wednesday, September 19 2012

Prévision, choix de modèle(s) et backtesting

Un billet pour reprendre (et compléter) le code vu hier, en cours de modèles de prévision. On a toujours pour but de prévoir l'évolution du trafic autoroutier, en France, sur l'autoroute A7.

> autoroute=read.table(
+ "http://freakonometrics.blog.free.fr/public/data/autoroute.csv",
+ header=TRUE,sep=";")
> a7=autoroute$a007
> A7=ts(a7,start = c(1989, 9), frequency = 12)

Comme je l'ai dit en cours, les constantes dans les modèles ARIMA, a complexifie inutilement les notations, donc pour faire des choses proches de ce qu'ai pu faire dans le cours, on va centrer la série, et considérer http://latex.codecogs.com/gif.latex?X_t-\overline%20X

> A7d=A7-mean(A7)
> plot(A7d,xlim=c(1989,2000))

On avait vu qu'un modèle http://latex.codecogs.com/gif.latex?AR(12) pouvait convenir, donc on va l'utiliser pour faire de la prévision.

>  fit=arima(A7d,order=c(p=12,d=0,q=0),
+      include.mean =FALSE)
>  summary(fit)
Series: A7d
ARIMA(12,0,0) with zero mean
 
Coefficients:
ar1      ar2      ar3     ar4      ar5     ar6
      0.0981  -0.0353  -0.0375  0.0530  -0.0934  0.0598
s.e.  0.0630   0.0623   0.0600  0.0606   0.0618  0.0568
ar7     ar8     ar9     ar10    ar11    ar12
      -0.0839  0.0278  0.0072  -0.1080  0.1578  0.7887
s.e.   0.0619  0.0629  0.0627   0.0629  0.0637  0.0623
 
sigma^2 estimated as 12823951:  log likelihood=-827.24
AIC=1678.48   AICc=1683.6   BIC=1710.23
 
In-sample error measures:
ME         RMSE          MAE          MPE
257.8363135 3581.0544202 2335.0686684    1.7258390
MAPE         MASE
34.7942409    0.2603139 

La fonction permettant de faire de la prévision est alors

> horizon=48
>
library(forecast) > A7hat=forecast(fit,h=horizon)

On a alors un ensemble de séries temporelles, la première étant la prévision proprement dite http://latex.codecogs.com/gif.latex?{}_T\widehat{X}_{T+h}, i.e. http://latex.codecogs.com/gif.latex?EL(X_{T+h}|\underline{X}_T), où http://latex.codecogs.com/gif.latex?EL(\cdot) désigne une espérance linéaire (ou en terme technique, une projection sur l'adhérence de l'ensemble des des combinaisons linéaires du passé). On a ensuite les intervalles de confiance estimés, à partir d'une hypothèse de normalité, en utilise un estimateur de http://latex.codecogs.com/gif.latex?{}_T\text{mse}_{T+h}=\mathbb{E}([{}_T\widehat{X}_{T+h}-X_{T+h}]^2|\underline{X}_T),

> A7hat
Point    Lo 80   Hi 80    Lo 95    Hi 95
Oct 1996    -7276.44 -11865.7 -2687.1 -14295.1  -257.70
Nov 1996   -10453.79 -15065.1 -5842.4 -17506.1 -3401.39
Dec 1996    -6970.67 -11583.5 -2357.8 -14025.3    84.02
Jan 1997   -13404.09 -18021.2 -8786.9 -20465.4 -6342.78

On peut alors représenter ces valeurs graphiquement, comme on l'a fait en cours

>  (TPS=tsp(A7hat$mean))
[1] 1996.750 2000.667   12.000
>  lines(A7hat$mean,col="red",lwd=2)
>  temps=seq(from=TPS[1],to=TPS[2],by=1/TPS[3])
>  lines(temps,A7hat$lower[,1],col="blue",lwd=2)
>  lines(temps,A7hat$lower[,2],col="purple",lwd=2)
>  lines(temps,A7hat$upper[,1],col="blue",lwd=2)
>  lines(temps,A7hat$upper[,2],col="purple",lwd=2)
>  abline(v=TPS[1]-.5/TPS[3],col="green")

En allant un peu plus loin, on peut faire plus joli, en utilisant des aires colorées afin d'illustrer une région de confiance

> plot(A7d,xlim=c(1989,2000))
> polygon(c(temps,rev(temps)),c(A7hat$lower[,2],
+ rev(A7hat$upper[,2])),col="yellow",border=NA)
> polygon(c(temps,rev(temps)),c(A7hat$lower[,1],
+ rev(A7hat$upper[,1])),col="orange",border=NA)
> lines(temps,A7hat$mean,col="red")
> abline(v=TPS[1]-.5/TPS[3],col="green")

Maintenant, si on revient sur notre série initiale, l'autocorrélogramme pourrait laisser penser qu'il y a une racine unité saisonnière,

> plot(acf(A7d,lag=36))

Donc on peut tenter de différencier à l'ordre 12,

> plot(acf(diff(A7d,12),lag=36))

On pourrait avoir l'impression que la série ainsi différenciée est un bruit blanc. On peut faire un test de Box-Pierce pour tester cette hypothèse,

> BP=function(h) Box.test(diff(A7d,12),lag=h,
+ type='Box-Pierce')$p.value
> plot(1:24,Vectorize(BP)(1:24),type='b',
+ ylim=c(0,1),col="red")
> abline(h=.05,lty=2,col="blue")

On pourrait valider ici l'hypothèse de bruit blanc. On va donc ajuster un (pur) modèle ARIMA saisonnier,

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

Si on compare les deux modèles,

> AIC(fit1)
[1] 1511.815
> AIC(fit2)
[1] 1245.374

Le second modèle ayant un critère d'Akaike plus faible, on aurait tendance  à préférer ce modèle au premier. Mais cette méthode de comparaison n'est pas forcément la plus pertinent, si l'objet est de faire une prévision à court terme (disons quelques mois). Pour cela, on va plutôt enlever les dernières observations, ajuster un modèle sur la série restreinte, et comparer les prévisions, avec ce qui a été réellement observé. Disons qu'on va faire du backtesting sur les 9 derniers mois,

> T=length(A7d)
> backtest=9
> subA7d=ts(A7d[-((T-backtest+1):T)],
+ start = c(1989, 9), frequency = 12) > fit1s=arima(subA7d,order=c(p=12,d=0,q=0), + include.mean =FALSE)

pour le premier modèle, et

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

pour le second. Ensuite, on peut faire de la prévision, et comparer les prévisions avec les données écartées pour l'estimation du modèle.,

> XA7p1=forecast(fit1s,h=backtest)
> XA7p2=forecast(fit2s,h=backtest)
> A7dobs=A7d[((T-backtest+1):T)]
> (M=cbind(observé=A7dobs,modèle1=as.numeric(XA7p1$mean),
+ modèle2=as.numeric(XA7p2$mean))) observé modèle1 modèle2 [1,] -13610.071 -12264.7138 -13558.071 [2,] -10108.071 -12236.6227 -10410.071 [3,] -8034.071 -9333.8365 -9504.071 [4,] 3899.929 4071.9313 4108.929 [5,] 2009.929 901.3745 2661.929 [6,] 5253.929 7610.9954 5866.929 [7,] 25357.929 29232.2317 30317.929 [8,] 35229.929 29641.5172 31737.929 [9,] 4341.929 5603.8008 4940.929

On peut alors retenir un critère classique, comme l'erreur quadratique, et donc comparer la somme des carrées des erreurs,

> sum((M[,1]-M[,2])^2)
[1] 62677238
> sum((M[,1]-M[,3])^2)
[1] 40253827

Le second modèle semble ici avoir été meilleur, sur les 9 derniers mois. On peut aussi visualiser les erreurs de prévision,

> (TPS=tsp(subA7d))
[1] 1989.667 1995.917   12.000
> temps=seq(TPS[2]+1/TPS[3],by=1/TPS[3],length=backtest)
> rangex=c(TPS[1],TPS[2]+backtest/TPS[3])
> plot(subA7d,type="l",xlim=rangex)
> abline(v=T-11.5,lty=2)
> polygon(c(temps,rev(temps)),c(A7dobs,rev(XA7p1$mean)),
+ col="yellow",border=NA) > lines(temps,A7dobs) > lines(temps,XA7p1$mean,col="red")

pour le premier modèle,

et

> polygon(c(temps,rev(temps)),c(A7dobs,rev(XA7p2$mean)),
+ col="yellow",border=NA) > lines(temps,XA7p2$mean,col="red")

pour le second,

Tuesday, September 18 2012

Examen final ACT6420

Chose promise, chose due: les deux bases de données qui seront utilisées pour l'examen sont désormais en ligne. 10 questions porteront sur l'analyse de la série temporelle suivante

serie=read.table("http://freakonometrics.free.fr/exam-ts.txt")
X=serie$x
XTS=serie=ts(X,start=c(1999,9),frequency=12)

(qui peut être modélisée soit comme un objet série temporel, soit comme un vecteur numérique).

Enfin, 10 questions porteront sur la base de données suivante

examen=read.table(
"http://freakonometrics.blog.free.fr/public/
data/basket-exam-v2.csv"
, header=TRUE,sep=";")

(si l'importation ne marche pas, veillez à ce que l'adresse du fichier tienne sur une ligne). Cette base devait servir initialement pour l'examen intra. Elle contient des matchs de baskets universitaires (aux États-Unis). Les variables sont les suivantes

  • EQUIPE le nom de l'équipe qui sert de référence
  • HALFDIFF le nombre de points de différence à la mi-temps (négatif signifie que l'EQUIPE perdait)
  • WINNER vaut 1 si l'équipe a gagné le match, 0 sinon,
  • TOTALPOINTSDIFF le nombre de points de différence à la fin du match
  • SECONDHALFPOINTS est le nombre de points marqués par l'EQUIPE pendant la seconde mi-temps
  • NODIFF_PTS_1000-LEFT la différence de points entre les équipes 10 minutes avant la fin du match (au milieu de la seconde période)
  • HOMETEAM vaut 0 si l'équipe joue à l'extérieur, 1 si elle joue à domicile
  • HALFWINNER vaut 1 si l'EQUIPE gagnait (supérieur ou nul) à la mi-temps
  • HALFLOOSER vaut 1 si l'EQUIPE perdait (inférieur ou nul) à la mi-temps
  • SECONDHALFPOINTSOTHER est le nombre de points marqués par l'adversaire au cours de la seconde mi-temps

La variable d'intérêt sera le nombre de points d'écart à la fin du match.

Enfin, 20 questions supplémentaires porteront sur la compréhension générale du cours. Les calculatrices seront autorisées. Maintenant, pour tous ceux qui ont encore du mal avec R, je recommande

  • “R pour les débutants”d'Emmanuel Paradis (PDF)
  • “Brise Glace-R” d'Andrew Robinson er Arnaud Schloesing (PDF)
  • “Introduction à la programmation en R” de Vincent Goulet (PDF)

(les questions d'analyse de sortie ne porteront pas sur la programmation en R, mais davantage sur l'analyse critique des sorties, et leur interprétation).

Monday, September 17 2012

Ordres d'un processus ARMA

Dans la méthodologie de Box & Jenkins, une étape qui arrive très rapidement est le choix des ordres du processus http://freakonometrics.blog.free.fr/public/perso6/ARMA__latex.gif, une fois que l'on validé l'hypothèse de stationnarité de la série, comme on l'a vu en cours la semaine passée. Considérons la série du trafic autoroutier,

> autoroute=read.table(
+"http://freakonometrics.blog.free.fr/public/data/autoroute.csv",
+ header=TRUE,sep=";")
> a7=autoroute$a007
> A7=ts(a7,start = c(1989, 9), frequency = 12)
Un outils pratique de sélections des ordres http://freakonometrics.blog.free.fr/public/perso6/platex.gif et http://freakonometrics.blog.free.fr/public/perso6/qlatex.gif dans un modèle http://freakonometrics.blog.free.fr/public/perso6/ARMAlatex.gif est la fonction d'autocorrélation étendue. La définition est donnée dans les notes de cours (Def. 223) à partir des statistiques proposées par Tsay & Tiao (1984),
>  EACF=eacf(A7)
AR/MA
  0 1 2 3 4 5 6 7 8 9 10 11 12 13
0 x o o x x x x x o o x  x  x  o
1 x x o o x x x o o o x  x  x  x
2 o x o o o o x o o o o  x  x  o
3 o x o x o o o o o o o  x  x  x
4 x x x x o o o o o o o  x  o  o
5 x x o o o x o o o o o  x  o  x
6 x x o o o x o o o o o  x  o  o
7 o x o o o x o o o o o  x  o  o
>  EACF
$eacf
           [,1]       [,2]        [,3]        [,4]
[1,]  0.6476234  0.2124105 -0.02413173 -0.24234535
[2,]  0.4889076  0.2797152 -0.02494135 -0.06094037
[3,]  0.1006541 -0.2285771  0.03514148  0.08000588
[4,]  0.1390240 -0.2788742  0.04386746  0.28307260
[5,] -0.5091680  0.3144899 -0.34572269  0.31450865
[6,] -0.4224571 -0.4877505  0.16054232 -0.09130728
[7,] -0.4731353 -0.4324857 -0.04847184 -0.10500350
[8,] -0.2129591 -0.4072901  0.09487899 -0.06493243
             [,5]        [,6]        [,7]        [,8]
[1,] -0.514330187 -0.61634046 -0.52314403 -0.28008661
[2,] -0.254912957 -0.28966664 -0.33963243 -0.21863077
[3,] -0.156624357 -0.01199786 -0.25116738  0.13079231
[4,] -0.183283544  0.03651508 -0.08711829  0.11626377
[5,] -0.190885091 -0.09786463 -0.09182557  0.08818875
[6,] -0.072904071  0.29271777 -0.09334712  0.01972648
[7,] -0.009873289  0.36909726  0.01698660 -0.03317456
[8,]  0.020485930  0.38342158  0.16981715 -0.02592442
             [,9]       [,10]       [,11]     [,12]
[1,] -0.082607058  0.15178887  0.56583179 0.8368975
[2,] -0.085026400  0.02731460  0.26158357 0.7844748
[3,]  0.001866994 -0.11658312  0.03038621 0.7026361
[4,]  0.025183793 -0.21608692  0.05660781 0.6674301
[5,] -0.006831894 -0.02514440 -0.07390257 0.4762563
[6,]  0.010058718 -0.03888613 -0.04382043 0.5338091
[7,]  0.032124905 -0.07022090 -0.04427400 0.4674165
[8,]  0.024189179 -0.20818201  0.01459933 0.4830369
          [,13]       [,14]
[1,]  0.5637439  0.17862571
[2,]  0.4530716  0.24413569
[3,]  0.2534178 -0.20160890
[4,]  0.2409861 -0.29462510
[5,] -0.1517324 -0.14763294
[6,] -0.1701182 -0.28771495
[7,] -0.1651515  0.05466457
[8,] -0.1403198 -0.04095030
 
$ar.max
[1] 8
 
$ma.ma
[1] 14
Les ronds dans la matrice désignent des valeurs non-significatives. Par défaut, le nombre de retards, pris en compte pour la composante autorégressive est faible, mais on peut l'augmenter.
> EACF=eacf(A7,13,13)
AR/MA
   0 1 2 3 4 5 6 7 8 9 10 11 12 13
0  x o o x x x x x o o x  x  x  o
1  x x o o x x x o o o x  x  x  x
2  o x o o o o x o o o o  x  x  o
3  o x o x o o o o o o o  x  x  x
4  x x x x o o o o o o o  x  o  o
5  x x o o o x o o o o o  x  o  x
6  x x o o o x o o o o o  x  o  o
7  o x o o o x o o o o o  x  o  o
8  o x o o o x o o o o o  x  o  o
9  x x o o x o x o o o o  x  o  o
10 x x o o x o x o o o o  x  o  o
11 x x o o x x x o o o x  x  o  o
12 x x o o o o o o o o o  o  o  o
13 x x o o o o o o o o o  o  o  o
> EACF$eacf
      [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9] [,10]
[1,]  0.64  0.21 -0.02 -0.24 -0.51 -0.61 -0.52 -0.28 -0.08  0.15
[2,]  0.48  0.27 -0.02 -0.06 -0.25 -0.28 -0.33 -0.21 -0.08  0.02
[3,]  0.10 -0.22  0.03  0.08 -0.15 -0.01 -0.25  0.13  0.00 -0.11
[4,]  0.13 -0.27  0.04  0.28 -0.18  0.03 -0.08  0.11  0.02 -0.21
[5,] -0.50  0.31 -0.34  0.31 -0.19 -0.09 -0.09  0.08  0.00 -0.02
[6,] -0.42 -0.48  0.16 -0.09 -0.07  0.29 -0.09  0.01  0.01 -0.03
[7,] -0.47 -0.43 -0.04 -0.10  0.00  0.36  0.01 -0.03  0.03 -0.07
[8,] -0.21 -0.40  0.09 -0.06  0.02  0.38  0.16 -0.02  0.02 -0.20
[9,] -0.14 -0.50 -0.10 -0.06  0.11  0.38 -0.01 -0.02 -0.08 -0.20
[10,] -0.29  0.48  0.00  0.18  0.27 -0.12  0.41  0.00  0.16  0.15
[11,]  0.24  0.48 -0.04  0.02  0.46  0.10  0.37 -0.10  0.05  0.16
[12,] -0.59  0.49 -0.18  0.05  0.31 -0.32  0.34 -0.16  0.07  0.16
[13,] -0.31  0.31 -0.12  0.16 -0.11  0.04 -0.04  0.12 -0.09  0.00
[14,]  0.47  0.26  0.11  0.13 -0.01 -0.01  0.00  0.07 -0.03  0.02
On peut aussi visualiser graphiquement les différentes valeurs, avec en ordonnée les ordres autoréregressifs (http://freakonometrics.blog.free.fr/public/perso6/platex.gif) et en abscisse les ordres moyenne mobile (http://freakonometrics.blog.free.fr/public/perso6/qlatex.gif). 
> library(RColorBrewer)
> CL=brewer.pal(6, "RdBu")
> ceacf=matrix(as.numeric(cut(EACF$eacf,
+ ((-3):3)/3,labels=1:6)),nrow(EACF$eacf),
+ ncol(EACF$eacf))
> for(i in 1:ncol(EACF$eacf)){
+ for(j in 1:nrow(EACF$eacf)){
+ polygon(c(i-1,i-1,i,i)-.5,c(j-1,j,j,j-1)-.5,
+ col=CL[ceacf[j,i]])
+ }}

Un bruit blanc est en bas à gauche (http://freakonometrics.blog.free.fr/public/perso6/platex.gif et http://freakonometrics.blog.free.fr/public/perso6/qlatex.gif nuls). Dans cette méthode, dite méthode des coins, on cherche un coin telle que dans le quadrant supérieur, les valeurs soient non-significatives (ternes sur le dessin ci-dessus). La figure ci-desssous correspond à l'analyse d'un modèle http://freakonometrics.blog.free.fr/public/perso6/ARMA126latex.gif

Les fortes valeurs positives sont en bleu foncé, les fortes valeurs négatives sont en rouge foncé. On peut aussi regarder la fonction suivante, qui utilise une identification de modèle par MINIC (Minimum Information Criterion) à l'aide du critère de Shcwarz (BIC, ou SBC Schwarz's Bayesian Criterion)
> armaselect(A7,nbmod=5)
p q      sbc
[1,] 12 1 1441.798
[2,] 13 0 1442.628
[3,] 12 0 1443.188
[4,] 12 2 1443.362
[5,] 14 0 1445.069
Enfin, une dernière fonction possible est évoquée dans la section 6.5. du livre de Cryer & Chan (2008),
> ARMA.SELECTION=
+ armasubsets(A7,nar=14,nma=14,ar.method='ols')
> plot(ARMA.SELECTION)

basé sur le critère de Schwarz.
Je vais me répéter, mais ces méthodes ne sont que des outils, histoire d'avoir des pistes si on ne sait trop dans quelle direction partir. Et compte tenu de la saisonnalité de la série, je serais pour ma part parti sur un modèle avec http://freakonometrics.blog.free.fr/public/perso6/platex.gif=12, histoire de voir si on ne pourrait pas avoir un modèle simple, et facilement interprétable en plus...

Wednesday, March 21 2012

La maudite constante des modèles ARIMA

Dans les modèles ARIMA, autant le dire tout de suite, les constantes c'est pénible ! Pour comprendre un peu mieux ce qui se passe, considérons ici un modèle ARMA avec constante. Ou pour commencer, juste un AR(1)
Supposons que la série  soit stationnaire, de moyenne . Alors en prenant l’espérance de part et d'autre dans l'expression précédente, , i.e. , ou . Plus généralement, avec un ARMA, , en prenant la encore l’espérance des deux cotés, on en déduit .

Simulons un processus AR(1),

Pour simuler des processus ARIMA, on pourrait utiliser la commande

> X=arima.sim(list(order=c(1,0,0),ar=1/3),n=1000)+2
> mean(X)
[1] 1.931767
mais comme le but est de comprendre ce qui se passe, autant faire les choses calmement,
> X=rep(NA,1010)
> X[1]=0
> for(t in 2:1010){X[t]=4/3+X[t-1]/3+rnorm(1)}
> X=X[-(1:10)]
> mean(X)
[1] 2.03397
I.e. ici le processus est de moyenne qui vaut 2.


Regardons maintenant ce que donnerait l'estimation du processus AR(1),

> arima(X, order = c(1, 0, 0))
 
Call:
arima(x = X, order = c(1, 0, 0))
 
Coefficients:
ar1  intercept
0.3738     2.0334
s.e.  0.0294     0.0487
 
sigma^2 estimated as 0.9318:  log likelihood = -1383.68
De manière un peu surprenante, le coefficient appelé intercept n'est pas la constante dans le modèle AR(1), mais la moyenne du processus. Autrement dit, R n'ajuste pas un processus
comme le laisserait penser l'intuition, mais un processus
Ces deux formes sont bien entendu équivalentes. Mais les coefficients estimés ne sont pas tout à fait ce que l'on attendait...
Plaçons nous maintenant dans le cas d'un processus non-stationnaire. L'extension naturelle serait de considérer un processus ARIMA(1,1,0), ou bien un processus tel que la différence soit un processus AR(1). Un processus ARIMA(1,1,0) avec une constante s'écrirait
en utilisant l'opérateur retard. Ceci fait penser à un modèle avec une tendance linéaire. Posons  (afin de se débarrasser de cette tendance). Alors
i.e.
ou encore
On note que  ou encore  suit alors un processus ARIMA(1,1,0) sans constante cette fois.
Afin de visualiser ce que donnerait l'inférence, simulons le processus suivant, qui est un processus AR(1) avec constante que l'on intègre:  avec
Peut-on retrouver les différents paramètres du processus avec R ?
Commençons (là encore) par simuler un tel processus,
> U=rep(NA,1010)
> U[1]=0
> for(t in 2:1010){U[t]=4/3+U[t-1]/3+rnorm(1)}
> U=U[-(1:10)]
> X=cumsum(U)
Sous R, on obtient l'estimation suivant si l'on tente de calibrer un modèle ARIMA(1,1,0)
> arima(X, order = c(1, 1, 0))
 
Call:
arima(x = X, order = c(1, 1, 0))
 
Coefficients:
ar1
0.8616
s.e.  0.0160
  
sigma^2 estimated as1.343:log likelihood = -1565.63

Ça ne convient pas du tout.... On peut tenter un processus AR(1) (avec constante) sur la série différenciée...

> arima(diff(X), order = c(1, 0, 0))
 
Call:
arima(x = diff(X), order = c(1, 0, 0))
 
Coefficients:
ar1  intercept
0.3564     2.0200
s.e.  0.0295     0.0486
 
sigma^2 estimated as 0.9782:  log likelihood = -1406.6
On progresse, sauf que comme auparavant, le terme qui est donne n'est pas la constante dans le modèle ARIMA, mais la moyenne du processus différencié. Mais cette fois, on a un interprétation, c'est que la constante est la pente de la tendance ! Si on estime la pente associée a , on récupère la même valeur:
> arima(X, order = c(1, 1, 0), xreg=1:length(X))
 
Call:
arima(x = X, order = c(1, 1, 0), xreg = 1:length(X))
 
Coefficients:
ar1  1:length(X)
0.3566       2.0519
s.e.  0.0296       0.0487
 
sigma^2 estimated as 0.9787:  log likelihood = -1406.82
Sur la figure ci-dessous, on retrouve le fait qu'en enlevant la tendance linéaire à  donneune série intégrée, sans constante (qui fait penser à une marche aléatoire)

Autrement dit

  • dans le cas d'une série stationnaire, la constante estimée n'est pas du tout la constante, mais la moyenne de la série temporelle
  • dans le cas d'une série non-stationnaire, la constante estimée dans la série différenciée a du sens, au sens ou il s'agit de la pente de la tendance (linéaire) du processus
Avant de conclure, une petite remarque. Quid de la prévision ? Si on commencer par une révision à l'aide du premier processus ARIMA, on obtient une prédiction pour avec un énorme intervalle de confiance (sans aucun bon sens)
> ARIMA1=arima(X, order = c(1, 1, 0))
> ARIMA2=arima(X, order = c(1, 1, 0), xreg=1:length(X))
> Xp1=predict(ARIMA1,20)
> Xp2=predict(ARIMA2,20,newxreg=
+ (length(X)+1):(length(X)+20))
> plot(960:1000,X[960:1000],xlim=c(960,1020),type="l")
> polygon(c(1001:1020,rev(1001:1020)),
+ c(Xp1$pred+2*Xp1$se,rev(Xp1$pred-2*Xp1$se)),
+ col=CL[3],border=NA)
> lines(1001:1020,Xp1$pred,col="red",lwd=2)
intervalle qui se visualise sur le graphique ci-dessous,

Si on regarde pour l'autre modèle,
> lines(1001:1020,Xp2$pred,col="blue",lwd=2)

Tous ceux qui auront reconnu ici des problèmes qui se posent dans la modélisation de la série http://freakonometrics.blog.free.fr/public/maths/viekt.png dans le modèle de Lee & Carter (1992) auront compris qu'il existe une vraie application a tout ce que je viens de raconter (je pense au particulier au commentaire poste par @ClaudeT en début de semaine) car la série des http://freakonometrics.blog.free.fr/public/maths/viekt.png ressemble très fortement à ce genre de série, non-stationnaire avec une tendance linéaire. Et qu'un modèle mal spécifié laisse à penser que l'incertitude est beaucoup beaucoup plus grande que ce qu'elle est vraiment (en plus de donner une prédiction surprenante à long terme !). Donc avant de faire tourner rapidement les codes, il vaut mieux regarder attentivement ce qu'ils font réellement...

Friday, February 24 2012

Sur le lissage exponentiel

Sur le lissage exponentiel, je pourrais renvoyer vers Gardner (1985), et la version remise au goût du jour, Gardner (2005). J'avais pris comme position, dans le cours, de présenter rapidement le lissage exponentiel, en notant qu'on les évoquerais à nouveau plus tard comme cas particulier des modèles ARIMA. Comme le note la dernière version, Gardner (2005), ce n'est pas forcément l'unique manière de voir, et le cours pourrait être orienté complètement autour des notions de lissage exponentiel (c'est d'ailleurs le point du vue de livre Hyndman, Koehler, Ord & Snyder (2008)), "when Gardner (2005) appeared, many believed that exponential smoothing should be disregarded because it was either a special case of ARIMA modeling or an ad hoc procedure with no statistical rationale. As McKenzie (1985) observed, this opinion was expressed in numerous references to my paper. Since 1985, the special case argument has been turned on its head, and today we know that exponential smoothing methods are optimal for a very general class of state-space models that is in fact broader than the ARIMA class."

N'étant toujours pas convaincu, j'évoquerais à nouveau le lissage exponentiel quand on présentera les modèles ARIMA. En attendant, un peu de code pour mieux comprendre ce qui est fait quand on fait du lissage exponentiel.

Commençons par le lissage exponentiel simple, i.e.

http://freakonometrics.blog.free.fr/public/perso5/lex03.gif

http://freakonometrics.blog.free.fr/public/perso5/le04.gif désigne un poids attribué à la nouvelle observation dans la fonction de lissage http://freakonometrics.blog.free.fr/public/perso5/lex05.gif. Le code pour lisser une série est alors assez simple,

> library(datasets)
> X=as.numeric(Nile)
> Lissage=function(a){
+  T=length(X)
+  L=rep(NA,T)
+  L[1]=X[1]
+  for(t in 2:T){L[t]=a*X[t]+(1-a)*L[t-1]}
+  return(L)
+ }
> plot(X,type="b",cex=.6)
> lines(Lissage(.2),col="red")

http://freakonometrics.blog.free.fr/public/perso5/lissage-exp-1.gif

Sur la figure ci-dessus, non visualise l'impact de http://freakonometrics.blog.free.fr/public/perso5/lisse16.gif sur le lissage. La prévision que l'on fait à la date http://freakonometrics.blog.free.fr/public/perso5/lex14.gif, pour un horizon http://freakonometrics.blog.free.fr/public/perso5/lisse12.gif est alors http://freakonometrics.blog.free.fr/public/perso5/lisse11.gif. Il est alors possible de voir le poids comme un paramètre et on va alors chercher le poids optimal. La stratégie classique est de minimiser l'erreur de prédiction commise à un horizon de 1

http://freakonometrics.blog.free.fr/public/perso5/lisse14.gif

> V=function(a){
+  T=length(X)
+  L=erreur=rep(NA,T)
+  erreur[1]=0
+  L[1]=X[1]
+  for(t in 2:T){
+  L[t]=a*X[t]+(1-a)*L[t-1]
+  erreur[t]=X[t]-L[t-1] }
+ return(sum(erreur^2))
+ }

Ici, on obtient comme poids optimal

> A=seq(0,1,by=.02)
> Ax=Vectorize(V)(A)
> plot(A,Ax,ylim=c(min(Ax),min(Ax)*1.05))
> optimize(V,c(0,.5))$minimum
[1] 0.246581

On notera que c'est ce que suggère la fonction de R,

> hw=HoltWinters(X,beta=FALSE,gamma=FALSE,l.start=X[1])
> hw
Holt-Winters exponential smoothing without trend an seasonal comp.
 
Call:
HoltWinters(x = X, beta = FALSE, gamma = FALSE, l.start = X[1])
 
Smoothing parameters:
alpha:  0.2465579
beta :  FALSE
gamma:  FALSE
 
Coefficients:
[,1]
a 805.0389
 
> plot(hw)
> points(2:(length(X)+1),Vectorize(Lissage)(.2465),col="blue")

Dans le cas du lissage exponentiel double, i.e.

http://freakonometrics.blog.free.fr/public/perso5/lex50.gif

et

http://freakonometrics.blog.free.fr/public/perso5/HW02.gif

Dans ce cas, la prédiction est http://freakonometrics.blog.free.fr/public/perso5/lisse13.gif. Le code pour faire un lissage double est là encore assez simple,
> Lissage=function(a,b){
+  T=length(X)
+  L=B=rep(NA,T)
+  L[1]=X[1]; B[1]=0
+  for(t in 2:T){
+  L[t]=a*X[t]+(1-a)*(L[t-1]+B[t-1])
+  B[t]=b*(L[t]-L[t-1])+(1-b)*B[t-1] }
+ return(L)
+ }
Sur la figure suivante, on visualise l'évolution du lissage en fonction de http://freakonometrics.blog.free.fr/public/perso5/lisse15.gif (http://freakonometrics.blog.free.fr/public/perso5/lisse16.gif étant ici fixé),
http://freakonometrics.blog.free.fr/public/perso5/lissage-exp-2.gif
(le lissage simple - avec le même poids http://freakonometrics.blog.free.fr/public/perso5/lisse16.gif - apparaissant ici en trait clair).

Méthodes de prévisions, suite des références

Pour la seconde partie du cours ACT6420 - méthodes de prévisions - on va aborder la modélisation des séries temporelles (la première partie portait sur l'analyse des données individuelles avec les références mentionnés dans un précédant billet).

On passera surtout du temps sur les modèles ARIMA, et le livre de référence sur le sujet sera le Time Series Analysis de James Hamilton (ou au moins les premiers chapitres). Ces chapitres serviront de références théoriques pour tout le cours. Pour les aspects plus pratiques, je renvoie vers Time Series Analysis and Its Applications With R Examples de Robert Shumway et David Stoffer, ou Introductory Time Series with R de Paul Cowpertwait et Andrew Metcalfe. Je voudrais aussi citer la page http://a-little-book-of-r-for-time-series.readthedocs.org qui propose beaucoup de codes R utile pour modéliser les séries temporelles. Sinon, le livre d'Edward Frees propose aussi deux (courts) chapitres sur les séries temporelles. Je vais engin retravailler dans les semaines à venir les notes de cours que j'avais écrites il y a quelques années .

Wednesday, December 15 2010

I really need to find hot (and sexy) topics

50 days ago (here), I was supposed to be very optimistic about the probability that I could reach a million viewed pages on that blog (over a bit more than two years). Unfortunately, the wind has changed and today, the probability is quite low...

 base=read.table("millionb.csv",sep=";",header=TRUE)
X1=cumsum(base$nombre)
base=read.table("million2b.csv",sep=";",header=TRUE)
X2=cumsum(base$nombre)
X=X1+X2
D=as.Date(as.character(base$date),"%m/%d/%Y") kt=which(D==as.Date("01/06/2010","%d/%m/%Y")) D0=as.Date("08/11/2008","%d/%m/%Y") D=D0+1:length(X1) P=rep(NA,(length(X)-kt)+1) for(h in 0:(length(X)-kt)){ model <- arima(X[1:(kt+h)],c(7 1,7),method="CSS")
forecast <- predict(model,200) u=max(D[1:kt+h])+1:300 k=which(u==as.Date("01/01/2011","%d/%m/%Y")) (P[h+1]=1-pnorm(1000000,forecast$pred[k],forecast$se[k])) } plot( D[length(D)-length(P)]+1:220,c(P,rep(NA,220-length(P))), ylab="Probability to reach 1,000,000",xlab="", type="l",col="red",ylim=c(0,1))
So, I guess my posts on multiple internal rates of return, or Young's inequality will have to wait next year... I really need to find some more sexy post to attract readers.. Challenge accepted !

Tuesday, December 7 2010

Statistique de l'assurance STT6705V, partie 12 bis

In the previous post (here) discussing forecasts of actuarial quantities, I did not mention much how to forecast the temporal component in the Lee-Carter model. Actually, many things can be done. Consider here some exponential smoothing techniques (see here for more details). For instance, with a simple model, with an additive error, no trend and no seasonal component,

or equivalently

Then http://freakonometrics.blog.free.fr/public/maths/ETS-NN-3.png, and if we want a confidence interval
> (ETS.ANN=ets(Y,model="ANN"))
ETS(A,N,N)
 
Call:
ets(y = Y, model = "ANN")
 
Smoothing parameters:
alpha = 0.7039
 
Initial states:
l = 75.0358
 
sigma: 11.4136
 
AIC AICc BIC
941.0089 941.1339 946.1991
> plot(forecast(ETS.ANN,h=100))
Here, a natural idea will probably be to take into account a linear trend on the series, i.e.

where

then the forecast we make is  i.e.
ETS.AAN=ets(Y,model="AAN")
plot(forecast(ETS.AAN,h=100))
It is also possible to ask for an "optimal" exponential smoothing model (without any specification)
> ETS
ETS(A,A,N)
 
Call:
ets(y = Y)
 
Smoothing parameters:
alpha = 0.6107
beta = 1e-04
 
Initial states:
l = 74.5622
b = -1.9812
 
sigma: 11.0094
 
AIC AICc BIC
937.8695 938.2950 948.2500
> plot(forecast(ETS,h=100))

Another idea can be to consider some autoregressive integrated moving average models (ARIMA), here with a linear trend

For instance, if we want only the integrated component, i.e.

then the code is
> ARIMA010T=Arima(Y,order=c(0,1,0),include.drift=TRUE)
Series: Y
ARIMA(0,1,0) with drift
 
Call: Arima(x = Y, order = c(0, 1, 0), include.drift = TRUE)
 
Coefficients:
drift
-2.0337
s.e. 1.1904
 
sigma^2 estimated as 138.9: log likelihood = -380.8
AIC = 765.6 AICc = 765.73 BIC = 770.77
> plot(forecast(ARIMA010T,h=100))
But note that any kind of ARIMA model can be considered, e.g. integrated with an autoregressive component (here of order 1)
or with also a moving average component (again of order 1)
Actually, note that, once again, we can ask for an automatic selection of the model,
> (autoARIMA=auto.arima(Y,allowdrift=TRUE))
Series: Y
ARIMA(0,1,1) with drift
 
Call: auto.arima(x = Y, allowdrift = TRUE)
 
Coefficients:
ma1 drift
-0.3868 -2.0139
s.e. 0.0970 0.6894
 
sigma^2 estimated as 122.3: log likelihood = -374.65
AIC = 755.29 AICc = 755.55 BIC = 763.05
> plot(forecast(autoARIMA,h=100))

Finally, it is possible to use also also some specific functions of the package, for instance to consider a random walk with a drift,
RWF=rwf(Y,h=100,drift=TRUE)
plot(RWF)
or Holt model (with a trend)
HOLT=holt(Y,h=100,damped=TRUE)
plot(HOLT)
And if all that is not enough, it is also possible to go further by changing the size of the series we use to fit the model. A question that naturally arises is the treatment of wars in our model: shouldn't we assume that the forecast should be based only on the last 50 years (and exclude wars from our analysis) ? In that case, for instance, the exponential smoothing technique gives
while the ARIMA procedure returns
And finally, with the Holt technique, we have
So, it looks like we have a lot of techniques that can be used to provide a forecast for the temporal component in the Lee-Carter model,
All those estimators can be used to estimate annuities of insurance contracts (as here),
BASEB=BASEB[,1:99]
BASEC=BASEC[,1:99]
AGE=AGE[1:99]
library(gnm)
D=as.vector(BASEB)
E=as.vector(BASEC)
A=rep(AGE,each=length(ANNEE))
Y=rep(ANNEE,length(AGE))
base=data.frame(D,E,A,Y,a=as.factor(A),
y=as.factor(Y))
LC2=gnm(D~a+Mult(a,y),offset=log(E),
family=poisson,data=base)
 
A=LC2$coefficients[1]+LC2$coefficients[2:99]
B=LC2$coefficients[100:198]
K0=LC2$coefficients[199:297]
Y=as.numeric(K01)
K1=c(K0,forecast(ets(Y),h=120)$mean)
K2=c(K0,forecast(auto.arima(Y,allowdrift=TRUE),h=120)$mean)
K3=c(K0,rwf(Y,h=120,drift=TRUE)$mean)
K4=c(K0,holt(Y,h=120,drift=TRUE)$mean)
K5=c(K0,forecast(ets(Y[50:99]),h=120)$mean)
K6=c(K0,forecast(auto.arima(Y[50:99],allowdrift=TRUE),h=120)$mean)
K7=c(K0,rwf(Y[50:99],h=120,drift=TRUE)$mean)
K8=c(K0,holt(Y[50:99],h=120,drift=TRUE)$mean)
 
MU=matrix(NA,length(A),length(K1))
MU1=MU2=MU3=MU4=MU5=MU6=MU7=MU8=MU
 
for(i in 1:length(A)){
for(j in 1:length(K1)){
MU1[i,j]=exp(A[i]+B[i]*K1[j])
MU2[i,j]=exp(A[i]+B[i]*K5[j])
MU3[i,j]=exp(A[i]+B[i]*K3[j])
MU4[i,j]=exp(A[i]+B[i]*K4[j])
MU5[i,j]=exp(A[i]+B[i]*K5[j])
MU6[i,j]=exp(A[i]+B[i]*K6[j])
MU7[i,j]=exp(A[i]+B[i]*K7[j])
MU8[i,j]=exp(A[i]+B[i]*K8[j])
}}
 
t=2000
x=40
s=seq(0,98-x-1)
 
Pxt1=cumprod(exp(-diag(MU1[x+1+s,t+s-1898])))
Pxt2=cumprod(exp(-diag(MU2[x+1+s,t+s-1898])))
Pxt3=cumprod(exp(-diag(MU3[x+1+s,t+s-1898])))
Pxt4=cumprod(exp(-diag(MU4[x+1+s,t+s-1898])))
Pxt5=cumprod(exp(-diag(MU5[x+1+s,t+s-1898])))
Pxt6=cumprod(exp(-diag(MU6[x+1+s,t+s-1898])))
Pxt7=cumprod(exp(-diag(MU7[x+1+s,t+s-1898])))
Pxt8=cumprod(exp(-diag(MU8[x+1+s,t+s-1898])))
 
r=.035
m=70
h=seq(0,21)
V1=1/(1+r)^(m-x+h)*Pxt1[m-x+h]
V2=1/(1+r)^(m-x+h)*Pxt2[m-x+h]
V3=1/(1+r)^(m-x+h)*Pxt3[m-x+h]
V4=1/(1+r)^(m-x+h)*Pxt4[m-x+h]
V5=1/(1+r)^(m-x+h)*Pxt5[m-x+h]
V6=1/(1+r)^(m-x+h)*Pxt6[m-x+h]
V7=1/(1+r)^(m-x+h)*Pxt7[m-x+h]
V8=1/(1+r)^(m-x+h)*Pxt8[m-x+h]
Hence, here the difference is significant,
> M=cbind(V1,V2,V3,V4,V5,V6,V7,V8)
> apply(M,2,sum)
V1       V2       V3       V4       V5       V6       V7       V8
4.389372 4.632793 4.406465 4.389372 4.632793 4.468934 4.482064 4.632793
or graphically,

Wednesday, October 27 2010

A million ? what are the odds...

50 days ago, I published a post, here, on forecasting techniques. I was wondering what could be the probability to have, by the end of this year, one million pages viewed (from Google Analytics) on this blog. Well, initially, it was on my blog at the Université de Rennes 1 (http://blogperso.univ-rennes1.fr/arthur.charpentier/), but since I transfered the blog, I had to revise my code. Initially, I had that kind of graphs,

and when I look at the cumulative distribution of the number of pages viewed on January first, I had

while for the distribution of the time I should read this million (the dual problem), I obtained

and I said that I should have around 35% chance to reach the million pages viewed by the end of this year.
Here is the updated graph, with the blog à Université de Rennes 1 (still in black) and the one here (in blue, where I add the two blogs together).

Actually, I decided to look at the evolution of the probability to reach the million by New Year's Eve...

The code looks like that,
> base=read.table("http://perso.univ-rennes1.fr/arthur.charpentier/million2.csv",
+ sep=";",header=TRUE)
> X2=cumsum(base$nombre)
> X=X1+X2
> kt=which(D==as.Date("01/06/2010","%d/%m/%Y"))
> D0=as.Date("08/11/2008","%d/%m/%Y")
> D=D0+1:length(X1)
> P=rep(NA,(length(X)-kt)+1)
> for(h in 0:(length(X)-kt)){
+ model  <- arima(X[1:(kt+h)],c(7 ,   # partie AR
+                     1,    # partie I
+                     7),method="CSS")   # partie MA
+ forecast <- predict(model,200)
+ u=max(D[1:kt+h])+1:300
+ k=which(u==as.Date("01/01/2011","%d/%m/%Y"))
+ (P[h+1]=1-pnorm(1000000,forecast$pred[k],forecast$se[k]))
+ }
It has been a bit tricky, since I wanted an automatic fit of the ARIMA process, meaning that I had to assess a priori the orders of the ARIMA process. And I had numerical problems, since we got non stationary AR part at least at one period of time considered.... So finally I used here the CSS method which uses conditional-sum-of-squares to find starting values in the optimization procedure.

Actually, if we consider a classical descritption of traders, it looks like I act as a trader (dealing with millions and forgetting about real people): it is the same here, I do not know what a million means, I cannot imagine 250,000 visitors looking at that blog... But I can still do the maths. Anyway, a million is huge when I start to think about it... but perhaps I should not... I cannot possibility imagine that so many people might find interesting my mathematical lucubration*....
* initially I was looking for the analogous of "élucubration" in French, meaning "divagation, absurd theory" (the proper translation might be "rantings" (here) , "ravings" (here) or "wild imagining" (everywhere else here or there)). When I asked Google for a possible translation (here), I got "lucubration" which means "composed by night; that which is produced by meditation in retirement". Well, it was not initially what I intended to say, but since I usually work on my blog during the night, when I got awake by one of the girls, I decided to keep this word.... At least, I learnt something today, appart for the code mentioned above....

Thursday, September 9 2010

Le million ! le milllion !

Hier soir (ou ce matin, je suis perdu avec ce décalage horaire) Christelle me demandait de parler un peu de prévision avec R (ici). Au lieu de renvoyer vers l'aide en ligne, penons un exemple pratique (et simple, voire si possible intéressant): la fréquentation d'un blog (en l'occurrence http://blogperso.univ-rennes1.fr/arthur.charpentier/). Considérons le nombre de pages vues, par jour, selon Google Analytics.
> base=read.table("http://perso.univ-rennes1.fr/arthur.charpentier/million.csv",
+ sep=";",header=TRUE
)> X=base$nombre
> D0=as.Date("08/11/2008","%d/%m/%Y")
> D=D0+1:length(X)
> plot(D,X)
> plot(acf(X,lag=90))
La série a l'allure suivante (oui, le compteur n'a été installé qu'il y a deux ans),

Les autocorrélations montre une forte saisonnalité hebdomadaire, avec moins de consultations en semaine que le week-end,

La question que l'on cherche à résoudre est "ai-je une chance d'atteindre le million de pages vues d'ici la fin de l'année ?".

On peut traduire cette question de deux manières
  • quelle est la probabilité que le 1er janvier, j'ai atteint le million de pages vues,
  • quelle est la probabilité que la date où le million de pages vues sera atteint soit avant le 1er janvier
Une fois formalisée la question, reste à faire un peu d'économétrie.
  • modélisation économétrique
En faisant simple et rapide, afin de prendre en compte la corrélation forte avec la semaine précédente, et le fait que l'on s'intéresse à la somme cumulée, on peut considérer un modèle ARIMA
  • avec un retard d'ordre 7 pour les composantes moyennes mobiles et autorégressives,
  • avec une série intégrée à l'ordre 1,
L'ajustement se fait de la manière suivante,
> X=cumsum(base$nombre)
> model  <- arima(X,c(7 ,   # partie AR
+                                             1,    # partie I
+                                             7))   # partie MA
et la prévision se fait via
> forecast <- predict(model,200)

Ensuite, ce n'est qu'une représentation graphique,
> u=max(D)+1:200
> polygon(  c(u,rev(u)),  c(forecast$pred - 1.96*forecast$se,
+  rev(forecast$pred + 1.96*forecast$se)), col = "yellow", border=NA)
> lines(u,forecast$pred,col="blue",lwd=2)
> lines(u,forecast$pred- 1.96*forecast$se,col="blue",lty=2)
> lines(u,forecast$pred+ 1.96*forecast$se,col="blue",lty=2)
> abline(v=as.Date("01/01/2011","%d/%m/%Y"),col="red")
> abline(h=1000000,col="red")


Afin de répondre à la question posée, on peut étudier les différentes probabilités envisagées,
  • probabilité que le 1er janvier, j'ai atteint le million de pages vues
Dans un premier temps, on utilise la normalité des prédictions (en supposant une normalité du bruit) pour obtenir la loi de la prédiction à une date quelconque,
> k=which(u==as.Date("01/01/2011","%d/%m/%Y"))
> x=seq(800000,1100000,by=100)
> y=dnorm(x,mean=forecast$pred[k],sd=forecast$se[k])
> plot(x,y,type="l",lwd=2)
> x0=x[x>=1000000]
> polygon(  c(x0,rev(x0)), c(y[x>=1000000],rep(0,length(x0))), col = "yellow",border=NA)
> lines(x,y,type="l",lwd=2)
> abline(v=1000000)

Le cas de la probabilité d'atteindre plus d'un million de visiteur le 1er janvier est alors
> k=which(u==as.Date("01/01/2011","%d/%m/%Y"))
> 1-pnorm(1000000,mean=forecast$pred[k],sd=forecast$se[k])
[1] 0.2604821

  • problème dual, 
Dans un second temps, on peut envisager une autre approche, consistant à se demander quelle pourrait être la loi de la date du jour où le million de pages vues sera atteint,
> P=rep(NA,300)
> for(k in 1:300){
+ P[k]=1-pnorm(1000000,mean=forecast$pred[k],sd=forecast$se[k])}
> plot(max(D)+1:300,P,type="l",lwd=2)
> x=max(D)+1:300
> x0=x[x<=as.Date("01/01/2011","%d/%m/%Y")]
> polygon(  c(x0,rev(x0)), c(P[x<=as.Date("01/01/2011","%d/%m/%Y")],rep(0,length(x0))),  col = "yellow", border=NA)
> lines(max(D)+1:300,P,type="l",lwd=2)
> abline(v=as.Date("01/01/2011","%d/%m/%Y"))

La probabilité que cette date soit antérieure au 1er janvier est alors
> P[u==as.Date("01/01/2011","%d/%m/%Y")]
[1] 0.2604821
(ce qui, fort heureusement, correspond à la probabilité calculée par le problème primal). Bref, j'ai 1 chance sur 4 d'atteindre le million de pages vues avant la nouvelle année....
  • sensibilité à un changement de modèle
C'est bien joli tout ça, mais ces calculs sont largement soumis à la contrainte du choix de modèle que j'ai fait arbitrairement au début... on peut se demander si en changeant de modèle, les résultats changent sensiblement, ou pas. Au lieu de tenter un autre modèle ARIMA (voire SARIMA), j'ai préféré changer la série de référence... et me focaliser sur 2010 uniquement. D'un côté j'enlève les premières semaines où le niveau de fréquentation était très faible, de l'autre, je donne un poids très important aux vacances d'été, i.e. la période juin-août, pendant laquelle les internautes semblent moins sensibles à la modélisation économétrique.
Si l'on modélise la fréquentation pour 2010, seulement, on obtient

avec les distributions suivantes, tout d'abord pour la densité du nombre de visiteur atteint au 1er janvier 2011,

et pour la fonction de répartition de la date où sera atteint le million de pages vues,
soit une probabilité de l'ordre de 35%, dans les deux cas,
> P[u==as.Date("01/01/2011","%d/%m/%Y")]
[1] 0.3531648
> k=which(u==as.Date("01/01/2011","%d/%m/%Y"))
> 1-pnorm(1000000,mean=forecast2$pred[k],sd=forecast2$se[k])
[1] 0.3531648


On obtient des probabilités relativement proches avec les deux modèles, et j'aurais envie de croire que l'objectif est envisageable.... "challenge accepted" comme dirait mon ami Barney Stinson (oui, c'est mon ami). Reste juste à trouver des sujets qui attireront du monde en cette fin d'année...

Saturday, May 2 2009

Différencier (indéfiniment) les séries temporelles ?

Dans un mail, un étudiant qui finissait son projet de séries temporelles m'a posé une question simple et intéressante  (que je me permets de reprendre ici): "quand on cherche à stationnariser une série, on a souvent  besoin de différencier deux fois, et ça marche tout le temps".
Effectivement, on peut toujours différencier une série, on finira bien par tomber sur une série stationnaire. Je retraduirais cette question sous la forme suivante "pourquoi cherche-t-on toujours à stationnariser les séries ?" ou encore "est-ce gênant de différencier une série ?".
Pour la première question, la réponse est simple: les seules séries que l'on sache modéliser sont les séries stationnaires. Les autocorrélations se calculent sur des séries stationnaires, par exemple, et la non-stationnarité n'est pas une notion simple à définir (tout comme la non-indépendance entre variables aléatoires).....etc.
Pour la seconde question, la réponse est simple: si on différencie une série pour mieux la modéliser, si on souhaite par la suite faire de la prévision, il conviendra d'intégrer la série modélisée. Or intégrer une série fait exploser l'intervalle de confiance. Autant faire ça sur un exemple simple, avec la série suivante, que l'on commence par supposer stationnaire, et que l'on prédit sur une trentaine de valeurs, avec un intervalle de confiance.
XXXMais si l'on suppose la série intégrée à l'ordre 1, on différencie puis on modélise par un processus stationnaire (c'est l'idée qui sous-tend les processus ARIMA, à savoir un processus ARMA intégré: en différenciant la série, on devrait tomber sur un processus ARMA. On peut rattacher ça à la notion de racine unité). La prédiction (avec l'intervalle de confiance) est alors présenté sur la série initiale. Le fait d'intégrer les erreurs fait que la variance de la prédiction augmente avec l'horizon.
Et si l'on continue, et que l'on différencie de fois la série avant de la modéliser, et que l'on intègre deux fois les erreurs, l'intervalle de confiance devient immense !
XXXBref, il convient de ne pas différencier si ce n'est pas vraiment indispensable ! L'outil pour vérifier que l'on n'a pas trop différencier est la fonction d'autocovariance inverse (c'est à dire la fonction d'autocovariance d'une série dont la densité spectrale serait l'inverse de la densité spectrale initiale*). Cette fonction fait partie des fonctions de base sous SAS, mais pas sous R... J'ai cherché un peu, mais sans succès. L'argument dans les forums est que cette fonction est redondante avec la fonction d'autocorrélation partielle. Et en effet, si la finalité est simplement de détecter les ordres AR ou MA d'une série stationnaire, alors effectivement, les deux fonctions s'utilisent de manière équivalente. Mais la fonction d'autocovariance inverse apporte plus d'information, en particulier afin de détecter si l'on n'a pas surdifférencié la série.
La sortie suivante présente la fonction d'autocorrélation (à gauche) et la fonction d'autocorrélation inverse (à droite) sur une série.
XXXXOn note que la série semble intégrée (en pratique, des autocorrélations très fortes et persistantes très longtemps se traduit par une suspission d'intégration). Si l'on différencie la série, on obtient les autocorrélogrammes suivants,
XXXEt si l'on différencie une nouvelle fois, on obtient les graphiques suivants,
XXXXBref, l'autocorrélogramme de droite semble caractéristique d'une série intégrée: on dira alors que l'on a surdifférencié le modèle ("en intégrant la série, elle est toujours stationnaires"...). Je parle un peu de tout ça dans mon polycopié de Dauphine, tome 1.

* J'avais parlé dans un ancien billet (ici) de l'importance de la densité spectrale en série temporelles. L'idée est que si f est une densité spectrale, 1/f peut également l'être. Le lien entre la fonction d'autocorrélation et la densité spectrale est donné par des théorèmes de Kolmogorov ou Wiener. Pour plus de détails, il y a des compléments dans le chapitre 1.2 du polycopié que j'avais fait à Dauphine (ici). On peut d'ailleurs noter que la fonction d'autocorrélation partielle d'un processus ARMA(p,q) est la fonction d'autocorrélation d'un processus ARMA(q,p) obtenu en permutant les deux polynômes d'opérateurs retard.