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.

  Natural language support but running in an English locale

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,75,74,91,101,103,106,102,105,105,100,95,96,98,99,92,84,81,72,89,96,91,88,90,98,87,100,100,104,107,105,102,98,106,97,101,100,93,94,96,96,98,102,95,85,84,82,87,77,90,90,94,97,96,93,93,93,97,100,95,97,103,102,93,99,100,97,104,102,103,100,90,90)
> par9 = '1'
> par8 = '2'
> par7 = '1'
> par6 = '3'
> par5 = '4'
> par4 = '0'
> par3 = '1'
> 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, 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.7736401 -0.09105228 -0.03848091 -0.9999598 0.2103627 -0.2469827
[2,] 0.7775099 -0.11384910  0.00000000 -1.0000012 0.2029542 -0.2446342
[3,] 0.7940361 -0.11646951  0.00000000 -0.9999969 0.0000000 -0.2470532
[4,] 0.7989715 -0.13530752  0.00000000 -1.0000004 0.0000000 -0.2501875
[5,] 0.7121531  0.00000000  0.00000000 -0.9999924 0.0000000 -0.2138747
[6,] 0.7158433  0.00000000  0.00000000 -1.0000002 0.0000000  0.0000000
[7,]        NA          NA          NA         NA        NA         NA
            [,7]
[1,] -0.24489937
[2,] -0.25495471
[3,] -0.06733045
[4,]  0.00000000
[5,]  0.00000000
[6,]  0.00000000
[7,]          NA

[[2]]
     [,1]    [,2]    [,3] [,4]    [,5]    [,6]    [,7]
[1,]    0 0.55989 0.78936    0 0.55341 0.06520 0.49415
[2,]    0 0.38568      NA    0 0.56073 0.06889 0.46893
[3,]    0 0.37756      NA    0      NA 0.05479 0.61712
[4,]    0 0.28906      NA    0      NA 0.05242      NA
[5,]    0      NA      NA    0      NA 0.09440      NA
[6,]    0      NA      NA    0      NA      NA      NA
[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.7736  -0.0911  -0.0385  -1.0000  0.2104  -0.2470  -0.2449
s.e.  0.1235   0.1554   0.1435   0.0442  0.3531   0.1317   0.3562

sigma^2 estimated as 29.65:  log likelihood = -226.22,  aic = 468.44

[[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.7736  -0.0911  -0.0385  -1.0000  0.2104  -0.2470  -0.2449
s.e.  0.1235   0.1554   0.1435   0.0442  0.3531   0.1317   0.3562

sigma^2 estimated as 29.65:  log likelihood = -226.22,  aic = 468.44

[[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.7775  -0.1138    0  -1.0000  0.2030  -0.2446  -0.255
s.e.  0.1231   0.1304    0   0.0439  0.3471   0.1323   0.350

sigma^2 estimated as 29.71:  log likelihood = -226.26,  aic = 466.51

[[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.7940  -0.1165    0  -1.0000     0  -0.2471  -0.0673
s.e.  0.1193   0.1311    0   0.0425     0   0.1264   0.1340

sigma^2 estimated as 29.92:  log likelihood = -226.43,  aic = 464.87

[[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.7990  -0.1353    0  -1.0000     0  -0.2502     0
s.e.  0.1188   0.1266    0   0.0424     0   0.1267     0

sigma^2 estimated as 30.05:  log likelihood = -226.56,  aic = 463.12

[[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.7122    0    0  -1.0000     0  -0.2139     0
s.e.  0.0876    0    0   0.0404     0   0.1261     0

sigma^2 estimated as 30.74:  log likelihood = -227.12,  aic = 462.24

[[3]][[7]]
NULL


$aic
[1] 468.4421 466.5140 464.8669 463.1191 462.2390 463.0081

Warning messages:
1: In log(s2) : NaNs produced
2: In log(s2) : NaNs produced
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
6: In arima(series, order = order, seasonal = seasonal, include.mean = include.mean,  :
  some AR parameters were fixed: setting transform.pars = FALSE
7: 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/1xrjx1196783158.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 = 73 
Frequency = 1 
 [1]   0.08699995 -10.81559327  -2.40776360  14.15831819  10.58083543
 [6]   4.67006448   5.58297633  -0.77727980   4.57202035  -0.11809026
[11]  -3.36336367  -0.85193406   3.17844344   3.16127998   2.88627801
[16]  -6.04977137  -7.44211798  -4.87878592 -12.38508568  10.76148142
[21]   6.24479119  -3.52595312  -2.92215858  -0.36911796   5.46531269
[26] -10.56384474   8.83576326   4.35727780   7.07573024   4.89891200
[31]   0.79327975   0.08887164  -0.36430936   6.82202193  -3.46686736
[36]   4.94337057   1.86074782  -4.36327258   0.78215764   1.70286292
[41]  -0.13870922   4.13549946   3.48147603  -4.13409943  -9.82123164
[46]  -4.88815275  -4.79533479   1.93390623 -11.83347411   8.74705254
[51]  -0.04397988   1.83290528   0.88656154  -0.94025245  -3.46440855
[56]   0.07585423  -2.80112502   5.48003445   3.57599200  -2.71226783
[61]   2.87973027   6.70880342   0.88714767  -6.89787760   5.52230352
[66]   3.04102113  -0.65636901   6.90878522   1.05606104   4.42058936
[71]  -0.44698853 -10.01183863  -0.16733477
> postscript(file="/var/www/html/rcomp/tmp/2qsx51196783158.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/3k3md1196783158.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/4dmoj1196783158.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/5shsb1196783158.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/6qycu1196783158.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/7fyie1196783158.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/82mb71196783158.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/995471196783158.tab") 
> 
> system("convert tmp/1xrjx1196783158.ps tmp/1xrjx1196783158.png")
> system("convert tmp/2qsx51196783158.ps tmp/2qsx51196783158.png")
> system("convert tmp/3k3md1196783158.ps tmp/3k3md1196783158.png")
> system("convert tmp/4dmoj1196783158.ps tmp/4dmoj1196783158.png")
> system("convert tmp/5shsb1196783158.ps tmp/5shsb1196783158.png")
> system("convert tmp/6qycu1196783158.ps tmp/6qycu1196783158.png")
> system("convert tmp/7fyie1196783158.ps tmp/7fyie1196783158.png")
> 
> 
> proc.time()
   user  system elapsed 
  4.996   1.832   5.338