x <- c(235.1
,280.7
,264.6
,240.7
,201.4
,240.8
,241.1
,223.8
,206.1
,174.7
,203.3
,220.5
,299.5
,347.4
,338.3
,327.7
,351.6
,396.6
,438.8
,395.6
,363.5
,378.8
,357
,369
,464.8
,479.1
,431.3
,366.5
,326.3
,355.1
,331.6
,261.3
,249
,205.5
,235.6
,240.9
,264.9
,253.8
,232.3
,193.8
,177
,213.2
,207.2
,180.6
,188.6
,175.4
,199
,179.6
,225.8
,234
,200.2
,183.6
,178.2
,203.2
,208.5
,191.8
,172.8
,148
,159.4
,154.5
,213.2
,196.4
,182.8
,176.4
,153.6
,173.2
,171
,151.2
,161.9
,157.2
,201.7
,236.4
,356.1
,398.3
,403.7
,384.6
,365.8
,368.1
,367.9
,347
,343.3
,292.9
,311.5
,300.9
,366.9
,356.9
,329.7
,316.2
,269
,289.3
,266.2
,253.6
,233.8
,228.4
,253.6
,260.1
,306.6
,309.2
,309.5
,271
,279.9
,317.9
,298.4
,246.7
,227.3
,209.1
,259.9
,266
,320.6
,308.5
,282.2
,262.7
,263.5
,313.1
,284.3
,252.6
,250.3
,246.5
,312.7
,333.2
,446.4
,511.6
,515.5
,506.4
,483.2
,522.3
,509.8
,460.7
,405.8
,375
,378.5
,406.8
,467.8
,469.8
,429.8
,355.8
,332.7
,378
,360.5
,334.7
,319.5
,323.1
,363.6
,352.1
,411.9
,388.6
,416.4
,360.7
,338
,417.2
,388.4
,371.1
,331.5
,353.7
,396.7
,447
,533.5
,565.4
,542.3
,488.7
,467.1
,531.3
,496.1
,444
,403.4
,386.3
,394.1
,404.1
,462.1
,448.1
,432.3
,386.3
,395.2
,421.9
,382.9
,384.2
,345.5
,323.4
,372.6
,376
,462.7
,487
,444.2
,399.3
,394.9
,455.4
,414
,375.5
,347
,339.4
,385.8
,378.8
,451.8
,446.1
,422.5
,383.1
,352.8
,445.3
,367.5
,355.1
,326.2
,319.8
,331.8
,340.9
,394.1
,417.2
,369.9
,349.2
,321.4
,405.7
,342.9
,316.5
,284.2
,270.9
,288.8
,278.8
,324.4
,310.9
,299
,273
,279.3
,359.2
,305
,282.1
,250.3
,246.5
,257.9
,266.5
,315.9
,318.4
,295.4
,266.4
,245.8
,362.8
,324.9
,294.2
,289.5
,295.2
,290.3
,272
,307.4
,328.7
,292.9
,249.1
,230.4
,361.5
,321.7
,277.2
,260.7
,251
,257.6
,241.8
,287.5
,292.3
,274.7
,254.2
,230
,339
,318.2
,287
,295.8
,284
,271
,262.7
,340.6
,379.4
,373.3
,355.2
,338.4
,466.9
,451
,422
,429.2
,425.9
,460.7
,463.6
,541.4
,544.2
,517.5
,469.4
,439.4
,549
,533
,506.1
,484
,457
,481.5
,469.5
,544.7
,541.2
,521.5
,469.7
,434.4
,542.6
,517.3
,485.7
,465.8
,447
,426.6
,411.6
,467.5
,484.5
,451.2
,417.4
,379.9
,484.7
,455
,420.8
,416.5
,376.3
,405.6
,405.8
,500.8
,514
,475.5
,430.1
,414.4
,538
,526
,488.5
,520.2
,504.4
,568.5
,610.6
,818
,830.9
,835.9
,782
,762.3
,856.9
,820.9
,769.6
,752.2
,724.4
,723.1
,719.5
,817.4
,803.3
,752.5
,689
,630.4
,765.5
,757.7
,732.2
,702.6
,683.3
,709.5
,702.2
,784.8
,810.9
,755.6
,656.8
,615.1
,745.3
,694.1
,675.7
,643.7
,622.1
,634.6
,588
,689.7
,673.9
,647.9
,568.8
,545.7
,632.6
,643.8
,593.1
,579.7
,546
,562.9
,572.5)
par9 = '1'
par8 = '2'
par7 = '1'
par6 = '3'
par5 = '12'
par4 = '1'
par3 = '1'
par2 = '0.5'
par1 = 'FALSE'
#'GNU S' R Code compiled by R2WASP v. 1.0.44 ()
#Author: Prof. Dr. P. Wessa
#To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/
#Source of accompanying publication: Office for Research, Development, and Education
#Technical description: Write here your technical program description (don't use hard returns!)
library(lattice)
if (par1 == 'TRUE') par1 <- TRUE
if (par1 == 'FALSE') par1 <- FALSE
par2 <- as.numeric(par2) #Box-Cox lambda transformation parameter
par3 <- as.numeric(par3) #degree of non-seasonal differencing
par4 <- as.numeric(par4) #degree of seasonal differencing
par5 <- as.numeric(par5) #seasonal period
par6 <- as.numeric(par6) #degree (p) of the non-seasonal AR(p) polynomial
par7 <- as.numeric(par7) #degree (q) of the non-seasonal MA(q) polynomial
par8 <- as.numeric(par8) #degree (P) of the seasonal AR(P) polynomial
par9 <- as.numeric(par9) #degree (Q) of the seasonal MA(Q) polynomial
armaGR <- function(arima.out, names, n){
try1 <- arima.out$coef
try2 <- sqrt(diag(arima.out$var.coef))
try.data.frame  <- data.frame(matrix(NA,ncol=4,nrow=length(names)))
dimnames(try.data.frame) <- list(names,c('coef','std','tstat','pv'))
try.data.frame[,1] <- try1
for(i in 1:length(try2)) try.data.frame[which(rownames(try.data.frame)==names(try2)[i]),2] <- try2[i]
try.data.frame[,3] <- try.data.frame[,1] / try.data.frame[,2]
try.data.frame[,4] <- round((1-pt(abs(try.data.frame[,3]),df=n-(length(try2)+1)))*2,5)
vector <- rep(NA,length(names))
vector[is.na(try.data.frame[,4])] <- 0
maxi <- which.max(try.data.frame[,4])
continue <- max(try.data.frame[,4],na.rm=TRUE) > .05
vector[maxi] <- 0
list(summary=try.data.frame,next.vector=vector,continue=continue)
}
arimaSelect <- function(series, order=c(13,0,0), seasonal=list(order=c(2,0,0),period=12), include.mean=F){
nrc <- order[1]+order[3]+seasonal$order[1]+seasonal$order[3]
coeff <- matrix(NA, nrow=nrc*2, ncol=nrc)
pval <- matrix(NA, nrow=nrc*2, ncol=nrc)
mylist <- rep(list(NULL), nrc)
names  <- NULL
if(order[1] > 0) names <- paste('ar',1:order[1],sep='')
if(order[3] > 0) names <- c( names , paste('ma',1:order[3],sep='') )
if(seasonal$order[1] > 0) names <- c(names, paste('sar',1:seasonal$order[1],sep=''))
if(seasonal$order[3] > 0) names <- c(names, paste('sma',1:seasonal$order[3],sep=''))
arima.out <- arima(series, order=order, seasonal=seasonal, include.mean=include.mean, method='ML')
mylist[[1]] <- arima.out
last.arma <- armaGR(arima.out, names, length(series))
mystop <- FALSE
i <- 1
coeff[i,] <- last.arma[[1]][,1]
pval [i,] <- last.arma[[1]][,4]
i <- 2
aic <- arima.out$aic
while(!mystop){
mylist[[i]] <- arima.out
arima.out <- arima(series, order=order, seasonal=seasonal, include.mean=include.mean, method='ML', fixed=last.arma$next.vector)
aic <- c(aic, arima.out$aic)
last.arma <- armaGR(arima.out, names, length(series))
mystop <- !last.arma$continue
coeff[i,] <- last.arma[[1]][,1]
pval [i,] <- last.arma[[1]][,4]
i <- i+1
}
list(coeff, pval, mylist, aic=aic)
}
arimaSelectplot <- function(arimaSelect.out,noms,choix){
noms <- names(arimaSelect.out[[3]][[1]]$coef)
coeff <- arimaSelect.out[[1]]
k <- min(which(is.na(coeff[,1])))-1
coeff <- coeff[1:k,]
pval  <- arimaSelect.out[[2]][1:k,]
aic   <- arimaSelect.out$aic[1:k]
coeff[coeff==0] <- NA
n <- ncol(coeff)
if(missing(choix)) choix <- k
layout(matrix(c(1,1,1,2,
3,3,3,2,
3,3,3,4,
5,6,7,7),nr=4),
widths=c(10,35,45,15),
heights=c(30,30,15,15))
couleurs <- rainbow(75)[1:50]#(50)
ticks <- pretty(coeff)
par(mar=c(1,1,3,1))
plot(aic,k:1-.5,type='o',pch=21,bg='blue',cex=2,axes=F,lty=2,xpd=NA)
points(aic[choix],k-choix+.5,pch=21,cex=4,bg=2,xpd=NA)
title('aic',line=2)
par(mar=c(3,0,0,0))
plot(0,axes=F,xlab='',ylab='',xlim=range(ticks),ylim=c(.1,1))
rect(xleft  = min(ticks) + (0:49)/50*(max(ticks)-min(ticks)),
xright = min(ticks) + (1:50)/50*(max(ticks)-min(ticks)),
ytop   = rep(1,50),
ybottom= rep(0,50),col=couleurs,border=NA)
axis(1,ticks)
rect(xleft=min(ticks),xright=max(ticks),ytop=1,ybottom=0)
text(mean(coeff,na.rm=T),.5,'coefficients',cex=2,font=2)
par(mar=c(1,1,3,1))
image(1:n,1:k,t(coeff[k:1,]),axes=F,col=couleurs,zlim=range(ticks))
for(i in 1:n) for(j in 1:k) if(!is.na(coeff[j,i])) {
if(pval[j,i]<.01)                            symb = 'green'
else if( (pval[j,i]<.05) & (pval[j,i]>=.01)) symb = 'orange'
else if( (pval[j,i]<.1)  & (pval[j,i]>=.05)) symb = 'red'
else                                         symb = 'black'
polygon(c(i+.5   ,i+.2   ,i+.5   ,i+.5),
c(k-j+0.5,k-j+0.5,k-j+0.8,k-j+0.5),
col=symb)
if(j==choix)  {
rect(xleft=i-.5,
xright=i+.5,
ybottom=k-j+1.5,
ytop=k-j+.5,
lwd=4)
text(i,
k-j+1,
round(coeff[j,i],2),
cex=1.2,
font=2)
}
else{
rect(xleft=i-.5,xright=i+.5,ybottom=k-j+1.5,ytop=k-j+.5)
text(i,k-j+1,round(coeff[j,i],2),cex=1.2,font=1)
}
}
axis(3,1:n,noms)
par(mar=c(0.5,0,0,0.5))
plot(0,axes=F,xlab='',ylab='',type='n',xlim=c(0,8),ylim=c(-.2,.8))
cols <- c('green','orange','red','black')
niv  <- c('0','0.01','0.05','0.1')
for(i in 0:3){
polygon(c(1+2*i   ,1+2*i   ,1+2*i-.5   ,1+2*i),
c(.4      ,.7      , .4        , .4),
col=cols[i+1])
text(2*i,0.5,niv[i+1],cex=1.5)
}
text(8,.5,1,cex=1.5)
text(4,0,'p-value',cex=2)
box()
residus <- arimaSelect.out[[3]][[choix]]$res
par(mar=c(1,2,4,1))
acf(residus,main='')
title('acf',line=.5)
par(mar=c(1,2,4,1))
pacf(residus,main='')
title('pacf',line=.5)
par(mar=c(2,2,4,1))
qqnorm(residus,main='')
title('qq-norm',line=.5)
qqline(residus)
residus
}
if (par2 == 0) x <- log(x)
if (par2 != 0) x <- x^par2
(selection <- arimaSelect(x, order=c(par6,par3,par7), seasonal=list(order=c(par8,par4,par9), period=par5)))
postscript(file="/var/www/html/rcomp/tmp/141851228837660.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) 
resid <- arimaSelectplot(selection)
dev.off()
resid
postscript(file="/var/www/html/rcomp/tmp/2otzc1228837660.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) 
acf(resid,length(resid)/2, main='Residual Autocorrelation Function')
dev.off()
postscript(file="/var/www/html/rcomp/tmp/3ygb11228837660.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) 
pacf(resid,length(resid)/2, main='Residual Partial Autocorrelation Function')
dev.off()
postscript(file="/var/www/html/rcomp/tmp/4crc11228837660.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) 
cpgram(resid, main='Residual Cumulative Periodogram')
dev.off()
postscript(file="/var/www/html/rcomp/tmp/50izw1228837660.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) 
hist(resid, main='Residual Histogram', xlab='values of Residuals')
dev.off()
postscript(file="/var/www/html/rcomp/tmp/6x0331228837660.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) 
densityplot(~resid,col='black',main='Residual Density Plot', xlab='values of Residuals')
dev.off()
postscript(file="/var/www/html/rcomp/tmp/7mj571228837661.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) 
qqnorm(resid, main='Residual Normal Q-Q Plot')
qqline(resid)
dev.off()
ncols <- length(selection[[1]][1,])
nrows <- length(selection[[2]][,1])-1

#Note: the /var/www/html/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
load(file="/var/www/html/rcomp/createtable")

a<-table.start()
a<-table.row.start(a)
a<-table.element(a,'ARIMA Parameter Estimation and Backward Selection', ncols+1,TRUE)
a<-table.row.end(a)
a<-table.row.start(a)
a<-table.element(a,'Iteration', header=TRUE)
for (i in 1:ncols) {
a<-table.element(a,names(selection[[3]][[1]]$coef)[i],header=TRUE)
}
a<-table.row.end(a)
for (j in 1:nrows) {
a<-table.row.start(a)
mydum <- 'Estimates ('
mydum <- paste(mydum,j)
mydum <- paste(mydum,')')
a<-table.element(a,mydum, header=TRUE)
for (i in 1:ncols) {
a<-table.element(a,round(selection[[1]][j,i],4))
}
a<-table.row.end(a)
a<-table.row.start(a)
a<-table.element(a,'(p-val)', header=TRUE)
for (i in 1:ncols) {
mydum <- '('
mydum <- paste(mydum,round(selection[[2]][j,i],4),sep='')
mydum <- paste(mydum,')')
a<-table.element(a,mydum)
}
a<-table.row.end(a)
}
a<-table.end(a)
table.save(a,file="/var/www/html/rcomp/tmp/8erb81228837661.tab") 
a<-table.start()
a<-table.row.start(a)
a<-table.element(a,'Estimated ARIMA Residuals', 1,TRUE)
a<-table.row.end(a)
a<-table.row.start(a)
a<-table.element(a,'Value', 1,TRUE)
a<-table.row.end(a)
for (i in (par4*par5+par3):length(resid)) {
a<-table.row.start(a)
a<-table.element(a,resid[i])
a<-table.row.end(a)
}
a<-table.end(a)
table.save(a,file="/var/www/html/rcomp/tmp/98hod1228837661.tab") 

system("convert tmp/141851228837660.ps tmp/141851228837660.png")
system("convert tmp/2otzc1228837660.ps tmp/2otzc1228837660.png")
system("convert tmp/3ygb11228837660.ps tmp/3ygb11228837660.png")
system("convert tmp/4crc11228837660.ps tmp/4crc11228837660.png")
system("convert tmp/50izw1228837660.ps tmp/50izw1228837660.png")
system("convert tmp/6x0331228837660.ps tmp/6x0331228837660.png")
system("convert tmp/7mj571228837661.ps tmp/7mj571228837661.png")