R version 2.9.0 (2009-04-17) Copyright (C) 2009 The R Foundation for Statistical Computing ISBN 3-900051-07-0 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(1149822 + ,1086979 + ,1276674 + ,1522522 + ,1742117 + ,1737275 + ,1979900 + ,2061036 + ,1867943 + ,1707752 + ,1298756 + ,1281814 + ,1281151 + ,1164976 + ,1454329 + ,1645288 + ,1817743 + ,1895785 + ,2236311 + ,2295951 + ,2087315 + ,1980891 + ,1465446 + ,1445026 + ,1488120 + ,1338333 + ,1715789 + ,1806090 + ,2083316 + ,2092278 + ,2430800 + ,2424894 + ,2299016 + ,2130688 + ,1652221 + ,1608162 + ,1647074 + ,1479691 + ,1884978 + ,2007898 + ,2208954 + ,2217164 + ,2534291 + ,2560312 + ,2429069 + ,2315077 + ,1799608 + ,1772590 + ,1744799 + ,1659093 + ,2099821 + ,2135736 + ,2427894 + ,2468882 + ,2703217 + ,2766841 + ,2655236 + ,2550373 + ,2052097 + ,1998055 + ,1920748 + ,1876694 + ,2380930 + ,2467402 + ,2770771 + ,2781340 + ,3143926 + ,3172235 + ,2952540 + ,2920877 + ,2384552 + ,2248987 + ,2208616 + ,2178756 + ,2632870 + ,2706905 + ,3029745 + ,3015402 + ,3391414 + ,3507805 + ,3177852 + ,3142961 + ,2545815 + ,2414007 + ,2372578 + ,2332664 + ,2825328 + ,2901478 + ,3263955 + ,3226738 + ,3610786 + ,3709274 + ,3467185 + ,3449646 + ,2802951 + ,2462530 + ,2490645 + ,2561520 + ,3067554 + ,3226951 + ,3546493 + ,3492787 + ,3952263 + ,3932072 + ,3720284 + ,3651555 + ,2914972 + ,2713514 + ,2703997 + ,2591373 + ,3163748 + ,3355137 + ,3613702 + ,3686773 + ,4098716 + ,4063517 + ,3551489 + ,3226663 + ,2656842 + ,2597484 + ,2572399 + ,2596631 + ,3165225 + ,3303145 + ,3698247 + ,3668631 + ,4130433 + ,4131400 + ,3864358 + ,3721110 + ,2892532 + ,2843451 + ,2747502 + ,2668775 + ,3018602 + ,3013392 + ,3393657 + ,3544233 + ,4075832 + ,4032923 + ,3734509 + ,3761285 + ,2970090 + ,2847849 + ,2741680 + ,2830639 + ,3257673 + ,3480085 + ,3843271 + ,3796961 + ,4337767 + ,4243630 + ,3927202 + ,3915296 + ,3087396 + ,2963792 + ,2955792 + ,2829925 + ,3281195 + ,3548011 + ,4059648 + ,3941175 + ,4528594 + ,4433151 + ,4145737 + ,4077132 + ,3198519 + ,3078660 + ,3028202 + ,2858642 + ,3398954 + ,3808883 + ,4175961 + ,4227542 + ,4744616 + ,4608012 + ,4295049 + ,4201144 + ,3353276 + ,3286851 + ,3169889 + ,3051720 + ,3695426 + ,3905501 + ,4296458 + ,4246247 + ,4921849 + ,4821446 + ,4425064 + ,4379099 + ,3472889 + ,3359160 + ,3200944 + ,3153170 + ,3741498 + ,3918719 + ,4403449 + ,4400407 + ,4847473 + ,4716136 + ,4297440 + ,4272253 + ,3271834 + ,3168388 + ,2911748 + ,2720999 + ,3199918 + ,3672623 + ,3892013 + ,3850845 + ,4532467 + ,4484739 + ,4014972 + ,3983758 + ,3158459 + ,3100569 + ,2935404 + ,2855719 + ,3465611 + ,3006985 + ,4095110 + ,4104793 + ,4730788 + ,4642726 + ,4246919 + ,4308117) > par9 = '1' > par8 = '2' > par7 = '1' > par6 = '0' > par5 = '12' > par4 = '1' > par3 = '1' > par2 = '0.3' > 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))) [[1]] [,1] [,2] [,3] [,4] [1,] -0.4394548 0.1137506 0.14622927 -0.8225080 [2,] -0.4396846 0.0000000 0.09223073 -0.7358951 [3,] -0.4485068 0.0000000 0.00000000 -0.6991125 [4,] NA NA NA NA [5,] NA NA NA NA [6,] NA NA NA NA [7,] NA NA NA NA [8,] NA NA NA NA [[2]] [,1] [,2] [,3] [,4] [1,] 0 0.40253 0.22797 0 [2,] 0 NA 0.38441 0 [3,] 0 NA NA 0 [4,] NA NA NA NA [5,] NA NA NA NA [6,] NA NA NA NA [7,] NA NA NA NA [8,] NA NA NA NA [[3]] [[3]][[1]] Call: arima(x = series, order = order, seasonal = seasonal, include.mean = include.mean, method = "ML") Coefficients: ma1 sar1 sar2 sma1 -0.4395 0.1138 0.1462 -0.8225 s.e. 0.0674 0.1356 0.1210 0.1154 sigma^2 estimated as 0.7068: log likelihood = -269.95, aic = 549.91 [[3]][[2]] Call: arima(x = series, order = order, seasonal = seasonal, include.mean = include.mean, method = "ML") Coefficients: ma1 sar1 sar2 sma1 -0.4395 0.1138 0.1462 -0.8225 s.e. 0.0674 0.1356 0.1210 0.1154 sigma^2 estimated as 0.7068: log likelihood = -269.95, aic = 549.91 [[3]][[3]] Call: arima(x = series, order = order, seasonal = seasonal, include.mean = include.mean, fixed = last.arma$next.vector, method = "ML") Coefficients: ma1 sar1 sar2 sma1 -0.4397 0 0.0922 -0.7359 s.e. 0.0679 0 0.1058 0.0836 sigma^2 estimated as 0.7115: log likelihood = -270.26, aic = 548.51 [[3]][[4]] NULL $aic [1] 549.9060 548.5109 547.2926 Warning messages: 1: In arima(series, order = order, seasonal = seasonal, include.mean = include.mean, : some AR parameters were fixed: setting transform.pars = FALSE 2: In arima(series, order = order, seasonal = seasonal, include.mean = include.mean, : some AR parameters were fixed: setting transform.pars = FALSE > postscript(file="/var/www/html/rcomp/tmp/1nume1292017783.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 = 226 Frequency = 1 [1] 0.0379864250 0.0161359186 0.0132666550 0.0131368507 0.0131500894 [6] 0.0108843072 0.0120935757 0.0114573201 0.0080239843 0.0052982733 [11] -0.0007661731 -0.0425484138 -0.2405999791 -0.5956928507 0.8277384948 [16] -0.4595685204 -0.7853315562 0.4729400165 0.9338745900 0.1682908503 [21] 0.0724055809 0.6685817726 -0.3915356056 -0.1955336503 0.5145979785 [26] -0.4852417929 1.1823752218 -1.4023651096 0.0359941483 -0.2893913077 [31] 0.0715581422 -0.7715044033 0.5766937645 0.0777745102 0.4427466949 [36] -0.0969652252 0.2300430587 -0.3818206585 0.6647691847 -0.6049015539 [41] -0.7708376521 -0.6467420250 -0.5663851516 -0.4541025300 0.4104754762 [46] 0.6198178582 0.5344385973 0.2859927008 -0.5958471148 0.6112229512 [51] 0.8137816346 -1.1765250444 -0.1348105528 0.1225577128 -1.1016034797 [56] -0.2309643668 0.4935433544 0.7403023420 1.0515502477 0.2459366298 [61] -0.9482934414 0.7999919473 1.0010124150 -0.2156791517 0.0984134310 [66] -0.1390116356 0.0621208384 -0.1725332508 -0.4900561801 0.8292651773 [71] 1.0298147989 -0.5200996812 -0.4392474085 0.7850345207 -0.1637154192 [76] -0.6263776199 -0.2171295947 -0.4798429944 -0.0882037218 0.4682407325 [81] -0.8998237696 0.3353036613 0.2631440568 -0.4615340691 -0.2957184950 [86] 0.4421725644 -0.0462375630 -0.4840417323 0.0433258853 -0.4161411675 [91] -0.2547407906 0.1174475129 0.1304304491 0.6794702362 0.2968901696 [96] -2.1665829753 -0.3587871650 1.3576088716 0.3473546956 0.4902571319 [101] -0.1413413299 -0.4728146768 0.1361541306 -0.7320403898 0.1144973908 [106] 0.1254448631 -0.4748180407 -0.4780644431 -0.1480220153 -0.7222822671 [111] 0.0150726027 0.4878720549 -0.5985079353 0.4315972592 -0.0054478733 [116] -0.6704965623 -2.2851969699 -3.0581420600 -0.7235766813 0.9404884528 [121] 0.2397004481 0.8330912348 0.5873050807 0.1145611464 0.5725297499 [126] -0.0107435552 0.1765665930 -0.0750572666 0.3633544913 0.2361502457 [131] -1.1213453751 0.4277121296 -0.5455974627 -0.5182290629 -2.0668573890 [136] -2.2291466402 -0.4716877103 0.9273240372 1.1576155311 0.0945570528 [141] 0.2702932208 1.6355487775 0.1149450277 -0.0414321086 -0.6377740033 [146] 0.9573044685 -0.4051930822 0.7354971225 0.1830027295 -0.5487332704 [151] 0.1655092475 -0.5936449165 -0.3181586499 0.5452750854 -0.1596237503 [156] -0.1317047900 0.4531804424 -0.7544602231 -0.5586095800 0.8758433956 [161] 1.3068925123 -0.5258068605 0.2300395056 -0.3498701973 0.1023380862 [166] 0.0361268561 -0.5318734556 -0.1598784432 -0.0518798559 -1.2437335946 [171] -0.0652646517 1.6526722997 0.2761619837 0.6605501104 0.0128350217 [176] -0.5411039223 -0.1644820818 -0.2153846186 -0.0789777050 0.4940705565 [181] -0.3904787529 -0.3928810898 0.7648257004 -0.0606869380 -0.3339717294 [186] -0.4048006366 0.5826299857 0.1038799362 -0.4121192847 0.0875178840 [191] -0.1865070901 -0.0665814025 -0.7719752897 0.1405062824 0.1658642628 [196] -0.5653383742 0.2603233976 0.1368637594 -0.8333266547 -0.6598331185 [201] -0.8097170321 -0.0070897385 -1.1061852256 -0.4715270230 -1.6316079370 [206] -1.7568602775 -1.1947494022 1.6168254262 -0.6731850720 -0.4965199317 [211] 0.8724388662 0.6610888213 -0.5145662111 -0.0523468568 0.2659553348 [216] 0.5324337447 0.0046657525 0.2138251800 0.7533343396 -5.6629310322 [221] 3.3199288153 1.6526880431 1.1236062489 0.4665437888 0.2239286086 [226] 0.8566039672 > postscript(file="/var/www/html/rcomp/tmp/2nume1292017783.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/www/html/rcomp/tmp/3xl4z1292017783.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/www/html/rcomp/tmp/4xl4z1292017783.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/www/html/rcomp/tmp/5xl4z1292017783.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/www/html/rcomp/tmp/6xl4z1292017783.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/www/html/rcomp/tmp/7qulk1292017783.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/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/8q5zy1292017783.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/90eg11292017783.tab") > > try(system("convert tmp/1nume1292017783.ps tmp/1nume1292017783.png",intern=TRUE)) character(0) > try(system("convert tmp/2nume1292017783.ps tmp/2nume1292017783.png",intern=TRUE)) character(0) > try(system("convert tmp/3xl4z1292017783.ps tmp/3xl4z1292017783.png",intern=TRUE)) character(0) > try(system("convert tmp/4xl4z1292017783.ps tmp/4xl4z1292017783.png",intern=TRUE)) character(0) > try(system("convert tmp/5xl4z1292017783.ps tmp/5xl4z1292017783.png",intern=TRUE)) character(0) > try(system("convert tmp/6xl4z1292017783.ps tmp/6xl4z1292017783.png",intern=TRUE)) character(0) > try(system("convert tmp/7qulk1292017783.ps tmp/7qulk1292017783.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 4.312 1.141 9.536