R version 2.6.0 (2007-10-03) Copyright (C) 2007 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(99.5,101.6,103.9,106.6,108.3,102,93.8,91.6,97.7,94.8,98,103.8,97.8,91.2,89.3,87.5,90.4,94.2,102.2,101.3,96,90.8,93.2,90.9,91.1,90.2,94.3,96,99,103.3,113.1,112.8,112.1,107.4,111,110.5,110.8,112.4,111.5,116.2,122.5,121.3,113.9,110.7,120.8,141.1,147.4,148,158.1,165,187,190.3,182.4,168.8,151.2,120.1,112.5,106.2,107.1,108.5,106.5,108.3,125.6,124,127.2,136.9,135.8,124.3,115.4,113.6,114.4,118.4,117,116.5,115.4,113.6,117.4,116.9,116.4,111.1,110.2,118.9,131.8,130.6,138.3,148.4,148.7,144.3,152.5,162.9,167.2,166.5,185.6) > par9 = '1' > par8 = '2' > par7 = '1' > par6 = '3' > par5 = '12' > par4 = '1' > par3 = '1' > par2 = '1' > par1 = 'TRUE' > #'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, ncol=nrc) + pval <- matrix(NA, nrow=nrc, 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) + 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.3124686 -0.03829222 0.1331288 0.2765602 0.005742307 -0.1827434 [2,] 0.3108679 -0.03821530 0.1331422 0.2792477 0.000000000 -0.1836193 [3,] 0.2569435 0.00000000 0.1226424 0.3322907 0.000000000 -0.1824489 [4,] 0.2366643 0.00000000 0.0000000 0.3664889 0.000000000 -0.1776800 [5,] 0.0000000 0.00000000 0.0000000 0.5401783 0.000000000 -0.1726916 [6,] 0.0000000 0.00000000 0.0000000 0.5671852 0.000000000 0.0000000 [7,] NA NA NA NA NA NA [,7] [1,] -0.9999719 [2,] -1.0000826 [3,] -1.0003418 [4,] -1.0000817 [5,] -0.9994337 [6,] -1.0000174 [7,] NA [[2]] [,1] [,2] [,3] [,4] [,5] [,6] [,7] [1,] 0.47788 0.88896 0.31151 0.52802 0.96718 0.15442 0.00951 [2,] 0.47580 0.88876 0.31213 0.51673 NA 0.14621 0.01238 [3,] 0.17449 NA 0.27337 0.07315 NA 0.14834 0.01827 [4,] 0.29206 NA NA 0.10441 NA 0.16064 0.02247 [5,] NA NA NA 0.00000 NA 0.16819 0.00061 [6,] NA NA NA 0.00000 NA NA 0.00001 [7,] NA 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 sar1 sar2 sma1 0.3125 -0.0383 0.1331 0.2766 0.0057 -0.1827 -1.0000 s.e. 0.4383 0.2734 0.1308 0.4365 0.1392 0.1272 0.3769 sigma^2 estimated as 36.1: log likelihood = -271.24, aic = 558.49 [[3]][[2]] Call: arima(x = series, order = order, seasonal = seasonal, include.mean = include.mean, method = "ML") Coefficients: ar1 ar2 ar3 ma1 sar1 sar2 sma1 0.3125 -0.0383 0.1331 0.2766 0.0057 -0.1827 -1.0000 s.e. 0.4383 0.2734 0.1308 0.4365 0.1392 0.1272 0.3769 sigma^2 estimated as 36.1: log likelihood = -271.24, aic = 558.49 [[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 sar1 sar2 sma1 0.3109 -0.0382 0.1331 0.2792 0 -0.1836 -1.0001 s.e. 0.4340 0.2724 0.1310 0.4289 0 0.1252 0.3914 sigma^2 estimated as 36.04: log likelihood = -271.24, aic = 556.49 [[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 sar1 sar2 sma1 0.2569 0 0.1226 0.3323 0 -0.1824 -1.0003 s.e. 0.1877 0 0.1113 0.1832 0 0.1251 0.4159 sigma^2 estimated as 36.05: log likelihood = -271.25, aic = 554.51 [[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 sar1 sar2 sma1 0.2367 0 0 0.3665 0 -0.1777 -1.0001 s.e. 0.2233 0 0 0.2234 0 0.1256 0.4305 sigma^2 estimated as 36.68: log likelihood = -271.85, aic = 553.7 [[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 sar1 sar2 sma1 0 0 0 0.5402 0 -0.1727 -0.9994 s.e. 0 0 0 0.0892 0 0.1243 0.2810 sigma^2 estimated as 37.31: log likelihood = -272.44, aic = 552.87 [[3]][[7]] NULL $aic [1] 558.4850 556.4868 554.5070 553.6989 552.8713 552.6789 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 3: In arima(series, order = order, seasonal = seasonal, include.mean = include.mean, : some AR parameters were fixed: setting transform.pars = FALSE 4: In arima(series, order = order, seasonal = seasonal, include.mean = include.mean, : some AR parameters were fixed: setting transform.pars = FALSE 5: 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/1k3701196518848.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 = 93 Frequency = 1 [1] 5.744633e-02 2.731736e-02 1.982719e-02 1.715940e-02 1.522607e-02 [6] 6.873607e-03 -1.747602e-03 -3.594743e-03 2.576015e-03 -4.407334e-04 [11] 2.657107e-03 -4.722780e-02 -3.590435e-01 -5.332106e+00 -3.799754e-01 [16] -2.910700e+00 2.387418e+00 5.744990e+00 8.182492e+00 -3.514174e+00 [21] -6.045085e+00 1.661621e+00 -1.456418e+00 -4.858110e+00 5.897980e+00 [26] -1.342839e+00 3.952992e+00 -8.695338e-01 9.243226e-01 3.079394e+00 [31] 4.841024e+00 -1.745234e+00 8.531366e-01 -8.074775e-01 1.099881e+00 [36] -1.782164e+00 1.975850e+00 8.605074e-01 -3.156896e+00 4.365236e+00 [41] 8.877868e-01 -1.335891e+00 -7.206780e+00 2.202080e+00 6.179701e+00 [46] 1.712912e+01 -6.671412e+00 2.588437e+00 8.728719e+00 2.326131e+00 [51] 1.775434e+01 -8.280980e+00 -5.607003e+00 -8.436007e+00 -1.000351e+01 [56] -2.029647e+01 1.541313e+00 -9.022811e+00 2.214638e+00 -9.698707e-01 [61] -2.801890e+00 2.742673e+00 8.610189e+00 -7.582973e+00 6.631233e+00 [66] 7.500787e+00 -3.392860e+00 -1.323154e+00 -6.255932e+00 4.839753e+00 [71] -4.394918e+00 5.073525e+00 -3.412800e+00 1.597872e+00 -5.753489e+00 [76] 2.578276e-01 2.636760e-01 -2.290833e+00 7.164618e-01 -1.734837e+00 [81] -3.021373e-01 6.802377e+00 5.111046e+00 -5.173361e+00 8.745823e+00 [86] 4.163083e+00 -5.877816e+00 -2.465420e+00 7.552222e+00 7.969056e+00 [91] 2.366862e+00 4.998432e+00 1.463923e+01 > postscript(file="/var/www/html/rcomp/tmp/2xxfx1196518848.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/3dsha1196518848.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/45wne1196518848.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/59y5s1196518848.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/6219a1196518848.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/733bp1196518848.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > qqnorm(resid, main='Residual Normal Q-Q Plot') > dev.off() null device 1 > ncols <- length(selection[[1]][1,]) > nrows <- length(selection[[2]][,1])-1 > 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/8z2l31196518848.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/9tsd11196518848.tab") > > system("convert tmp/1k3701196518848.ps tmp/1k3701196518848.png") > system("convert tmp/2xxfx1196518848.ps tmp/2xxfx1196518848.png") > system("convert tmp/3dsha1196518848.ps tmp/3dsha1196518848.png") > system("convert tmp/45wne1196518848.ps tmp/45wne1196518848.png") > system("convert tmp/59y5s1196518848.ps tmp/59y5s1196518848.png") > system("convert tmp/6219a1196518848.ps tmp/6219a1196518848.png") > system("convert tmp/733bp1196518848.ps tmp/733bp1196518848.png") > > > proc.time() user system elapsed 8.601 1.278 9.273