Home » date » 2010 » Dec » 14 »

ws 9 7

*The author of this computation has been verified*

R Software Module: /rwasp_arimabackwardselection.wasp (opens new window with default values)

Title produced by software: ARIMA Backward Selection

Date of computation: Tue, 14 Dec 2010 19:14:46 +0000

Cite this page as follows:

Statistical Computations at FreeStatistics.org, Office for Research Development and Education, URL http://www.freestatistics.org/blog/date/2010/Dec/14/t1292353978u36ub0ze50d9rzj.htm/, Retrieved Fri, 24 May 2013 20:04:51 +0000

Original text written by user:

IsPrivate?

No (this computation is public)

User-defined keywords:

System-generated keywords (parent):

t1291051354cz4s5yz0dopvtx8 (pk = 102963)

Estimated Impact

Dataseries X:

» Textfile « » CSV « » Stem and Leaf « » Histogram « » Kernel Density « » Harrell-Davis Quantiles « » Central Tendency « » Variability «

Output produced by software:

Summary of computational transaction
Raw Input	view raw input (R code)
Raw Output	view raw output of R engine
Computing time	10 seconds
R Server	'Sir Ronald Aylmer Fisher' @ 193.190.124.24

ARIMA Parameter Estimation and Backward Selection
Iteration	ar1	ar2	ar3	ma1	sar1	sar2	sma1
Estimates ( 1 )	-0.3995	0.0639	0.1858	0.1506	-0.7323	-0.3971	0.183
(p-val)	(0.4771 )	(0.7441 )	(0.1433 )	(0.7911 )	(0.2421 )	(0.1575 )	(0.7882 )
Estimates ( 2 )	-0.2546	0.0995	0.1779	0	-0.7335	-0.3965	0.1899
(p-val)	(0.0521 )	(0.4567 )	(0.1662 )	(NA )	(0.2436 )	(0.1556 )	(0.7814 )
Estimates ( 3 )	-0.2494	0.1067	0.1698	0	-0.5607	-0.3248	0
(p-val)	(0.0543 )	(0.4148 )	(0.1758 )	(NA )	(0 )	(0.0447 )	(NA )
Estimates ( 4 )	-0.271	0	0.1451	0	-0.5702	-0.3405	0
(p-val)	(0.0345 )	(NA )	(0.2359 )	(NA )	(0 )	(0.0326 )	(NA )
Estimates ( 5 )	-0.2532	0	0	0	-0.574	-0.3376	0
(p-val)	(0.0476 )	(NA )	(NA )	(NA )	(0 )	(0.033 )	(NA )
Estimates ( 6 )	NA	NA	NA	NA	NA	NA	NA
(p-val)	(NA )	(NA )	(NA )	(NA )	(NA )	(NA )	(NA )
Estimates ( 7 )	NA	NA	NA	NA	NA	NA	NA
(p-val)	(NA )	(NA )	(NA )	(NA )	(NA )	(NA )	(NA )
Estimates ( 8 )	NA	NA	NA	NA	NA	NA	NA
(p-val)	(NA )	(NA )	(NA )	(NA )	(NA )	(NA )	(NA )
Estimates ( 9 )	NA	NA	NA	NA	NA	NA	NA
(p-val)	(NA )	(NA )	(NA )	(NA )	(NA )	(NA )	(NA )
Estimates ( 10 )	NA	NA	NA	NA	NA	NA	NA
(p-val)	(NA )	(NA )	(NA )	(NA )	(NA )	(NA )	(NA )
Estimates ( 11 )	NA	NA	NA	NA	NA	NA	NA
(p-val)	(NA )	(NA )	(NA )	(NA )	(NA )	(NA )	(NA )
Estimates ( 12 )	NA	NA	NA	NA	NA	NA	NA
(p-val)	(NA )	(NA )	(NA )	(NA )	(NA )	(NA )	(NA )
Estimates ( 13 )	NA	NA	NA	NA	NA	NA	NA
(p-val)	(NA )	(NA )	(NA )	(NA )	(NA )	(NA )	(NA )

Estimated ARIMA Residuals
Value
867.887509505211
-2250.28069676076
33618.3570398550
9954.34467606629
354.19173153525
18882.4063857050
20229.4310873144
268402.416151264
-113346.926118834
-45016.3941980673
35069.8612498888
58531.0958011903
-77256.3771257806
-31473.5945834065
-52391.0075182321
32854.9847993680
101107.732844000
-176275.960400893
79531.8844457594
-176414.251618811
151290.579592969
167731.594072401
143237.122503809
80251.9664432193
118735.726196695
75035.8258514422
19198.3084889299
-36364.5639798847
-36170.5787601385
-109567.395076199
-100783.336067031
-149267.403955552
38947.3511523581
58613.0601346168
16074.4602556389
-41563.0049451506
-15970.596469234
-47563.9548284742
59595.3577509129
65897.8405210807
-166489.283208111
46312.3270186926
-15952.872342146
-87780.652231175
134744.172684142
75232.812196889
24408.7558701518
-15406.1404274264
-3766.75352115813
27197.2239798236
-46777.2890164611
-82472.8212389017
-35154.718472515
-46946.8699812614
-43641.5364342392
-54920.7084749253
54905.4038614514
-10509.5840596962
-13706.8046613529
-42347.6087745922
-28990.4687416657

Charts produced by software:

http://www.freestatistics.org/blog/date/2010/Dec/14/t1292353978u36ub0ze50d9rzj/1dwg81292354074.png (opens in new window)

http://www.freestatistics.org/blog/date/2010/Dec/14/t1292353978u36ub0ze50d9rzj/1dwg81292354074.ps (opens in new window)

Click here to open pdf file.

http://www.freestatistics.org/blog/date/2010/Dec/14/t1292353978u36ub0ze50d9rzj/255ft1292354074.png (opens in new window)

http://www.freestatistics.org/blog/date/2010/Dec/14/t1292353978u36ub0ze50d9rzj/255ft1292354074.ps (opens in new window)

Click here to open pdf file.

http://www.freestatistics.org/blog/date/2010/Dec/14/t1292353978u36ub0ze50d9rzj/355ft1292354074.png (opens in new window)

http://www.freestatistics.org/blog/date/2010/Dec/14/t1292353978u36ub0ze50d9rzj/355ft1292354074.ps (opens in new window)

Click here to open pdf file.

http://www.freestatistics.org/blog/date/2010/Dec/14/t1292353978u36ub0ze50d9rzj/4gxww1292354074.png (opens in new window)

http://www.freestatistics.org/blog/date/2010/Dec/14/t1292353978u36ub0ze50d9rzj/4gxww1292354074.ps (opens in new window)

Click here to open pdf file.

http://www.freestatistics.org/blog/date/2010/Dec/14/t1292353978u36ub0ze50d9rzj/5gxww1292354074.png (opens in new window)

http://www.freestatistics.org/blog/date/2010/Dec/14/t1292353978u36ub0ze50d9rzj/5gxww1292354074.ps (opens in new window)

Click here to open pdf file.

http://www.freestatistics.org/blog/date/2010/Dec/14/t1292353978u36ub0ze50d9rzj/6gxww1292354074.png (opens in new window)

http://www.freestatistics.org/blog/date/2010/Dec/14/t1292353978u36ub0ze50d9rzj/6gxww1292354074.ps (opens in new window)

Click here to open pdf file.

http://www.freestatistics.org/blog/date/2010/Dec/14/t1292353978u36ub0ze50d9rzj/7gxww1292354074.png (opens in new window)

http://www.freestatistics.org/blog/date/2010/Dec/14/t1292353978u36ub0ze50d9rzj/7gxww1292354074.ps (opens in new window)

Click here to open pdf file.

Parameters (Session):

par1 = FALSE ; par2 = 1 ; par3 = 0 ; par4 = 1 ; par5 = 12 ; par6 = 3 ; par7 = 1 ; par8 = 2 ; par9 = 1 ;

Parameters (R input):

par1 = FALSE ; par2 = 1 ; par3 = 0 ; par4 = 1 ; par5 = 12 ; par6 = 3 ; par7 = 1 ; par8 = 2 ; par9 = 1 ;

R code (references can be found in the software module):

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)))
bitmap(file='test1.png')
resid <- arimaSelectplot(selection)
dev.off()
resid
bitmap(file='test2.png')
acf(resid,length(resid)/2, main='Residual Autocorrelation Function')
dev.off()
bitmap(file='test3.png')
pacf(resid,length(resid)/2, main='Residual Partial Autocorrelation Function')
dev.off()
bitmap(file='test4.png')
cpgram(resid, main='Residual Cumulative Periodogram')
dev.off()
bitmap(file='test5.png')
hist(resid, main='Residual Histogram', xlab='values of Residuals')
dev.off()
bitmap(file='test6.png')
densityplot(~resid,col='black',main='Residual Density Plot', xlab='values of Residuals')
dev.off()
bitmap(file='test7.png')
qqnorm(resid, main='Residual Normal Q-Q Plot')
qqline(resid)
dev.off()
ncols <- length(selection[[1]][1,])
nrows <- length(selection[[2]][,1])-1
load(file='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='mytable.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='mytable1.tab')