R version 2.13.0 (2011-04-13) Copyright (C) 2011 The R Foundation for Statistical Computing ISBN 3-900051-07-0 Platform: i486-pc-linux-gnu (32-bit) R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > x <- c(396 + ,297 + ,559 + ,967 + ,270 + ,143 + ,1562 + ,109 + ,371 + ,656 + ,511 + ,655 + ,465 + ,525 + ,885 + ,497 + ,1436 + ,612 + ,865 + ,385 + ,567 + ,639 + ,963 + ,398 + ,410 + ,966 + ,801 + ,892 + ,513 + ,469 + ,683 + ,643 + ,535 + ,625 + ,264 + ,992 + ,238 + ,818 + ,937 + ,70 + ,507 + ,260 + ,503 + ,927 + ,1269 + ,537 + ,910 + ,532 + ,345 + ,918 + ,1635 + ,330 + ,557 + ,1178 + ,740 + ,452 + ,218 + ,764 + ,255 + ,454 + ,866 + ,574 + ,1276 + ,379 + ,825 + ,798 + ,663 + ,1069 + ,921 + ,858 + ,711 + ,503 + ,382 + ,464 + ,717 + ,690 + ,462 + ,657 + ,385 + ,577 + ,619 + ,479 + ,817 + ,752 + ,430 + ,451 + ,537 + ,519 + ,1000 + ,637 + ,465 + ,437 + ,711 + ,299 + ,248 + ,1162 + ,714 + ,905 + ,649 + ,512 + ,472 + ,905 + ,786 + ,489 + ,479 + ,617 + ,925 + ,351 + ,1144 + ,669 + ,707 + ,458 + ,214 + ,599 + ,572 + ,897 + ,819 + ,720 + ,273 + ,508 + ,506 + ,451 + ,699 + ,407 + ,465 + ,245 + ,370 + ,316 + ,603 + ,154 + ,229 + ,577 + ,192 + ,617 + ,411 + ,975 + ,146 + ,705 + ,184 + ,200 + ,274 + ,502 + ,382 + ,964 + ,537 + ,438 + ,369 + ,417 + ,276 + ,514 + ,822 + ,389 + ,466 + ,1255 + ,694 + ,1024 + ,400 + ,397 + ,350 + ,719 + ,1277 + ,356 + ,457 + ,1402 + ,600 + ,480 + ,595 + ,436 + ,230 + ,651 + ,1367 + ,564 + ,716 + ,747 + ,467 + ,671 + ,861 + ,319 + ,612 + ,433 + ,434 + ,503 + ,85 + ,564 + ,824 + ,74 + ,259 + ,69 + ,535 + ,239 + ,438 + ,459 + ,426 + ,288 + ,498 + ,454 + ,376 + ,225 + ,555 + ,252 + ,208 + ,130 + ,481 + ,389 + ,565 + ,173 + ,278 + ,609 + ,422 + ,445 + ,387 + ,339 + ,181 + ,245 + ,384 + ,212 + ,399 + ,229 + ,224 + ,203 + ,333 + ,384 + ,636 + ,185 + ,93 + ,581 + ,248 + ,304 + ,344 + ,407 + ,170 + ,312 + ,507 + ,224 + ,340 + ,168 + ,443 + ,204 + ,367 + ,210 + ,335 + ,364 + ,178 + ,206 + ,279 + ,387 + ,490 + ,238 + ,343 + ,232 + ,530 + ,291 + ,67 + ,397 + ,467 + ,178 + ,175 + ,299 + ,154 + ,106 + ,189 + ,194 + ,135 + ,201 + ,207 + ,280 + ,260 + ,227 + ,239 + ,333 + ,428 + ,230 + ,292 + ,350 + ,186 + ,326 + ,155 + ,75 + ,361 + ,261 + ,299 + ,300 + ,450 + ,183 + ,238 + ,165 + ,234 + ,176 + ,329) > par9 = '1' > par8 = '0' > par7 = '1' > par6 = '0' > par5 = '1' > par4 = '0' > par3 = '0' > par2 = '-0.1' > par1 = 'FALSE' > 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 <- 7 #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))) [[1]] [,1] [,2] [1,] 0.7477302 0.8538630 [2,] 0.0000000 0.9555866 [3,] NA NA [4,] NA NA [[2]] [,1] [,2] [1,] 0 0 [2,] NA 0 [3,] NA NA [4,] NA NA [[3]] [[3]][[1]] Call: arima(x = series, order = order, seasonal = seasonal, include.mean = include.mean, method = "ML") Coefficients: ma1 sma1 0.7477 0.8539 s.e. 0.0263 0.0326 sigma^2 estimated as 0.03575: log likelihood = 66.17, aic = -126.35 [[3]][[2]] Call: arima(x = series, order = order, seasonal = seasonal, include.mean = include.mean, method = "ML") Coefficients: ma1 sma1 0.7477 0.8539 s.e. 0.0263 0.0326 sigma^2 estimated as 0.03575: log likelihood = 66.17, aic = -126.35 $aic [1] -126.3478 130.2748 > postscript(file="/var/wessaorg/rcomp/tmp/1iavh1324395415.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > resid <- arimaSelectplot(selection) > dev.off() null device 1 > resid Time Series: Start = 1 End = 289 Frequency = 1 [1] 0.334877243 0.209741944 0.249449058 0.199345888 0.285234174 [6] 0.250629359 0.099760718 0.230584398 0.073669579 0.171696758 [11] 0.124021314 0.118002683 0.122128127 0.196441759 0.111465358 [16] 0.246884714 0.106204204 0.233079813 0.150588711 0.251709036 [21] 0.104177680 0.228119346 0.060604168 0.259458519 0.092711967 [26] 0.165300160 0.090334420 0.189570006 0.136068223 0.240343788 [31] 0.084451627 0.218470986 0.168326577 0.214115932 0.191913881 [36] 0.123573377 0.196749159 0.139926934 0.161041409 0.250955924 [41] 0.062309320 0.229155619 0.138436807 0.156000717 0.129426972 [46] 0.209891828 0.036082366 0.293086791 0.102678863 0.165966897 [51] 0.132179293 0.250348727 0.083498293 0.264673009 0.046892696 [56] 0.233397866 0.201215213 0.145788448 0.168516425 0.185372716 [61] 0.091837856 0.251860358 0.071578631 0.179204075 0.124637265 [66] 0.182771511 0.120360900 0.210923544 0.074799537 0.230979542 [71] 0.146989444 0.206083087 0.162137241 0.200423335 0.111738483 [76] 0.237848096 0.118670323 0.161445948 0.161141903 0.139386961 [81] 0.147335845 0.205910349 0.083464700 0.200274646 0.181890753 [86] 0.166141883 0.187206224 0.180331921 0.096780079 0.249114517 [91] 0.130404581 0.163991453 0.138137750 0.196397314 0.155948145 [96] 0.179304379 0.110005962 0.153720952 0.185116828 0.174822388 [101] 0.153735250 0.132846147 0.161504951 0.209134015 0.181494664 [106] 0.134150730 0.137457807 0.210857711 0.125250432 0.205376110 [111] 0.083748389 0.190859919 0.211589243 0.166289331 0.137904719 [116] 0.162023498 0.134874878 0.214418291 0.193844903 0.088935583 [121] 0.192997260 0.174495657 0.162601241 0.208123192 0.116300903 [126] 0.187558149 0.213618934 0.181065657 0.119629536 0.264591617 [131] 0.101450614 0.221486548 0.191080510 0.081009960 0.196261755 [136] 0.137962687 0.202089350 0.112374913 0.255717642 0.092962626 [141] 0.309759793 0.086025709 0.244389374 0.059672174 0.263740845 [146] 0.057042057 0.268421855 0.022477108 0.282025716 0.061212231 [151] 0.258351213 0.094349909 0.253310018 0.034931624 0.303140304 [156] 0.018181373 0.303357261 0.063193052 0.263906690 0.044148372 [161] 0.264502226 0.076807133 0.275511431 0.007904171 0.273909072 [166] 0.068866359 0.270209565 0.088493805 0.279907095 0.029594762 [171] 0.280978513 0.081706467 0.223438844 0.074266208 0.237226491 [176] 0.048711727 0.268341711 0.102393169 0.200691577 0.151928747 [181] 0.225144022 0.118518740 0.359618251 0.001605944 0.251031001 [186] 0.225801929 0.146982891 0.255664503 0.097423452 0.122734649 [191] 0.221573408 0.160725448 0.072581633 0.243685373 0.043022553 [196] 0.263782209 0.188450195 0.173347223 0.123264484 0.318491791 [201] 0.093839607 0.352131079 0.023243599 0.204111746 0.109682334 [206] 0.299364164 -0.004855058 0.246831257 0.001202806 0.297887720 [211] 0.139234166 0.230364099 0.096717395 0.317571508 0.106408553 [216] 0.347098513 0.034758845 0.245721050 0.112742807 0.273865288 [221] 0.021753925 0.241647107 -0.020607192 0.357432362 0.136285425 [226] 0.174100237 0.140171033 0.266326114 0.138264726 0.308256477 [231] 0.075815857 0.161834741 0.179733813 0.216834741 0.079247246 [236] 0.251705489 0.004008246 0.322996794 0.125925483 0.234887553 [241] 0.083575728 0.285886428 0.116320566 0.335863462 0.039939799 [246] 0.207490865 0.102138437 0.280851474 0.050321400 0.260555672 [251] -0.021833843 0.334823642 0.203711921 0.177685701 0.102958729 [256] 0.296339341 0.120428912 0.327740732 0.087278361 0.174317827 [261] 0.179920159 0.254610361 0.103155456 0.219247756 0.066009695 [266] 0.236095083 0.192356701 0.172543248 0.117019639 0.221304619 [271] 0.127032255 0.289203511 0.106858037 0.161779851 0.201883977 [276] 0.199597133 0.190981624 0.256807616 0.034876358 0.271274460 [281] 0.156297909 0.172771850 0.114338680 0.217955633 0.074367453 [286] 0.350791021 0.063337928 0.242263618 0.131658001 > postscript(file="/var/wessaorg/rcomp/tmp/2gprd1324395415.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > acf(resid,length(resid)/2, main='Residual Autocorrelation Function') > dev.off() null device 1 > postscript(file="/var/wessaorg/rcomp/tmp/3wgem1324395415.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > pacf(resid,length(resid)/2, main='Residual Partial Autocorrelation Function') > dev.off() null device 1 > postscript(file="/var/wessaorg/rcomp/tmp/42wz11324395415.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > cpgram(resid, main='Residual Cumulative Periodogram') > dev.off() null device 1 > postscript(file="/var/wessaorg/rcomp/tmp/533d21324395415.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > hist(resid, main='Residual Histogram', xlab='values of Residuals') > dev.off() null device 1 > postscript(file="/var/wessaorg/rcomp/tmp/6bj201324395415.ps",horizontal=F,onefile=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() null device 1 > postscript(file="/var/wessaorg/rcomp/tmp/7kl591324395415.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > qqnorm(resid, main='Residual Normal Q-Q Plot') > qqline(resid) > dev.off() null device 1 > ncols <- length(selection[[1]][1,]) > nrows <- length(selection[[2]][,1])-1 > > #Note: the /var/wessaorg/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/wessaorg/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/wessaorg/rcomp/tmp/8g4zz1324395415.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/wessaorg/rcomp/tmp/9evhu1324395415.tab") > > try(system("convert tmp/1iavh1324395415.ps tmp/1iavh1324395415.png",intern=TRUE)) character(0) > try(system("convert tmp/2gprd1324395415.ps tmp/2gprd1324395415.png",intern=TRUE)) character(0) > try(system("convert tmp/3wgem1324395415.ps tmp/3wgem1324395415.png",intern=TRUE)) character(0) > try(system("convert tmp/42wz11324395415.ps tmp/42wz11324395415.png",intern=TRUE)) character(0) > try(system("convert tmp/533d21324395415.ps tmp/533d21324395415.png",intern=TRUE)) character(0) > try(system("convert tmp/6bj201324395415.ps tmp/6bj201324395415.png",intern=TRUE)) character(0) > try(system("convert tmp/7kl591324395415.ps tmp/7kl591324395415.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 2.544 0.562 3.112