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(492,436,694,1137,380,179,2354,111,740,595,809,693,738,1184,713,1729,844,1298,514,689,837,1330,491,622,1332,1043,1082,636,586,1170,973,721,863,343,1278,1186,1334,652,284,1273,1518,715,671,486,1022,2084,330,658,1385,930,620,218,840,255,454,1149,684,1190,1079,883,1331,1159,1217,946,579,474,626,843,893,633,873,385,729,774,769,996,1194,575,725,706,665,1259,653,694,437,822,458,1545,987,1051,838,703,613,1128,967,617,654,805,1355,1456,878,887,663,214,733,830,1174,1068,413,946,657,690,156,779,192,461,1213,146,866,200,1290,715,514,697,276,752,1021,481,1626,884,1187,488,403,977,1525,551,1807,723,632,898,621,1606,811,716,1001,732,1024,831,0,85,0,0,0,0,773,1128,0,0,74,259,69,301,0,668)
> par9 = '0'
> par8 = '2'
> par7 = '0'
> par6 = '3'
> par5 = '12'
> par4 = '0'
> par3 = '0'
> par2 = '1'
> 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]       [,5]
 [1,] 0.2876443 0.3411959 0.2924562 -0.06400004 0.03020621
 [2,] 0.2891454 0.3393500 0.2970033 -0.06736088 0.00000000
 [3,] 0.2813827 0.3399263 0.2950578  0.00000000 0.00000000
 [4,]        NA        NA        NA          NA         NA
 [5,]        NA        NA        NA          NA         NA
 [6,]        NA        NA        NA          NA         NA
 [7,]        NA        NA        NA          NA         NA
 [8,]        NA        NA        NA          NA         NA
 [9,]        NA        NA        NA          NA         NA
[10,]        NA        NA        NA          NA         NA

[[2]]
         [,1]  [,2]    [,3]    [,4]    [,5]
 [1,] 0.00019 1e-05 0.00016 0.47281 0.74755
 [2,] 0.00017 1e-05 0.00010 0.44619      NA
 [3,] 0.00022 1e-05 0.00011      NA      NA
 [4,]      NA    NA      NA      NA      NA
 [5,]      NA    NA      NA      NA      NA
 [6,]      NA    NA      NA      NA      NA
 [7,]      NA    NA      NA      NA      NA
 [8,]      NA    NA      NA      NA      NA
 [9,]      NA    NA      NA      NA      NA
[10,]      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     sar1    sar2
      0.2876  0.3412  0.2925  -0.0640  0.0302
s.e.  0.0752  0.0733  0.0757   0.0889  0.0937

sigma^2 estimated as 203242:  log likelihood = -1235.8,  aic = 2483.61

[[3]][[2]]

Call:
arima(x = series, order = order, seasonal = seasonal, include.mean = include.mean, 
    method = "ML")

Coefficients:
         ar1     ar2     ar3     sar1    sar2
      0.2876  0.3412  0.2925  -0.0640  0.0302
s.e.  0.0752  0.0733  0.0757   0.0889  0.0937

sigma^2 estimated as 203242:  log likelihood = -1235.8,  aic = 2483.61

[[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     sar1  sar2
      0.2891  0.3393  0.2970  -0.0674     0
s.e.  0.0750  0.0732  0.0743   0.0882     0

sigma^2 estimated as 203365:  log likelihood = -1235.86,  aic = 2481.71

[[3]][[4]]
NULL

[[3]][[5]]
NULL


$aic
[1] 2483.607 2481.711 2480.288

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/wessaorg/rcomp/tmp/1xvb51324463784.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 = 164 
Frequency = 1 
  [1]  257.419157   41.721668  269.391296  644.371197 -309.351357 -517.838762
  [7] 1835.583776 -737.416419 -139.167537 -350.193388  355.671118   43.899608
 [13]   91.362136  494.817411  -68.726051  945.123125 -270.676050  220.243119
 [19] -537.593339 -200.833949   68.136751  677.547734 -358.464828 -217.379252
 [25]  596.334408  334.308774  137.904245 -365.655600 -291.656848  480.584805
 [31]  202.404918 -141.577730  -17.911636 -392.933446  646.077090  428.948777
 [37]  495.289836 -495.486715 -699.783726  544.701112  841.382360 -209.050753
 [43] -412.322833 -410.355452  439.854408 1394.628636 -718.483648 -418.268491
 [49]  494.485721  173.488211 -362.123056 -649.588792  348.279663 -262.186895
 [55]    1.572289  654.665413  151.705215  563.412692  110.098716  -66.153112
 [61]  387.334201  168.061451  146.778001 -240.866725 -432.188613 -392.464182
 [67]   13.548651  375.098836  304.253855  -65.713618  147.439473 -349.879016
 [73]  157.410104  183.655317  194.782427  281.372682  384.740402 -361.948227
 [79] -141.483322  -31.204832   63.999086  605.262170 -137.187496 -142.965711
 [85] -350.202044  277.863694 -121.735113 1023.671384  168.677575   82.614725
 [91] -269.282912 -192.698533 -183.827171  504.514140  214.168421 -235.507115
 [97] -211.769150  137.248693  708.019359  664.410300 -232.424149 -256.310211
[103] -341.351658 -552.213843  170.109292  379.730773  636.797860  213.857930
[109] -553.352031  123.512703  -25.581349   96.544406 -563.723134  286.867664
[115] -312.930505   58.458390  795.489511 -394.721614  317.113498 -444.772276
[121]  858.508205   24.706082 -194.877891  -73.594798 -349.192840  303.528106
[127]  483.280993 -144.995897  969.854193  -80.788364  255.293880 -669.128455
[133] -343.178894  343.468995  948.017370 -346.397101  818.962002 -420.333865
[139] -320.031165  -76.972023   -6.078593  930.453951 -114.864892 -290.916899
[145]   14.594241  -18.219240  324.722850  -33.771825 -748.597471 -530.729360
[151] -295.226510  -33.342758  -29.816059   62.914951  764.188332  887.789650
[157] -585.659798 -615.150847 -243.505717  236.876456  -85.238242  137.422449
[163] -205.652642  543.419430
> postscript(file="/var/wessaorg/rcomp/tmp/2rqmt1324463784.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/3j7vj1324463784.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/4mkst1324463784.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/5vqrb1324463784.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/6j7e91324463784.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/7w7qu1324463784.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/873951324463784.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/9zkar1324463784.tab") 
> 
> try(system("convert tmp/1xvb51324463784.ps tmp/1xvb51324463784.png",intern=TRUE))
character(0)
> try(system("convert tmp/2rqmt1324463784.ps tmp/2rqmt1324463784.png",intern=TRUE))
character(0)
> try(system("convert tmp/3j7vj1324463784.ps tmp/3j7vj1324463784.png",intern=TRUE))
character(0)
> try(system("convert tmp/4mkst1324463784.ps tmp/4mkst1324463784.png",intern=TRUE))
character(0)
> try(system("convert tmp/5vqrb1324463784.ps tmp/5vqrb1324463784.png",intern=TRUE))
character(0)
> try(system("convert tmp/6j7e91324463784.ps tmp/6j7e91324463784.png",intern=TRUE))
character(0)
> try(system("convert tmp/7w7qu1324463784.ps tmp/7w7qu1324463784.png",intern=TRUE))
character(0)
> 
> 
> proc.time()
   user  system elapsed 
  4.839   0.989   5.839