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(87.28 + ,87.09 + ,86.92 + ,87.59 + ,90.72 + ,90.69 + ,90.3 + ,89.55 + ,88.94 + ,88.41 + ,87.82 + ,87.07 + ,86.82 + ,86.4 + ,86.02 + ,85.66 + ,85.32 + ,85.00 + ,84.67 + ,83.94 + ,82.83 + ,81.95 + ,81.19 + ,80.48 + ,78.86 + ,69.47 + ,68.77 + ,70.06 + ,73.95 + ,75.8 + ,77.79 + ,81.57 + ,83.07 + ,84.34 + ,85.1 + ,85.25 + ,84.26 + ,83.63 + ,86.44 + ,85.3 + ,84.1 + ,83.36 + ,82.48 + ,81.58 + ,80.47 + ,79.34 + ,82.13 + ,81.69 + ,80.7 + ,79.88 + ,79.16 + ,78.38 + ,77.42 + ,76.47 + ,75.46 + ,74.48 + ,78.27 + ,80.7 + ,79.91 + ,78.75 + ,77.78 + ,81.14 + ,81.08 + ,80.03 + ,78.91 + ,78.01 + ,76.9 + ,75.97 + ,81.93 + ,80.27 + ,78.67 + ,77.42 + ,76.16 + ,74.7 + ,76.39 + ,76.04 + ,74.65 + ,73.29 + ,71.79 + ,74.39 + ,74.91 + ,74.54 + ,73.08 + ,72.75 + ,71.32 + ,70.38 + ,70.35 + ,70.01 + ,69.36 + ,67.77 + ,69.26 + ,69.8 + ,68.38 + ,67.62 + ,68.39 + ,66.95 + ,65.21 + ,66.64 + ,63.45 + ,60.66 + ,62.34 + ,60.32 + ,58.64 + ,60.46 + ,58.59 + ,61.87 + ,61.85 + ,67.44 + ,77.06 + ,91.74 + ,93.15 + ,94.15 + ,93.11 + ,91.51 + ,89.96 + ,88.16 + ,86.98 + ,88.03 + ,86.24 + ,84.65 + ,83.23 + ,81.7 + ,80.25 + ,78.8 + ,77.51 + ,76.2 + ,75.04 + ,74.00 + ,75.49 + ,77.14 + ,76.15 + ,76.27 + ,78.19 + ,76.49 + ,77.31 + ,76.65 + ,74.99 + ,73.51 + ,72.07 + ,70.59 + ,71.96 + ,76.29 + ,74.86 + ,74.93 + ,71.9 + ,71.01 + ,77.47 + ,75.78 + ,76.6 + ,76.07 + ,74.57 + ,73.02 + ,72.65 + ,73.16 + ,71.53 + ,69.78 + ,67.98 + ,69.96 + ,72.16 + ,70.47 + ,68.86 + ,67.37 + ,65.87 + ,72.16 + ,71.34 + ,69.93 + ,68.44 + ,67.16 + ,66.01 + ,67.25 + ,70.91 + ,69.75 + ,68.59 + ,67.48 + ,66.31 + ,64.81 + ,66.58 + ,65.97 + ,64.7 + ,64.7 + ,60.94 + ,59.08 + ,58.42 + ,57.77 + ,57.11 + ,53.31 + ,49.96 + ,49.4 + ,48.84 + ,48.3 + ,47.74 + ,47.24 + ,46.76 + ,46.29 + ,48.9 + ,49.23 + ,48.53 + ,48.03 + ,54.34 + ,53.79 + ,53.24 + ,52.96 + ,52.17 + ,51.7 + ,58.55 + ,78.2 + ,77.03 + ,76.19 + ,77.15 + ,75.87 + ,95.47 + ,109.67 + ,112.28 + ,112.01 + ,107.93 + ,105.96 + ,105.06 + ,102.98 + ,102.2 + ,105.23 + ,101.85 + ,99.89 + ,96.23 + ,94.76 + ,91.51 + ,91.63 + ,91.54 + ,85.23 + ,87.83 + ,87.38 + ,84.44 + ,85.19 + ,84.03 + ,86.73 + ,102.52 + ,104.45 + ,106.98 + ,107.02 + ,99.26 + ,94.45 + ,113.44 + ,157.33 + ,147.38 + ,171.89 + ,171.95 + ,132.71 + ,126.02 + ,121.18 + ,115.45 + ,110.48 + ,117.85 + ,117.63 + ,124.65 + ,109.59 + ,111.27 + ,99.78 + ,98.21 + ,99.2 + ,97.97 + ,89.55 + ,87.91 + ,93.34 + ,94.42 + ,93.2 + ,90.29 + ,91.46 + ,89.98 + ,88.35 + ,88.41 + ,82.44 + ,79.89 + ,75.69 + ,75.66 + ,84.5 + ,96.73 + ,87.48 + ,82.39 + ,83.48 + ,79.31 + ,78.16 + ,72.77 + ,72.45 + ,68.46 + ,67.62 + ,68.76 + ,70.07 + ,68.55 + ,65.3 + ,58.96 + ,59.17 + ,62.37 + ,66.28 + ,55.62 + ,55.23 + ,55.85 + ,56.75 + ,50.89 + ,53.88 + ,52.95 + ,55.08 + ,53.61 + ,58.78 + ,61.85 + ,55.91 + ,53.32 + ,46.41 + ,44.57 + ,50.00 + ,50.00 + ,53.36 + ,46.23 + ,50.45 + ,49.07 + ,45.85 + ,48.45 + ,49.96 + ,46.53 + ,50.51 + ,47.58 + ,48.05 + ,46.84 + ,47.67 + ,49.16 + ,55.54 + ,55.82 + ,58.22 + ,56.19 + ,57.77 + ,63.19 + ,54.76 + ,55.74 + ,62.54 + ,61.39 + ,69.6 + ,79.23 + ,80.00 + ,93.68 + ,107.63 + ,100.18 + ,97.3 + ,90.45 + ,80.64 + ,80.58 + ,75.82 + ,85.59 + ,89.35 + ,89.42 + ,104.73 + ,95.32 + ,89.27 + ,90.44 + ,86.97 + ,79.98 + ,81.22 + ,87.35 + ,83.64 + ,82.22 + ,94.4 + ,102.18) > par9 = '0' > par8 = '0' > par7 = '0' > par6 = '3' > par5 = '12' > par4 = '1' > par3 = '1' > par2 = '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 <- as.numeric(par5) #seasonal period > par6 <- as.numeric(par6) #degree (p) of the non-seasonal AR(p) polynomial > par6 <- 3 > par7 <- as.numeric(par7) #degree (q) of the non-seasonal MA(q) polynomial > par7 <- 3 > 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] [,5] [,6] [1,] -0.9313761 -0.83096029 -0.5599887 1.020754064 1.02576504 0.9948550 [2,] 0.0000000 0.06946220 -0.7165064 0.006315295 0.00632898 1.0000052 [3,] 0.0000000 0.07340186 -0.7163801 0.009420422 0.00000000 0.9965082 [4,] 0.0000000 -0.02856430 0.2747690 0.000000000 0.00000000 -0.2007743 [5,] 0.0000000 0.00000000 0.2163869 0.000000000 0.00000000 -0.1415992 [6,] 0.0000000 0.00000000 0.0724933 0.000000000 0.00000000 0.0000000 [7,] 0.0000000 0.00000000 0.0000000 0.000000000 0.00000000 0.0000000 [8,] NA NA NA NA NA NA [9,] NA NA NA NA NA NA [10,] NA NA NA NA NA NA [11,] NA NA NA NA NA NA [12,] NA NA NA NA NA NA [[2]] [,1] [,2] [,3] [,4] [,5] [,6] [1,] 0 0.00000 0.00000 0.00000 0.00000 0.00000 [2,] NA 0.12399 0.00000 0.69687 0.69777 0.00000 [3,] NA 0.09498 0.00000 0.56643 NA 0.00000 [4,] NA 0.60489 0.21978 NA NA 0.38488 [5,] NA NA 0.40863 NA NA 0.58741 [6,] NA NA 0.17782 NA NA NA [7,] NA NA NA NA NA NA [8,] NA NA NA NA NA NA [9,] NA NA NA NA NA NA [10,] NA NA NA NA NA NA [11,] NA NA NA NA NA NA [12,] NA NA NA NA NA NA [[3]] [[3]][[1]] Call: arima(x = series, order = order, seasonal = seasonal, include.mean = include.mean, method = "ML") Coefficients: ar1 ar2 ar3 ma1 ma2 ma3 -0.9314 -0.8310 -0.5600 1.0208 1.0258 0.9949 s.e. 0.0548 0.0632 0.0484 0.0246 0.0230 0.0181 sigma^2 estimated as 49.56: log likelihood = -1170.97, aic = 2355.95 [[3]][[2]] Call: arima(x = series, order = order, seasonal = seasonal, include.mean = include.mean, method = "ML") Coefficients: ar1 ar2 ar3 ma1 ma2 ma3 -0.9314 -0.8310 -0.5600 1.0208 1.0258 0.9949 s.e. 0.0548 0.0632 0.0484 0.0246 0.0230 0.0181 sigma^2 estimated as 49.56: log likelihood = -1170.97, aic = 2355.95 [[3]][[3]] Call: arima(x = series, order = order, seasonal = seasonal, include.mean = include.mean, fixed = last.arma$next.vector, method = "ML") Coefficients: ar1 ar2 ar3 ma1 ma2 ma3 0 0.0695 -0.7165 0.0063 0.0063 1.0000 s.e. 0 0.0450 0.0413 0.0162 0.0163 0.0238 sigma^2 estimated as 53.97: log likelihood = -1185.47, aic = 2382.93 [[3]][[4]] Call: arima(x = series, order = order, seasonal = seasonal, include.mean = include.mean, fixed = last.arma$next.vector, method = "ML") Coefficients: ar1 ar2 ar3 ma1 ma2 ma3 0 0.0734 -0.7164 0.0094 0 0.9965 s.e. 0 0.0438 0.0420 0.0164 0 0.0214 sigma^2 estimated as 54.12: log likelihood = -1185.56, aic = 2381.12 [[3]][[5]] Call: arima(x = series, order = order, seasonal = seasonal, include.mean = include.mean, fixed = last.arma$next.vector, method = "ML") Coefficients: ar1 ar2 ar3 ma1 ma2 ma3 0 -0.0286 0.2748 0 0 -0.2008 s.e. 0 0.0552 0.2235 0 0 0.2308 sigma^2 estimated as 63.84: log likelihood = -1210.01, aic = 2428.03 [[3]][[6]] Call: arima(x = series, order = order, seasonal = seasonal, include.mean = include.mean, fixed = last.arma$next.vector, method = "ML") Coefficients: ar1 ar2 ar3 ma1 ma2 ma3 0 0 0.2164 0 0 -0.1416 s.e. 0 0 0.2616 0 0 0.2607 sigma^2 estimated as 63.88: log likelihood = -1210.14, aic = 2426.27 [[3]][[7]] Call: arima(x = series, order = order, seasonal = seasonal, include.mean = include.mean, fixed = last.arma$next.vector, method = "ML") Coefficients: ar1 ar2 ar3 ma1 ma2 ma3 0 0 0.0725 0 0 0 s.e. 0 0 0.0537 0 0 0 sigma^2 estimated as 63.92: log likelihood = -1210.24, aic = 2424.48 $aic [1] 2355.947 2382.934 2381.121 2428.026 2426.272 2424.476 2424.294 There were 11 warnings (use warnings() to see them) > postscript(file="/var/www/html/rcomp/tmp/1onmp1261123784.ps",horizontal=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 = 359 Frequency = 1 [1] 5.039112e-02 2.238841e-02 1.451297e-02 1.140820e-02 1.193177e-02 [6] 9.900212e-03 8.114944e-03 6.391791e-03 5.102235e-03 4.086515e-03 [11] 3.150427e-03 -4.813126e-02 -3.109582e-01 -2.293197e-01 -2.093654e-01 [16] -1.027195e+00 -3.453327e+00 -2.747764e-01 1.346681e-01 2.715518e-01 [21] -4.789769e-01 -3.543496e-01 -1.714499e-01 7.624665e-02 -1.344627e+00 [26] -8.957676e+00 -3.228997e-01 1.749316e+00 4.880265e+00 2.193198e+00 [31] 2.200386e+00 4.203353e+00 2.452690e+00 1.981816e+00 1.193055e+00 [36] 6.707925e-01 4.741394e-01 8.649810e+00 3.447656e+00 -2.475671e+00 [41] -5.725041e+00 -2.844452e+00 -2.693841e+00 -4.311009e+00 -2.422242e+00 [46] -2.191944e+00 2.369269e+00 -4.007925e-01 1.739839e-01 -3.371614e-01 [51] -3.487229e+00 3.600000e-01 2.537737e-01 4.590137e-02 -1.560976e-01 [56] -9.739839e-02 4.915224e+00 3.569424e+00 -3.574201e+00 -1.075217e+00 [61] -2.380762e-01 4.439526e+00 7.121952e-01 -2.714499e-01 -4.630220e-01 [66] 2.154418e-03 -8.042681e-02 6.159893e-02 2.166375e+00 -4.082751e+00 [71] -8.136247e-01 -2.473105e-01 6.497623e-03 -4.761280e+00 1.756524e+00 [76] 7.210231e-01 7.941774e-02 -5.868633e-01 -4.407453e-01 3.549573e+00 [81] -5.406653e+00 1.318272e+00 -1.159014e-01 1.314364e+00 -2.635164e-01 [86] 5.098509e-01 -1.786694e+00 2.232386e-02 7.023035e-01 -1.053115e-01 [91] 2.989275e+00 -2.113645e+00 -1.923327e+00 -6.067550e-01 2.379336e+00 [96] -9.693630e-01 -2.817276e-01 2.208340e+00 -3.079532e+00 -2.427527e+00 [101] 2.158191e+00 -2.009212e-01 -2.992391e+00 1.111091e+00 -4.188279e-01 [106] 4.269804e+00 -8.827914e-01 7.062622e+00 1.106713e+01 1.330727e+01 [111] 4.090372e+00 2.966476e+00 -3.680536e+00 8.653079e-02 -1.447496e-01 [116] -3.422818e+00 6.595528e-01 -2.239424e+00 -1.507574e+00 -7.230020e+00 [121] -1.087834e+01 -1.608169e+01 -2.339498e+00 -1.649674e+00 9.251165e-01 [126] 4.973309e-01 5.676086e-01 7.781233e-01 2.648977e+00 5.717276e-01 [131] 7.449051e-01 1.516443e+00 3.296504e+00 -2.279946e-01 2.146036e+00 [136] 5.478724e-01 -3.576761e-01 -3.345598e-01 -3.372697e-01 -4.131775e-01 [141] -1.076761e-01 2.700298e+00 -4.081029e-01 -4.130080e-02 -5.144282e+00 [146] 8.418971e-01 5.643625e+00 -6.711581e-01 2.421280e+00 5.411378e-01 [151] 1.466811e-02 -2.497834e-01 -1.808869e+00 -3.815650e+00 -1.949255e-01 [156] -1.693862e+00 1.506924e+00 2.884499e+00 -4.128062e+00 -8.916677e-02 [161] -2.638056e+00 -6.511785e-01 2.842171e-14 8.016159e+00 -3.804064e-01 [166] -1.920000e+00 -4.283475e-01 5.026220e-01 7.891871e-01 -7.501491e-01 [171] 1.425928e+00 4.828794e-01 5.036450e-01 2.741598e-01 2.915785e-01 [176] -7.822622e+00 2.562453e+00 7.760772e-01 7.847229e-01 1.092242e+00 [181] -2.667995e+00 -3.115949e+00 -4.412791e+00 6.992075e-01 7.247292e-01 [186] -2.376829e+00 -2.216972e+00 9.037533e-01 -2.134993e+00 2.280354e-01 [191] 6.418563e-01 -3.310906e-01 3.274925e+00 1.338530e+00 3.306247e+00 [196] 7.422220e-01 -1.407657e-01 3.062947e+00 9.588957e+00 1.289973e-02 [201] -2.292279e-01 -4.402853e-01 -2.307249e-01 2.927507e-02 7.311152e+00 [206] 2.013667e+01 -3.782175e+00 -1.701376e+00 2.014347e-01 -5.059753e-01 [211] 1.337482e+01 1.462966e+01 3.216545e+00 -9.534360e-01 -4.359276e+00 [216] -1.729079e+00 -7.750725e+00 -2.149150e+01 4.987400e-01 4.431823e+00 [221] -2.764720e+00 -7.082724e-01 -2.354055e+01 -1.535538e+01 -5.810705e+00 [226] 2.076194e+00 5.125970e+00 -3.915189e+00 3.471728e+00 1.340752e+00 [231] -1.845379e+00 -2.533727e+00 2.101836e+00 4.816586e+00 1.961528e+01 [236] 3.239065e+00 5.442181e+00 -1.489995e+00 -7.916477e+00 1.080989e+00 [241] 1.639580e+01 4.489602e+01 -7.118740e+00 2.257183e+01 -1.994353e+00 [246] -4.143182e+01 -2.420244e+01 -6.858442e+00 -5.219631e+00 -3.380350e+00 [251] 1.562078e+01 5.188795e+00 -1.160681e+01 -6.004682e+01 1.129726e+01 [256] -3.513226e+01 2.643480e+00 3.938690e+01 8.069759e+00 -3.461836e+00 [261] 1.173594e+00 1.000419e+01 -6.030474e+00 -1.296498e+00 -1.068393e+01 [266] 1.668598e+01 -3.087507e+00 1.057986e+01 4.534336e-01 -6.730921e+00 [271] -2.034784e+00 4.101836e+00 2.114553e+00 3.505691e+00 1.084408e+01 [276] -8.146714e+00 -2.427202e+00 -8.883004e-01 -2.107879e+00 6.380354e-01 [281] -5.444201e+00 5.845007e+00 -1.474797e+00 3.755089e+00 7.604128e-01 [286] -7.425610e+00 -1.399358e+01 5.915183e+00 -7.041254e-01 1.167830e-01 [291] 6.935040e+00 5.150617e+00 -5.206206e+00 -6.042757e-01 4.243184e+00 [296] 2.122040e+00 -6.994925e+00 1.345806e+00 4.638616e-01 5.887453e+00 [301] 4.748211e+00 4.917229e+00 -5.200140e-01 -1.020304e+01 7.710433e+00 [306] -6.510576e+00 -1.745941e+00 3.944979e+00 6.332656e+00 5.483335e-01 [311] -6.528395e+00 1.665189e+00 6.317748e-02 -7.940542e+00 -6.215110e-01 [316] 7.443476e+00 -2.317812e-01 1.092407e+01 -1.630075e+00 -4.899106e+00 [321] -1.999452e+00 -2.450982e+00 8.979567e+00 2.247717e+00 1.843408e+00 [326] 4.995108e+00 -4.786586e+00 -5.033889e-02 8.442588e+00 -1.207436e+01 [331] 3.904925e+00 5.688434e+00 9.596419e-01 7.096551e+00 7.681117e+00 [336] -5.614350e+00 1.286500e+01 1.095990e+01 -5.013313e+00 -5.431410e+00 [341] -1.310730e+01 -9.870863e-01 -7.166799e-01 -1.067051e+01 1.102004e+01 [346] -4.374607e+00 -8.721977e+00 1.374837e+01 -2.276740e+01 -1.930696e+01 [351] 7.565947e+00 1.083870e+00 1.309866e+00 1.042511e+01 6.232771e+00 [356] 1.060149e+00 -1.199105e+01 7.971266e+00 7.633882e+00 > postscript(file="/var/www/html/rcomp/tmp/2yos41261123784.ps",horizontal=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/32ywd1261123784.ps",horizontal=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/4yxg91261123784.ps",horizontal=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/5y8tv1261123784.ps",horizontal=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/6uz1g1261123784.ps",horizontal=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/70pqw1261123784.ps",horizontal=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/8s5981261123784.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/9zpy21261123784.tab") > > try(system("convert tmp/1onmp1261123784.ps tmp/1onmp1261123784.png",intern=TRUE)) character(0) > try(system("convert tmp/2yos41261123784.ps tmp/2yos41261123784.png",intern=TRUE)) character(0) > try(system("convert tmp/32ywd1261123784.ps tmp/32ywd1261123784.png",intern=TRUE)) character(0) > try(system("convert tmp/4yxg91261123784.ps tmp/4yxg91261123784.png",intern=TRUE)) character(0) > try(system("convert tmp/5y8tv1261123784.ps tmp/5y8tv1261123784.png",intern=TRUE)) character(0) > try(system("convert tmp/6uz1g1261123784.ps tmp/6uz1g1261123784.png",intern=TRUE)) character(0) > try(system("convert tmp/70pqw1261123784.ps tmp/70pqw1261123784.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 6.734 1.105 8.893