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

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


