# Inference for ARMA(p,q) Time Series

As we mentioned in our previous post, as soon as we have a moving average part, inference becomes more complicated. Again, to illustrate, we do not need a two general model. Consider, here, some  process,

where  is some white noise, and assume further that .

> theta=.7
> phi=.5
> n=1000
> Z=rep(0,n)
> set.seed(1)
> e=rnorm(n)
> for(t in 2:n) Z[t]=phi*Z[t-1]+e[t]+theta*e[t-1]
> Z=Z[800:1000]
> plot(Z,type="l") • A two step procedure

To start with something simple, assume that we did miss the moving average component,

The estimator of  – by least squares – is not longer consistent. But still. We can still compute it

> base=data.frame(Y=Z[2:n],X=Z[1:(n-1)])
> regression=lm(Y~0+X,data=base)
> summary(regression)

Call:
lm(formula = Y ~ 0 + X, data = base)

Residuals:
Min      1Q  Median      3Q     Max
-3.2445 -0.7909  0.0626  0.9707  3.0685

Coefficients:
Estimate Std. Error t value Pr(>|t|)
X  0.69571    0.05101   13.64   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.225 on 199 degrees of freedom
(799 observations deleted due to missingness)
Multiple R-squared:  0.4832,	Adjusted R-squared:  0.4806
F-statistic:   186 on 1 and 199 DF,  p-value: < 2.2e-16

and then, we cancompute the autocorrelation of the noise,

> n=200
> cor(residuals(regression)[2:n],residuals(regression)[1:(n-1)])
 0.2663076

or more formally, use Durbin-Watson estimator, to get autocorrelation of the noise (and some significance test)

> library(car)
> durbinWatsonTest(regression)
lag Autocorrelation D-W Statistic p-value
1       0.2656555       1.46323       0
Alternative hypothesis: rho != 0

The point, here, is that we would like to assume that

meaning that  should be some  process. And

i.e.  is a root of this quadratic problem,

> polyroot(c(1,-1/cor(residuals(regression)[2:n],residuals(regression)[1:(n-1)]),1))
 0.2884681+0i 3.4665883+0i

Here, we do have two positive roots. I would go for the one smaller than one, in order to be able to invert the polynomial, if necessary…

•  Use of the empirical autocorrelation function

An alternative might be to use properties of the autocorrelation function,

and

Again, we have a set of three equations, with three unknown parameters. Numerically, it is possible to find some roots. If we run the code, we get

> v=c(as.numeric(acf(Z)$acf[2:3]),1)*var(Z) > library(rootSolve) > seteq=function(x){ + F1=v-x^2*(x^2+2*x*x+1)/(1-x^2) + F2=v-(x*v+x*x^2) + F3=v-x*v + return(c(F1,F2,F3))} > multiroot(f=seteq,start=c(.1,.1,1))$root
  3.643734 -3.188145  1.427759

$f.root  1.371170e-11 -3.714573e-11 0.000000e+00$iter
 8

$estim.precis  1.695248e-11 Here, we have a situation… • Use of least square techniques We can use, here, the algorithm described in the context of processes. > V=function(p){ + phi=p + theta=p + u=rep(0,length(Z)) + for(t in 2:length(Z)) u[t]=Z[t]-phi*Z[t-1]-theta*u[t-1] + return(sum(u^2)) + } > p=optim(par=c(.1,.1),V)$par
 0.3637783 0.7773845
> coef=c(p,sqrt(V(p)/(length(Z))))

which is not so bad. Actually, if we run that procecure on 1,000 samples, we get the following output • Use of maximum likelihood techniques

Last, but not least, one more time, we can use (global) maximum likelihood techniques, since the process is a Gaussian process (all finite dimensional vector will have a joint Gaussian distribution) if we assume that the noise is Gaussian.

> library(mnormt)
> GlobalLogLik=function(A,TS){
+ n=length(TS)
+ phi=A;  theta=A
+ sigma=A
+ SIG=matrix(0,n,n)
+ rho=rep(0,n)
+ rho=sigma^2*(theta^2+2*phi*theta+1)/(1-phi^2)
+ rho=phi*rho+theta*sigma^2
+ for(h in 3:n) rho[h]=phi*rho[h-1]
+ for(i in 1:n){for(j in 1:n){
+ SIG[i,j]=rho[abs(i-j)+1]}}
+ return(dmnorm(TS,rep(0,n),SIG,log=TRUE))}
> LogL=function(A) -GlobalLogLik(A,TS=Z)
> optim(c(.1,.1,1),LogL)$par  0.3890991 0.7672036 1.0731340 It works fine, one more time. But maybe we got lucky here. We’ve seen in the post on autoregressive time series that the algorithm might fell if the time series is not stationary. In order to avoid such problems, we can consider a constraint optimization problem, where we simply recall that , > U=matrix(c(1,-1,0,0,0,0),2,3) > C=-c(.999,.999) > constrOptim(c(.1,.1,1),LogL,grad=NULL,ui=U,ci=C)$par
 0.3890991 0.7672036 1.0731340

$value  300.1956$counts
118       NA

$convergence  0$message
NULL

$outer.iterations  2$barrier.value
 -1.536358e-05

If we run that algorithm 1,000 times, on simulated time series (with the same parameters), we get # Introduction aux processus 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,

# Prévision des vente de voitures au Québec

Après discussion, la date limite pour me renvoyer le second devoir est fixée à vendredi 5 octobre midi, par courriel. Je veux un fichier pdf par équipe, avec une page de garde indiquant le nom des membres du groupe, et le mot clé retenu.

Maintenant, je vais en profiter pour mettre en ligne le code tapé de matin, en cours. On avait travaillé sur la modélisation de la série des ventesde voiture, au Québec. La série est en ligne,

X=read.table(
"http://freakonometrics.blog.free.fr/public/
data/car-sales-quebec.csv",
Xt=ts(X[,2],start=c(1960,1),frequency=12)

On peut regarder l’évolution de cette série temporelles (car les dessins, et la visualisation, c’est important),

plot(Xt) On note une tendance linéaire (ou qu’on pourrait supposer linéaire), que l’on va estimer,

X=as.numeric(Xt)
temps=1:length(X)
base=data.frame(temps,X)
reg=lm(X~temps)

et que l’on pourra représenter

plot(temps,X,type="l",xlim=c(0,120))
abline(reg,col="red")

Si on veut une série stationnaire (et c’est effectivement ce que l’on cherche), on va retrancher cette tendance à notre série, et travailler sur la série résiduelle. La série résiduelle est ici

Y=X-predict(reg)
plot(temps,Y,type="l") que l’on devrait pouvoir supposer stationnaire. Classiquement, on regarde les autocorrélations,

acf(Y,36,lwd=4) où on repère un beau cycle annuel, mais on va surtout utiliser les autocorrélations partielles pour identifier l’ordre de la composante autorégressive.

pacf(Y,36,lwd=4) La 12ème est non-nulle, et tous les autres ensuite sont significativement nulles (ou presque). On va donc tenter un .

> fit.ar12=arima(Y,order=c(12,0,0))
> fit.ar12
Series: Y
ARIMA(12,0,0) with non-zero mean

Coefficients:
ar1     ar2      ar3      ar4     ar5      ar6     ar7
0.1975  0.0832  -0.1062  -0.1212  0.1437  -0.1051  0.0319
s.e.  0.0838  0.0809   0.0826   0.0843  0.0850   0.0833  0.0854
ar9     ar10    ar11    ar12  intercept
-0.0332  -0.0616  0.2635  0.4913  -148.3180
s.e.   0.0853   0.0840  0.0840  0.0841   384.5095

sigma^2 estimated as 2177974:  log likelihood=-946.75
AIC=1921.51   AICc=1926.03   BIC=1959.06

La douzième composante  semble significative, si on repense au test de Student,

u=seq(-6,6,by=.1)
plot(u,dnorm(u),type="l")
points(0.4913/0.0841,0,pch=19,col="red")

En revanche, la huitième  ne l’est pas

points(-0.1018/0.0847,0,pch=19,col="blue") (voilà pour le petit retour sur la lecture des tests). Vérifions maintenant que le bruit résiduel est bien blanc… Si on regarde les autocorrélations

acf(residuals(fit.ar12),36,lwd=4) Mais ce n’est pas un test de bruit blanc. En particulier, si une autocorrélation semblerait significative, on pourrait accepter que – globalement – les autocorrélations soient très proches de 0. On va alors faire des tests de Box-Pierce

> Box.test(residuals(fit.ar12),lag=12,
+  type='Box-Pierce')

Box-Pierce test

data:  residuals(fit.ar12)
X-squared = 7.7883, df = 12, p-value = 0.8014

Pour l’interpétation du test, c’est toujours pareil: on a une statistique de l’ordre de 7, et la loi sous-jacente (la somme des carrés des 12 premières autocorrélations du bruit) est une loi du chi-deux, à 12 degrés de liberté,

u=seq(0,30,by=.1)
plot(u,dchisq(u,df=12),type="l")
points(7.7883,0,pch=19,
col="red") On accepte l’hypothèse que les 12 premières autocorrélations soient nulles. Si on va un peu plus loin,

> Box.test(residuals(fit.ar12),lag=18,
+  type='Box-Pierce')

Box-Pierce test

data:  residuals(fit.ar12)
X-squared = 20.3861, df = 18, p-value = 0.3115

on a la même conclusion, Pour rappels, on peut retrouver à la main la p-value,

> 1-pchisq(20.3861,df=18)
 0.3115071

Faisons notre petite dessin de toutes les p-value, pour s’assurer que le bruit est bien un bruit blanc,

 BP=function(h) Box.test(residuals(fit.ar12),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") Cette fois c’est bon, on tient notre premier modèle, un . Si on veut aller plus loin, on peut regarder ce que donnerait de la sélection automatique, > library(caschrono) > armaselect(Y,nbmod=5) p q sbc [1,] 14 1 1635.214 [2,] 12 1 1635.645 [3,] 15 1 1638.178 [4,] 12 3 1638.297 [5,] 12 4 1639.232 On peut être tenté d’aller voir ce qui se passerait en rajoutant une composante moyenne mobile à notre modèle précédant, > fit.arma12.1=arima(Y,order=c(12,0,1)) > fit.arma12.1 Series: Y ARIMA(12,0,1) with non-zero mean Coefficients: ar1 ar2 ar3 ar4 ar5 ar6 ar7 0.0301 0.1558 -0.0941 -0.1461 0.1063 -0.0688 -0.002 s.e. 0.1235 0.0854 0.0757 0.0784 0.0807 0.0774 0.080 ar9 ar10 ar11 ar12 ma1 intercept -0.0646 -0.0798 0.2538 0.5786 0.2231 -131.3495 s.e. 0.0802 0.0766 0.0751 0.0861 0.1393 368.8156 sigma^2 estimated as 2127759: log likelihood=-945.65 AIC=1921.31 AICc=1926.52 BIC=1961.54 Le nouveau coefficient, est à peine significatif. Éventuellement avec un seuil à 10%… Pourquoi pas? Si on regarde les résidus, sans grande surprise, on a toujours un bruit blanc, encore plus blanc qu’auparavant (en violet sur le dessin ci-dessous) BP=function(h) Box.test( residuals(fit.arma12.1),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")

BP=function(h) Box.test(residuals(fit.ar12),lag=h,
type='Ljung-Box')$p.value lines(1:24,Vectorize(BP)(1:24),col="purple", type="b") On peut aussi aller voir parmi les modèles proposés, en particulier le modèle . > fit.ar14= + arima(Y,order=c(14,0,0),method="CSS") > fit.ar14 Series: Y ARIMA(14,0,0) with non-zero mean Coefficients: ar1 ar2 ar3 ar4 ar5 ar6 ar7 0.2495 0.2105 -0.0584 -0.1569 0.1282 -0.1152 0.0268 s.e. 0.0956 0.0972 0.0854 0.0830 0.0838 0.0840 0.0847 ar9 ar10 ar11 ar12 ar13 ar14 intercept -0.0327 -0.1116 0.2649 0.5887 -0.1575 -0.1572 80.5 s.e. 0.0855 0.0851 0.0853 0.0886 0.1031 0.0999 338.9 sigma^2 estimated as 2218612: part log likelihood=-942.31 C’est un peu tiré par les cheveux, mais on pourrait accepter l’hypothèse que soit significativement non-nuls. Mais on est encore limite… Allez, on l’accepte aussi dans notre gang de modèle. On a finalement trois modèles. Si on fait un peu de backtesting, sur les 12 derniers mois, T=length(Y) backtest=12 subY=Y[-((T-backtest+1):T)] subtemps=1:(T-backtest) plot(temps,Y,type="l") lines(subtemps,subY,lwd=4) fit.ar12.s=arima(subY, order=c(p=12,d=0,q=0),method="CSS") fit.arma12.1.s=arima(subY, order=c(p=12,d=0,q=1),method="CSS") fit.ar14.s=arima(subY, order=c(p=14,d=0,q=0),method="CSS") p.ar12=predict(fit.ar12.s,12) pred.ar12=as.numeric(p.ar12$pred)
p.arma12.1=predict(fit.arma12.1.s,12)
pred.arma12.1=as.numeric(p.arma12.1$pred) p.ar14=predict(fit.ar14.s,12) pred.ar14=as.numeric(p.ar14$pred)

on obtient les prévisions suivantes (avec les valeurs observées dans la première colonne)

> (M=cbind(observé=Y[((T-backtest+1):T)],
+  modèle1=pred.ar12,
+  modèle2=pred.arma12.1,
+  modèle3=pred.ar14))
observé    modèle1    modèle2    modèle3
97  -4836.2174 -5689.3331 -5885.4486 -6364.2471
98  -3876.4199 -4274.0391 -4287.2193 -4773.8116
99   1930.3776  1817.8411  2127.9915  2290.1460
100  3435.1751  4089.3598  3736.1110  4039.4150
101  7727.9726  6998.9829  7391.6694  7281.4797
102  2631.7701  3456.8819  3397.5478  4230.5324
103  -509.4324 -2128.6315 -2268.9672 -2258.7216
104 -1892.6349 -3877.7609 -3694.9409 -3620.4798
105 -4310.8374 -3384.0905 -3430.4090 -2881.4942
106  2564.9600  -504.6883  -242.5018   183.2891
107 -1678.2425 -1540.9904 -1607.5996  -855.7677
108 -4362.4450 -3927.4772 -3928.0626 -3718.3922
>  sum((M[,1]-M[,2])^2)
 19590931
>  sum((M[,1]-M[,3])^2)
 17293716
>  sum((M[,1]-M[,4])^2)
 21242230

I.e. on aurait envie de retenir le second modèle, . On va maintenant l’utiliser pour faire un peu de prévision,

library(forecast)
fit.arma12.1=
arima(Y,order=c(12,0,1))
fit.arma12.1
PREDARMA=forecast(fit.arma12.1,12)
plot(Y,xlim=c(1,120),type="l")
temps=T+1:12
polygon(c(temps,rev(temps)),c(PREDARMA$lower[,2], rev(PREDARMA$upper[,2])),col="yellow",border=NA)
polygon(c(temps,rev(temps)),c(PREDARMA$lower[,1], rev(PREDARMA$upper[,1])),col="orange",border=NA)
lines(temps,PREDARMA$mean,col="red") Cela dit, on peut aussi aller beaucoup plus loin dans la prévision, PREDARMA=forecast(fit.arma12.1,120) plot(Y,xlim=c(1,210),type="l") temps=T+1:120 polygon(c(temps,rev(temps)),c(PREDARMA$lower[,2],
rev(PREDARMA$upper[,2])),col="yellow",border=NA) polygon(c(temps,rev(temps)),c(PREDARMA$lower[,1],
rev(PREDARMA$upper[,1])),col="orange",border=NA) lines(temps,PREDARMA$mean,col="red") Bon, on y est presque, car on a modélisé , la série obtenue en enlevant la tendance linéaire. Pour remonter sur , on va rajouter la tendance à la prévision faite auparavant,

X=as.numeric(Xt)
temps=1:length(X)
plot(temps,X,type="l",xlim=c(0,210),
ylim=c(5000,30000))
base=data.frame(temps,X)
reg=lm(X~temps)
abline(reg,col="red")
PREDTENDANCE=predict(reg,newdata=
data.frame(temps=T+1:120))
temps=T+1:120
polygon(c(temps,rev(temps)),c(PREDTENDANCE+
PREDARMA$lower[,2],rev(PREDTENDANCE+PREDARMA$upper[,2])),
col="yellow",border=NA)
polygon(c(temps,rev(temps)),c(PREDTENDANCE+
PREDARMA$lower[,1],rev(PREDTENDANCE+PREDARMA$upper[,1])),
col="orange",border=NA)
lines(temps,PREDTENDANCE+PREDARMA$mean,col="red") # 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 , 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 et dans un modèle 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
 8

$ma.ma  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 ( ) et en abscisse les ordres moyenne mobile ( ).

> library(RColorBrewer)
> CL=brewer.pal(6, "RdBu")
> ceacf=matrix(as.numeric(cut(EACF$eacf, +1,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 ( et 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  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 deCryer & 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 =12, histoire de voir si on ne pourrait pas avoir un modèle simple, et facilement interprétable en plus… 1. -3):3)/3,labels=1:6 [] # Inference and autoregressive processes Consider a (stationary) autoregressive process, say of order 2, for some white noise with variance . Here is a code to generate such a process, > phi1=.5 > phi2=-.4 > sigma=1.5 > set.seed(1) > n=240 > WN=rnorm(n,sd=sigma) > Z=rep(NA,n) > Z[1:2]=rnorm(2,0,1) > for(t in 3:n){Z[t]=phi1*Z[t-1]+phi2*Z[t-2]+WN[t]} Here, we have to estimate two sets of parameters: the autoregressive coefficients , and the variance of the innovation process . There are (at least) three techniques to estimate those parameters. • using least square regression A natural idea is to see here a regression model, and thus, if we consider a matrix formulation, Here we can run (conditional) ordinary least squares estimation, > base=data.frame(Y=Z[3:n],X1=Z[2:(n-1)],X2=Z[1:(n-2)]) > regression=lm(Y~0+X1+X2,data=base) > summary(regression) Call: lm(formula = Y ~ 0 + X1 + X2, data = base) Residuals: Min 1Q Median 3Q Max -4.3491 -0.8890 -0.0762 0.9601 3.6105 Coefficients: Estimate Std. Error t value Pr(>|t|) X1 0.45107 0.05924 7.615 6.34e-13 *** X2 -0.41454 0.05924 -6.998 2.67e-11 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.449 on 236 degrees of freedom Multiple R-squared: 0.2561, Adjusted R-squared: 0.2497 F-statistic: 40.61 on 2 and 236 DF, p-value: 6.949e-16 > regression$coefficients
X1         X2
0.4510703 -0.4145365
> summary(regression)$sigma  1.449276 • using Yule-Walker equations As we’ve seen in class, we can easily get the following equations for the autocovariance functions, which can also be written So we just have to solve a simple linear system of equations. Note that if we divide by the variance , those equations can be written in terms of the autocorrelation functions The code is the following > rho1=cor(Z[1:(n-1)],Z[2:n]) > rho2=cor(Z[1:(n-2)],Z[3:n]) > A=matrix(c(1,rho1,rho1,1),2,2) > b=matrix(c(rho1,rho2),2,1) > (PHI=solve(A,b)) [,1] [1,] 0.4517579 [2,] -0.4155920 Now, we need to extract the estimated innovation process, from this set of parameters (note that it could be possible to include the variance term in Yule-Walker equations, to get a three dimensional linear equation) > estWN=base$Y-(PHI*base$X1+PHI*base$X2)
> sd(estWN)
 1.445706

This estimator is probably not the best one (we can take into account that we’ve lost two degrees of freedom), but as a starting point, let us consider this one.

• using (conditional) likelihood estimators

Finally, we can assume some distribution for the innovation process. Thestandard model is a Gaussian model, i.e. In that case, the conditional log likelihood (conditional since we set the first two observations here) is

> CondLogLik=function(A,TS){
+ phi1=A;  phi2=A
+ sigma=A	; L=0
+ for(t in 3:length(TS)){
+ L=L+dnorm(TS[t],mean=phi1*TS[t-1]+
+ phi2*TS[t-2],sd=sigma,log=TRUE)}
+ return(-L)}

Now, we can run standard optimization procedures,

> LogL=function(A) CondLogLik(A,TS=Z)
> optim(c(0,0,1),LogL)
$par  0.4509685 -0.4144938 1.4430930$value
 425.0164

$counts function gradient 88 NA$convergence
 0

$message NULL Here, our three estimators are rather close. Actually, if we generate 1,000 time series (of size 240), those are the Box-plots of our three estimators, for the first order autoregressive coefficient for the second one, and finally for the standard deviation of the innovation process All those estimators behave nicely, and are rather close. Note that they all might be biased, but they are consistent (see Davidson and MacKinnon for instance, in their book, for more details). # 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.931767 mais comme le but est de comprendre ce qui se passe, autant faire les choses calmement, > X=rep(NA,1010) > X=0 > for(t in 2:1010){X[t]=4/3+X[t-1]/3+rnorm(1)} > X=X[-(1:10)] > mean(X)  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=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:loglikelihood = -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,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 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 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…

# The return period of the 2003 heat wave

The paper on “on the return period of the 2003 heat wave” has been published online here, before Christmas. It should be published soon by Climatic Change, http://link.springer.com/article/...

Extremal events are difficult to model since it is difficult to characterize formally those events. The 2003 heat wave in Europe was not characterized by very high temperatures, but mainly the fact that night temperature were no cool enough for a long period of time. Hence, simulation of several models (either with heavy tailed noise or long range dependence) yield different estimations for the return period of that extremal event.

To go further on the impact on mortality (which was not the aim of the paper), there is a paper in Nature (here). 