Home » date » 2010 » Dec » 14 »

*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 18:42:32 +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/t1292352171v2retoqa1ltgmm0.htm/, Retrieved Mon, 20 May 2013 06:29:42 +0000   Original text written by user:   IsPrivate? No (this computation is public)   User-defined keywords:   System-generated keywords (parent): t122876443062pvsyb6eh8ar26 (pk = 30802) Estimated Impact 41   Dataseries X: » Textfile « » CSV « » Stem and Leaf « » Histogram « » Kernel Density « » Harrell-Davis Quantiles « » Central Tendency « » Variability « 655362 873127 1107897 1555964 1671159 1493308 2957796 2638691 1305669 1280496 921900 867888 652586 913831 1108544 1555827 1699283 1509458 3268975 2425016 1312703 1365498 934453 775019 651142 843192 1146766 1652601 1465906 1652734 2922334 2702805 1458956 1410363 1019279 936574 708917 885295 1099663 1576220 1487870 1488635 2882530 2677026 1404398 1344370 936865 872705 628151 953712 1160384 1400618 1661511 1495347 2918786 2775677 1407026 1370199 964526 850851 683118 847224 1073256 1514326 1503734 1507712 2865698 2788128 1391596 1366378 946295 859626   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 19 seconds R Server 'George Udny Yule' @ 72.249.76.132 ARIMA Parameter Estimation and Backward Selection Iteration ar1 ar2 ar3 ma1 sar1 sar2 sma1 Estimates ( 1 ) -1.0256 -0.6136 -0.1521 -0.9983 -0.7658 -0.4072 0.2762 (p-val) (0 ) (8e-04 ) (0.2388 ) (0 ) (0.1818 ) (0.0904 ) (0.6575 ) Estimates ( 2 ) -1.0179 -0.5968 -0.1499 -0.9982 -0.5142 -0.3081 0 (p-val) (0 ) (8e-04 ) (0.2452 ) (0 ) (2e-04 ) (0.0732 ) (NA ) Estimates ( 3 ) -0.9402 -0.4461 0 -1.0016 -0.5369 -0.3152 0 (p-val) (0 ) (2e-04 ) (NA ) (0 ) (1e-04 ) (0.0658 ) (NA ) Estimates ( 4 ) -0.9742 -0.4468 0 -1.0014 -0.4618 0 0 (p-val) (0 ) (2e-04 ) (NA ) (0 ) (3e-04 ) (NA ) (NA ) Estimates ( 5 ) NA NA NA NA NA NA NA (p-val) (NA ) (NA ) (NA ) (NA ) (NA ) (NA ) (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 5541.25053290481 -26717.2631373883 -11096.8938495013 14084.4649923033 13010.9061967 223945.398198488 -251634.409606684 -119098.276978991 50755.3492216522 85286.7486171569 -118295.357417467 -27949.5984842587 -8831.6318054938 81876.7986759022 111213.435799922 -209439.267485816 102709.407732824 -142057.309692732 224215.15055235 159976.984307709 50329.5963208744 -88286.3930182652 -2243.52417341197 -45436.8202021278 -95395.997962044 -110290.858203990 -51662.9063146404 -76675.0534672818 -36681.2975568886 -46055.5753885475 161247.512431262 127419.962009706 12952.7536599314 -67934.4670124602 -4593.51767474543 -22020.9682773320 92522.2675639386 69956.299209602 -198425.794955626 73598.859296121 30371.3686477884 -61250.7054894422 147406.287862065 69233.3971561549 -42776.0057326601 -77570.6822775456 -16129.2369725288 23547.726738058 -59135.8641103752 -75223.804557123 18080.0748124718 5755.79845291289 4421.23825175518 -9909.98704405827 105634.782959724 4954.18782811865 -15609.1018257305 -37191.883874967 -11.6995275141632   Charts produced by software: http://www.freestatistics.org/blog/date/2010/Dec/14/t1292352171v2retoqa1ltgmm0/1fyyx1292352132.png (opens in new window) http://www.freestatistics.org/blog/date/2010/Dec/14/t1292352171v2retoqa1ltgmm0/1fyyx1292352132.ps (opens in new window) Click here to open pdf file. http://www.freestatistics.org/blog/date/2010/Dec/14/t1292352171v2retoqa1ltgmm0/2fyyx1292352132.png (opens in new window) http://www.freestatistics.org/blog/date/2010/Dec/14/t1292352171v2retoqa1ltgmm0/2fyyx1292352132.ps (opens in new window) Click here to open pdf file. http://www.freestatistics.org/blog/date/2010/Dec/14/t1292352171v2retoqa1ltgmm0/3q7g01292352132.png (opens in new window) http://www.freestatistics.org/blog/date/2010/Dec/14/t1292352171v2retoqa1ltgmm0/3q7g01292352132.ps (opens in new window) Click here to open pdf file. http://www.freestatistics.org/blog/date/2010/Dec/14/t1292352171v2retoqa1ltgmm0/4q7g01292352132.png (opens in new window) http://www.freestatistics.org/blog/date/2010/Dec/14/t1292352171v2retoqa1ltgmm0/4q7g01292352132.ps (opens in new window) Click here to open pdf file. http://www.freestatistics.org/blog/date/2010/Dec/14/t1292352171v2retoqa1ltgmm0/5q7g01292352132.png (opens in new window) http://www.freestatistics.org/blog/date/2010/Dec/14/t1292352171v2retoqa1ltgmm0/5q7g01292352132.ps (opens in new window) Click here to open pdf file. http://www.freestatistics.org/blog/date/2010/Dec/14/t1292352171v2retoqa1ltgmm0/6q7g01292352132.png (opens in new window) http://www.freestatistics.org/blog/date/2010/Dec/14/t1292352171v2retoqa1ltgmm0/6q7g01292352132.ps (opens in new window) Click here to open pdf file. http://www.freestatistics.org/blog/date/2010/Dec/14/t1292352171v2retoqa1ltgmm0/7jhf31292352132.png (opens in new window) http://www.freestatistics.org/blog/date/2010/Dec/14/t1292352171v2retoqa1ltgmm0/7jhf31292352132.ps (opens in new window) Click here to open pdf file.   Parameters (Session): par1 = FALSE ; par2 = 1 ; par3 = 2 ; par4 = 1 ; par5 = 12 ; par6 = 3 ; par7 = 1 ; par8 = 2 ; par9 = 1 ;   Parameters (R input): par1 = FALSE ; par2 = 1 ; par3 = 2 ; 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')