Identification of ARMA processes

Last week (in the MAT8181 course) in order to identify the orders of an ARMA process, we’ve seen the eacf method, and I mentioned the scan method, introduced in Tsay and Tiao (1985). The code below – to produce the output of the scan procedure – has been adapted from an old code by Steve Chen (where I included a visualization of the p-values, with the following colors)

The procedure was described in the course, last Thursday,

```arma.scan=function(z,ar.max=15,ma.max=15,alpha=0.01)
{
ym=function(z,t,m){return(z[t:(t-m)])}
n=length(z)
z=z - mean(z)
cmax=ma.max + 1
rmax=ar.max + 1
corref=matrix(0,nrow=rmax,ncol=cmax)
cmj.table=matrix(0,nrow=rmax,ncol=cmax)
pv=matrix(0,nrow=rmax,ncol=cmax)
mark=matrix(rep("X",(rmax)*(cmax)),nrow=rmax,ncol=cmax)
Rnames=paste("AR",0:(ar.max),sep="-")
Cnames=paste("MA",0:(ma.max),sep="-")
rownames(corref)=Rnames
colnames(corref)=Cnames
rownames(cmj.table)=Rnames
colnames(cmj.table)=Cnames
rownames(pv)=Rnames
colnames(pv)=Cnames
rownames(mark)=Rnames
colnames(mark)=Cnames
for (m in 0:ar.max)
{
m1=m+1
for (j in 0:ma.max)
{
j1=j+1
if (m == 0 && j != 0)
{
racf=acf(z,plot=FALSE)\$acf[1:(j+1)]
lamb=racf[j+1]^2
corref[m1,j]=round(lamb,4)
dmj=1 + 2*sum(racf[1:j]^2)
cmj=-1*(n-m-j)*log(1.0 - lamb/dmj)
pvalue =pchisq(cmj,1,lower.tail=FALSE)
pv[m1,j]=round(pvalue,4)
cmj.table[m1,j]=round(cmj,4)
mark[m1,j]=ifelse(pvalue > alpha,"O","X")
}
else if (m != 0 && j == 0)
{
racf=pacf(z,plot=FALSE)\$acf[1:(m+1)]
lamb=racf[m+1]^2
corref[m1,j1]=round(lamb,4)
dmj = 1
cmj=-1*(n-m-j)*log(1.0 - lamb/dmj)
pvalue =pchisq(cmj,1,lower.tail=FALSE)
pv[m1,j1]=round(pvalue,4)
cmj.table[m1,j1]=round(cmj,4)
mark[m1,j1]=ifelse(pvalue > alpha,"O","X")
}
else
{
mat1=matrix(0,nrow=m1,ncol=m1)
mat2=matrix(0,nrow=m1,ncol=m1)
mat3=matrix(0,nrow=m1,ncol=m1)
mat4=matrix(0,nrow=m1,ncol=m1)
for (t in (j+m+2):n)
{
tj1=t-j-1
ym1=ym(z,tj1,m)
ym2=ym(z,t,m)

mat1=mat1 + as.matrix(ym1)%*%ym1
mat2=mat2 + as.matrix(ym1)%*%ym2
mat3=mat3 + as.matrix(ym2)%*%ym2
mat4=mat4 + as.matrix(ym2)%*%ym1
}
b1=solve(mat1)%*%mat2
b2=solve(mat3)%*%mat4
A=b2%*%b1
eig <-eigen(A)
eig.val <-eig\$values
eig.val=Re(eig.val)
eig.len=length(eig.val)
eig.vector=eig\$vectors
lamb=min(eig.val)
eig.vector0=eig.vector[,which.min(eig.val)]
eig.vector0 = eig.vector0/eig.vector0[1]
resid=(1:n)*0
for (t in (j+m+1):n)
{
z0=z[seq(t,t-m,-1)]
resid[t]=sum(z0 * eig.vector0)
}
jm1=j + m + 1
rx=Re(resid[jm1:n])
racf=acf(rx,plot=FALSE)\$acf[1:j]
dmj=1 + 2*sum(racf^2)
cmj=-1*(n-m-j)*log(1.0 - lamb/dmj)
pvalue =pchisq(cmj,df=1,lower.tail=FALSE)
corref[m1,j1]=round(lamb,4)
pv[m1,j1]=round(pvalue,4)
cmj.table[m1,j1]=round(cmj,4)
mark[m1,j1]=ifelse(pvalue > alpha,"O","X")
}
}
}

cat("\n\nSCAN: Smallest CANonical Correlation Method for ARIMA(p,d,q)\n\n")
cat("Estimates of Squared Canonical Correlation \n\n")
print(corref)
cat("\n\nC(m,j)\n\n")
print(cmj.table)
cat("\n\nChi-Square(1) Test p-value\n\n")
print(pv)
cat("\nSCAN Matrix \n\n")
print(mark)

plot(0:1,0:1,col="white",xlim=c(0,nrow(pv)-1),ylim=c(0,ncol(pv)-1),axes=FALSE,xlab="AR",ylab="MA")
axis(1); axis(2)
library(RColorBrewer)
CL=brewer.pal(6, "RdBu")[c(1,2,3,5)]
cpv=matrix(as.numeric(cut(as.vector(pv),c(-1,.01,.05,.1,2))),nrow(pv),ncol(pv))
for(i in 1:nrow(pv)){
for(j in 1:ncol(pv)){
polygon(c(i-1,i-1,i,i)-.5,c(j-1,j,j,j-1)-.5,
col=CL[cpv[i,j]])
}}
}```

Consider the following simulated time series,

```> s=arima.sim(n=200,model=list(ar=c(0,0,0,.4,0,0,0,.5),ma=c(0,0,1)))
> plot(s,type="l")```

The output is here

```> arma.scan(s,6,6)

SCAN: Smallest CANonical Correlation Method for ARIMA(p,d,q)

Estimates of Squared Canonical Correlation

MA-0   MA-1   MA-2   MA-3   MA-4   MA-5   MA-6
AR-0 0.0614 0.0104 0.1862 0.3516 0.0971 0.0128 0.0000
AR-1 0.0302 0.0294 0.1501 0.0943 0.0855 0.0127 0.0385
AR-2 0.3070 0.2781 0.2140 0.0006 0.1589 0.1884 0.2243
AR-3 0.1627 0.0037 0.1927 0.2311 0.1379 0.0207 0.0376
AR-4 0.2087 0.3947 0.3653 0.3075 0.1502 0.1364 0.1013
AR-5 0.1677 0.1219 0.0110 0.0263 0.0332 0.0350 0.0044
AR-6 0.0114 0.0485 0.0561 0.0427 0.0009 0.0089 0.0308

C(m,j)

MA-0    MA-1    MA-2    MA-3   MA-4   MA-5    MA-6
AR-0  4.1161  0.6585 12.0315 20.6512 4.5388 0.5620  0.0000
AR-1  6.1127  1.9499  9.9356  4.9145 4.7219 0.4642  1.9015
AR-2 72.6011 19.1679 14.3512  0.0337 7.9668 9.6479 11.4573
AR-3 34.9724  0.2386 10.1620 13.4082 6.7875 0.8725  1.4071
AR-4 45.8691 27.5070 19.1422 20.2835 7.3339 5.5374  3.5874
AR-5 35.7981  8.0498  0.6280  1.3543 1.8470 1.7930  0.2338
AR-6  2.2147  3.1466  3.5990  1.9904 0.0511 0.4816  1.6440

Chi-Square(1) Test p-value

MA-0   MA-1   MA-2   MA-3   MA-4   MA-5   MA-6
AR-0 0.0425 0.4171 0.0005 0.0000 0.0331 0.4534 0.0000
AR-1 0.0134 0.1626 0.0016 0.0266 0.0298 0.4957 0.1679
AR-2 0.0000 0.0000 0.0002 0.8543 0.0048 0.0019 0.0007
AR-3 0.0000 0.6252 0.0014 0.0003 0.0092 0.3503 0.2355
AR-4 0.0000 0.0000 0.0000 0.0000 0.0068 0.0186 0.0582
AR-5 0.0000 0.0046 0.4281 0.2445 0.1741 0.1806 0.6287
AR-6 0.1367 0.0761 0.0578 0.1583 0.8212 0.4877 0.1998

SCAN Matrix

MA-0 MA-1 MA-2 MA-3 MA-4 MA-5 MA-6
AR-0 "O"  "O"  "X"  "X"  "O"  "O"  "X"
AR-1 "O"  "O"  "X"  "O"  "O"  "O"  "O"
AR-2 "X"  "X"  "X"  "O"  "X"  "X"  "X"
AR-3 "X"  "O"  "X"  "X"  "X"  "O"  "O"
AR-4 "X"  "X"  "X"  "X"  "X"  "O"  "O"
AR-5 "X"  "X"  "O"  "O"  "O"  "O"  "O"
AR-6 "O"  "O"  "O"  "O"  "O"  "O"  "O"```

with the following graph

Of course, it is possible to ask for larger values,

`> arma.scan(s,12,12)`

The graph is now

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",
> 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
[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 () 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 []