Home » date » 2009 » Dec » 04 »

*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: Fri, 04 Dec 2009 07:52:27 -0700   Cite this page as follows: Statistical Computations at FreeStatistics.org, Office for Research Development and Education, URL http://www.freestatistics.org/blog/date/2009/Dec/04/t1259938451456sbedp4x0hfsa.htm/, Retrieved Mon, 20 May 2013 06:15:47 +0000   Original text written by user:   IsPrivate? No (this computation is public)   User-defined keywords:   System-generated keywords (parent): t1259333643wfudtum7qa0jfpl (pk = 60855) Estimated Impact 41   Dataseries X: » Textfile « » CSV « » Stem and Leaf « » Histogram « » Kernel Density « » Harrell-Davis Quantiles « » Central Tendency « » Variability « 455 456 517 525 523 519 509 512 519 517 510 509 501 507 569 580 578 565 547 555 562 561 555 544 537 543 594 611 613 611 594 595 591 589 584 573 567 569 621 629 628 612 595 597 593 590 580 574 573 573 620 626 620 588 566 557 561 549 532 526 511 499 555 565 542 527 510 514 517 508 493 490 469 478 528 534 518 506 502 516 528 533 536 537 524 536 587 597 581 564   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 11 seconds R Server 'Gwilym Jenkins' @ 72.249.127.135 ARIMA Parameter Estimation and Backward Selection Iteration ar1 ar2 ar3 ma1 sar1 sar2 sma1 Estimates ( 1 ) 0.4935 0.0811 0.1339 -0.3806 0.2603 -0.1671 -0.7128 (p-val) (0.4176 ) (0.6039 ) (0.4744 ) (0.538 ) (0.3849 ) (0.346 ) (0.0565 ) Estimates ( 2 ) 0.6843 0 0.1166 -0.549 0.2275 -0.1856 -0.6738 (p-val) (0.1095 ) (NA ) (0.5358 ) (0.2153 ) (0.4385 ) (0.2836 ) (0.0554 ) Estimates ( 3 ) 0.9097 0 0 -0.764 0.251 -0.1916 -1.4161 (p-val) (0 ) (NA ) (NA ) (0 ) (0.3601 ) (0.2537 ) (0.0392 ) Estimates ( 4 ) 0.9255 0 0 -0.7829 0 -0.2754 -0.4342 (p-val) (0 ) (NA ) (NA ) (0 ) (NA ) (0.0288 ) (0.0042 ) 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 -0.0745156793299108 0.0663929577506382 -0.0406261158725139 0.0290516062801986 -0.00677617991735565 -0.120303035474756 -0.0884683200501621 0.0872577839056256 0.00471520051978586 0.0212330558570349 0.0205634498903666 -0.131497714198925 0.0355228346858527 0.0256692441598289 -0.171271626992804 0.113290258701022 0.0639798075234621 0.120365117506391 -0.0184123966156435 -0.0914572992045392 -0.145975433798118 0.0130436514982875 0.0376814989807520 -0.0182317663091494 0.0325564497129637 -0.0282744745341756 -0.0610862395531608 -0.0700926502438601 0.00394887132382057 -0.157538574303815 0.0149228151799967 0.0409189674652656 -0.0248275529568095 0.0143712407963754 -0.0362779014450171 0.0487239154363810 0.097526855974283 -0.0357701892407696 -0.146993520377629 -0.0106756514755384 -0.0366975140112644 -0.239514192000376 -0.0312398056855772 -0.121629436121252 0.120921197262675 -0.0955017162080598 -0.0728496219231542 0.0541999120080228 -0.128535923222649 -0.154248267436466 0.171727796378330 0.0663859935067707 -0.258803180604707 0.131870180425831 0.0574389396703112 0.143087609195562 -0.00274337572772709 -0.0222467357479128 -0.0455175435107223 0.0671414078409635 -0.155371697749837 0.248978128480308 -0.082192697269453 -0.0650579057031216 -0.0437437865089546 -0.000840045079394987 0.193914753728192 0.138576586132039 0.122876147829280 0.111044113809799 0.155718998385065 0.00944240860297473 -0.0591282157649024 0.0173833631705559 -0.0769783747924267 -0.0044840780664122 -0.132728047946625 -0.0529043475195727   Charts produced by software: http://www.freestatistics.org/blog/date/2009/Dec/04/t1259938451456sbedp4x0hfsa/1yquq1259938335.png (opens in new window) http://www.freestatistics.org/blog/date/2009/Dec/04/t1259938451456sbedp4x0hfsa/1yquq1259938335.ps (opens in new window) Click here to open pdf file. http://www.freestatistics.org/blog/date/2009/Dec/04/t1259938451456sbedp4x0hfsa/27q9k1259938335.png (opens in new window) http://www.freestatistics.org/blog/date/2009/Dec/04/t1259938451456sbedp4x0hfsa/27q9k1259938335.ps (opens in new window) Click here to open pdf file. http://www.freestatistics.org/blog/date/2009/Dec/04/t1259938451456sbedp4x0hfsa/364hx1259938335.png (opens in new window) http://www.freestatistics.org/blog/date/2009/Dec/04/t1259938451456sbedp4x0hfsa/364hx1259938335.ps (opens in new window) Click here to open pdf file. http://www.freestatistics.org/blog/date/2009/Dec/04/t1259938451456sbedp4x0hfsa/4dde71259938335.png (opens in new window) http://www.freestatistics.org/blog/date/2009/Dec/04/t1259938451456sbedp4x0hfsa/4dde71259938335.ps (opens in new window) Click here to open pdf file. http://www.freestatistics.org/blog/date/2009/Dec/04/t1259938451456sbedp4x0hfsa/589i41259938335.png (opens in new window) http://www.freestatistics.org/blog/date/2009/Dec/04/t1259938451456sbedp4x0hfsa/589i41259938335.ps (opens in new window) Click here to open pdf file. http://www.freestatistics.org/blog/date/2009/Dec/04/t1259938451456sbedp4x0hfsa/6ow001259938335.png (opens in new window) http://www.freestatistics.org/blog/date/2009/Dec/04/t1259938451456sbedp4x0hfsa/6ow001259938335.ps (opens in new window) Click here to open pdf file. http://www.freestatistics.org/blog/date/2009/Dec/04/t1259938451456sbedp4x0hfsa/7x3i01259938335.png (opens in new window) http://www.freestatistics.org/blog/date/2009/Dec/04/t1259938451456sbedp4x0hfsa/7x3i01259938335.ps (opens in new window) Click here to open pdf file.   Parameters (Session): par1 = FALSE ; par2 = 0.5 ; par3 = 1 ; par4 = 1 ; par5 = 12 ; par6 = 3 ; par7 = 1 ; par8 = 2 ; par9 = 1 ;   Parameters (R input): par1 = FALSE ; par2 = 0.5 ; par3 = 1 ; 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')