Séries temporelles sous R - Simon Bussy

1 S ries temporelles sous R. V. Lefieux EXEMPLES DE SERIES temporelles . Les s ries apparaissent dans l'ordre du cours. S rie sunspot : nombre annuel de t ches solaires de 1790 1970. plot( ,xlab="t",ylab="Sunspots"). 150. Sunspots 100. 50. 0. 1700 1750 1800 1850 1900 1950. t Bruit blanc gaussien de loi N (0, 32 ). Pour les simulations effectu es dans ce document, on fixe arbitrairement la racine (seed) . 1789. (1789). plot(ts(rnorm(100,sd=3),start=1,end=100) ,xlab="t",ylab="Bruit blanc gaussien de variance 9"). abline(h=0). 1. 6. Bruit blanc gaussien de variance 9. 4. 2. 0. 2. 4. 6. 0 20 40 60 80 100. t S rie uspop : population des Etats-Unis, en millions, de 1790 1990 (Pas de temps d cennal).

2 Plot(uspop,xlab="t",ylab="Uspop"). 200. 150. Uspop 100. 50. 0. 1800 1850 1900 1950. t 2. S rie airpass : nombre mensuel de passagers a riens, en milliers, de janvier 1949. d cembre 1960. S rie Brute plot(AirPassengers,xlab="t",ylab="Airpas s"). 600. 500. 400. Airpass 300. 200. 100. 1950 1952 1954 1956 1958 1960. t Logarithme de la s rie airpass plot(log(AirPassengers),xlab="t",ylab="A irpass"). 3. Airpass 1950 1952 1954 1956 1958 1960. t S rie beer : production mensuelle de bi re en Australie, en m galitres, de janvier 1956 f vrier 1991. beer= ("../ ",header=F,dec=".",sep=","). beer=ts(beer[,2],start=1956,freq=12). plot(beer,xlab="t",ylab="Beer").

3 200. 150. Beer 100. 1960 1970 1980 1990. t 4. S rie lynx : nombre annuel de lynx captur s au Canada, de 1821 1934. plot(lynx,xlab="t",ylab="Lynx"). 6000. 4000. Lynx 2000. 0. 1820 1840 1860 1880 1900 1920. t Sauf mention contraire, on travaille dans la suite sur la s rie temporelle airpass. x=AirPassengers y=log(x). CHAPITRE 1 : DECOMPOSITION SAISONNIERE. D composition saisonni re l'aide de la r gression lin aire Cr ation des bases tendancielle et saisonni re t=1:144. for (i in 1:12). {. su=rep(0,times=12). su[i]=1. s=rep(su,times=12). assign(paste("s",i,sep=""),s). }. 5. R gression lin aire reg=lm(y~t+s1+s2+s3+s4+s5+s6+s7+s8+s9+s1 0+s11+s12-1).

4 Summary(reg). ##. ## Call: ## lm(formula = y ~ t + s1 + s2 + s3 + s4 + s5 + s6 + s7 + s8 +. ## s9 + s10 + s11 + s12 - 1). ##. ## Residuals: ## Min 1Q Median 3Q Max ## ##. ## Coefficients: ## Estimate Std. Error t value Pr(>|t|). ## t <2e-16 **. ## s1 <2e-16 **. ## s2 <2e-16 **. ## s3 <2e-16 **. ## s4 <2e-16 **. ## s5 <2e-16 **. ## s6 <2e-16 **. ## s7 <2e-16 **. ## s8 <2e-16 **. ## s9 <2e-16 **. ## s10 <2e-16 **. ## s11 <2e-16 **. ## s12 <2e-16 **. ## --- ## Signif. codes: 0 '**' '**' '*' '.' ' ' 1. ##. ## Residual standard error: on 131 degrees of freedom ## Multiple R-squared: , Adjusted R-squared: ## F-statistic: +04 on 13 and 131 DF, p-value: < reg$coefficients ## t s1 s2 s3 s4 s5 s6.

5 ## ## s7 s8 s9 s10 s11 s12. ## a=mean(reg$coefficients[2:13]). b=reg$coefficients[1]. c=reg$coefficients[2:13]-mean(reg$coeffi cients[2:13]). Calcul de la s rie corrig e des variations saisonni res 6. y_cvs=y-(c[1]*s1+c[2]*s2+c[3]*s3+c[4]*s4 +c[5]*s5+c[6]*s6+c[7]*s7+c[8]*s8+c[9]*s9 +c[10]*s10+c[11]*s11+c[1. x_cvs=exp(y_cvs). (x,x_cvs,xlab="t",ylab="Airpass",col=c(1 ,2),lwd=c(1,2)). legend("topleft",legend=c("X","X_CVS"),c ol=c(1,2),lwd=c(1,2)). 600. X. X_CVS. 500. 400. Airpass 300. 200. 100. 1950 1952 1954 1956 1958 1960. t D composition saisonni re l'aide des moyennes mobiles On utilise les moyennes mobiles M2 12 et M3 dans la premi re tape de l'algorithme X11.)]

6 M2_12=function(x){. y=(1/12)*filter(x,c( ,rep(1,times=11), )). return(y). }. m3=function(x){. y=(1/3)*filter(x,rep(1,times=3)). return(y). }. H. On utiliserait les moyennes mobiles M13 et M5 dans la deuxi me tape de l'algorithme X11. m13h=function(x){. y=(1/16796)*filter(x,c(-325,-468,0,1100, 2475,3600,4032,3600,2475,1100,0,-468,-32 5)). return(y). }. m5=function(x){. y=(1/5)*filter(x,rep(1,times=5)). return(y). }. 7. Le premier jeu d'estimation donne : t1=m2_12(y). sig1=y-t1. s1=m3(m3(sig1)). shat1=s1-m2_12(s1). ycvs1=y-shat1. xcvs1=exp(ycvs1). (x,xcvs1,col=c(1,2),lwd=c(1,2)). legend("topleft",legend=c("X","X_CVS"),c ol=c(1,2),lwd=c(1,2)). 600.

7 X. X_CVS. 500. 400. 300. 200. 100. 1950 1952 1954 1956 1958 1960. Time Il faudrait effectuer les 4 tapes suivantes et compl ter les donn es limin es par des moyennes mobiles asym triques. Notons qu'il est galement possible d'utiliser la librairie (compl te) X12. D composition saisonni re l'aide de la fonction decompose (x,type="multiplicative"). $figure ## [1] ## [8] plot( ). 8. observed Decomposition of multiplicative time series 400. 100. 450. trend 300. random seasonal 1950 1952 1954 1956 1958 1960. Time CHAPITRE 1 : LISSAGE EXPONENTIEL. On utilise la librairie forecast. library(forecast). Lissage exponentiel simple les=ets(y,model="ANN").

8 (les,12). plot( ). 9. Forecasts from ETS(A,N,N). 1950 1952 1954 1956 1958 1960 1962. Lissage exponentiel double led=ets(x,model="MMN"). (led,12). plot( ). Forecasts from ETS(M,M,N). 800. 600. 400. 200. 1950 1952 1954 1956 1958 1960 1962. 10. M thode de Holt-Winters hw=ets(x,model="MMM"). (hw,12). plot( ). Forecasts from ETS(M,Md,M). 700. 500. 300. 100. 1950 1952 1954 1956 1958 1960 1962. CHAPITRE 2. Blancheur On utilise la librairie caschrono. library(caschrono).. On obtient sur un bruit blanc gaussien de loi N 0, 32 : (1789). (rnorm(100,sd=3),start=1,end=100). ( ,nlag=c(5,10,20),type="Ljung-Box",decim= 5). ## Retard p-value ## [1,] 5 ## [2,] 10 ## [3,] 20 On peut galement visualiser ses autocorr logrammes empiriques simple et partiel.

9 11. acf2y( , ). Time series: ACF. PACF. 5 10 15 20. Lag ## LAG ACF1 PACF. ## [1,] 1 ## [2,] 2 ## [3,] 3 ## [4,] 4 ## [5,] 5 ## [6,] 6 ## [7,] 7 ## [8,] 8 ## [9,] 9 ## [10,] 10 ## [11,] 11 ## [12,] 12 ## [13,] 13 ## [14,] 14 ## [15,] 15 ## [16,] 16 ## [17,] 17 ## [18,] 18 ## [19,] 19 ## [20,] 20 On obtient pour un processus AR(1), Xt = 1 + t o Var (Xt ) = 32 : 12. (1789). (n=100,list(ar= ),sd=3). ( ,nlag=c(1,5,10,20),type="Ljung-Box",deci m=5). ## Retard p-value ## [1,] 1 0. ## [2,] 5 0. ## [3,] 10 0. ## [4,] 20 0. P riodogramme On utilise la librairie TSA. library(TSA). On obtient pour la s rie lynx : (x,ylab="Periodogramme"). 8e+05. Periodogramme 4e+05.

10 0e+00. Frequency On peut ensuite d terminer pour quelle fr quence le p riodogramme est maximal, etc. $freq[ ( ( $spec))]*114. ## [1] 13. CHAPITRE 4 : PROCESSUS AR, MA & ARMA. On utilise la librairie caschrono. library(caschrono). Autocorr logrammes simple et partiel d'un processus AR. On obtient pour un processus AR(1), Xt = 1 + t o Var (Xt ) = 32 : (1789). (n=100,list(ar= ),sd=3). plot( ,xlab="t",ylab="X",main="AR(1):phi1= ; cart-type=3"). abline(h=0,lty=2). AR(1):phi1= ; cart type=3. 5. X. 0. 5. 0 20 40 60 80 100. t acf2y( , ). 14. Time series: ACF. PACF. 5 10 15 20. Lag ## LAG ACF1 PACF. ## [1,] 1 ## [2,] 2 ## [3,] 3 ## [4,] 4 ## [5,] 5 ## [6,] 6 ## [7,] 7 ## [8,] 8 ## [9,] 9 ## [10,] 10 ## [11,] 11 ## [12,] 12 ## [13,] 13 ## [14,] 14 ## [15,] 15 ## [16,] 16 ## [17,] 17 ## [18,] 18 ## [19,] 19 ## [20,] 20 On obtient pour un processus AR(1), Xt = 1 + t o Var (Xt ) = 32 : (1789).