[This article was first published on Freakonometrics » R-english, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

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 language="(1,2,3,5)"][/c]
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

To leave a comment for the author, please follow the link and comment on their blog: Freakonometrics » R-english.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

# Never miss an update! Subscribe to R-bloggers to receive e-mails with the latest R posts.(You will not see this message again.)

Click here to close (This popup will not appear again)