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(83.87
+ ,84.23
+ ,84.61
+ ,84.82
+ ,85.04
+ ,85.06
+ ,84.93
+ ,84.98
+ ,85.23
+ ,85.30
+ ,85.33
+ ,85.55
+ ,85.70
+ ,85.88
+ ,86.04
+ ,86.07
+ ,86.31
+ ,86.38
+ ,86.35
+ ,86.55
+ ,86.70
+ ,86.74
+ ,86.85
+ ,86.95
+ ,86.80
+ ,87.01
+ ,87.17
+ ,87.43
+ ,87.66
+ ,87.68
+ ,87.59
+ ,87.65
+ ,87.72
+ ,87.70
+ ,87.71
+ ,87.80
+ ,87.62
+ ,87.84
+ ,88.17
+ ,88.47
+ ,88.58
+ ,88.57
+ ,88.55
+ ,88.68
+ ,88.79
+ ,88.85
+ ,88.95
+ ,89.27
+ ,89.09
+ ,89.42
+ ,89.72
+ ,89.85
+ ,89.96
+ ,90.25
+ ,90.20
+ ,90.27
+ ,90.78
+ ,90.79
+ ,90.98
+ ,91.25
+ ,90.75
+ ,91.01
+ ,91.50
+ ,92.09
+ ,92.56
+ ,92.66
+ ,92.38
+ ,92.38
+ ,92.66
+ ,92.69
+ ,92.59
+ ,92.98
+ ,92.98
+ ,93.15
+ ,93.65
+ ,94.06
+ ,94.24
+ ,94.24
+ ,94.11
+ ,94.16
+ ,94.43
+ ,94.67
+ ,94.60
+ ,95.00
+ ,94.84
+ ,95.26
+ ,95.81
+ ,95.92
+ ,95.85
+ ,95.90
+ ,95.80
+ ,96.00
+ ,96.34
+ ,96.43
+ ,96.48
+ ,96.75
+ ,96.51
+ ,96.69
+ ,97.28
+ ,97.69
+ ,98.08
+ ,98.09
+ ,97.92
+ ,98.06
+ ,98.23
+ ,98.57
+ ,98.53
+ ,98.92
+ ,98.42
+ ,98.73
+ ,99.32
+ ,99.73
+ ,100.00
+ ,100.08
+ ,100.02
+ ,100.26
+ ,100.71
+ ,100.95
+ ,100.75
+ ,101.03
+ ,100.64
+ ,100.93
+ ,101.41
+ ,102.07
+ ,102.42
+ ,102.53
+ ,102.43
+ ,102.60
+ ,102.65
+ ,102.74
+ ,102.82
+ ,103.21
+ ,102.75
+ ,103.09
+ ,103.71
+ ,104.30
+ ,104.58
+ ,104.71
+ ,104.44
+ ,104.57
+ ,104.95
+ ,105.49
+ ,106.03
+ ,106.48
+ ,106.25
+ ,106.70
+ ,107.60
+ ,108.05
+ ,108.72
+ ,109.17
+ ,109.08
+ ,109.04
+ ,109.34
+ ,109.37
+ ,108.96
+ ,108.77
+ ,108.11
+ ,108.67
+ ,109.05
+ ,109.43
+ ,109.62
+ ,109.85
+ ,109.34
+ ,109.65
+ ,109.69
+ ,109.91
+ ,110.09)
> par9 = '1'
> par8 = '0'
> par7 = '0'
> par6 = '2'
> par5 = '12'
> par4 = '1'
> par3 = '1'
> par2 = '1'
> par1 = 'TRUE'
> 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(4) #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.2424187 0.04215668 -0.7906195 0.08248815 0.1713431 -0.08242956
 [2,] 0.2519555 0.00000000 -0.7988225 0.09658820 0.1586200 -0.07511134
 [3,] 0.2483080 0.00000000 -0.8124274 0.12005870 0.1072339  0.00000000
 [4,] 0.2543724 0.00000000 -0.7615930 0.00000000 0.1818585  0.00000000
 [5,]        NA         NA         NA         NA        NA          NA
 [6,]        NA         NA         NA         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

[[2]]
         [,1]    [,2] [,3]    [,4]    [,5]    [,6]
 [1,] 0.00347 0.61082    0 0.52685 0.20894 0.45102
 [2,] 0.00190      NA    0 0.45429 0.24406 0.49066
 [3,] 0.00224      NA    0 0.36630 0.36169      NA
 [4,] 0.00172      NA    0      NA 0.03461      NA
 [5,]      NA      NA   NA      NA      NA      NA
 [6,]      NA      NA   NA      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     sma1    sma2    sma3     sma4
      0.2424  0.0422  -0.7906  0.0825  0.1713  -0.0824
s.e.  0.0817  0.0827   0.0949  0.1301  0.1358   0.1091

sigma^2 estimated as 0.02668:  log likelihood = 56.01,  aic = -98.02

[[3]][[2]]

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

Coefficients:
         ar1     ar2     sma1    sma2    sma3     sma4
      0.2424  0.0422  -0.7906  0.0825  0.1713  -0.0824
s.e.  0.0817  0.0827   0.0949  0.1301  0.1358   0.1091

sigma^2 estimated as 0.02668:  log likelihood = 56.01,  aic = -98.02

[[3]][[3]]

Call:
arima(x = series, order = order, seasonal = seasonal, include.mean = include.mean, 
    fixed = last.arma$next.vector, method = "ML")

Coefficients:
         ar1  ar2     sma1    sma2    sma3     sma4
      0.2520    0  -0.7988  0.0966  0.1586  -0.0751
s.e.  0.0798    0   0.0939  0.1288  0.1357   0.1087

sigma^2 estimated as 0.02671:  log likelihood = 55.88,  aic = -99.76

[[3]][[4]]

Call:
arima(x = series, order = order, seasonal = seasonal, include.mean = include.mean, 
    fixed = last.arma$next.vector, method = "ML")

Coefficients:
         ar1  ar2     sma1    sma2    sma3  sma4
      0.2483    0  -0.8124  0.1201  0.1072     0
s.e.  0.0799    0   0.0941  0.1325  0.1172     0

sigma^2 estimated as 0.02676:  log likelihood = 55.64,  aic = -101.28

[[3]][[5]]
NULL

[[3]][[6]]
NULL


$aic
[1]  -98.01983  -99.76022 -101.28055 -102.44651

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
> postscript(file="/var/www/html/rcomp/tmp/1t9vj1261926007.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 = 167 
Frequency = 1 
  [1]  0.048422353  0.021933958  0.014680268  0.011127190  0.009075384
  [6]  0.007569532  0.006361095  0.005609541  0.005220641  0.004763735
 [11]  0.004358350 -0.043399390 -0.290165383 -0.133989394 -0.135263227
 [16] -0.096884149  0.049420676  0.034201267  0.066892801  0.095750911
 [21] -0.106474182 -0.005026549  0.065444566 -0.112902922 -0.228325238
 [26] -0.006314236 -0.091771832  0.148338025 -0.029969236 -0.021707992
 [31] -0.001221853 -0.053594987 -0.108040222 -0.039664544 -0.036369818
 [36] -0.057911574 -0.169834349  0.009777775  0.095430839  0.116374909
 [41] -0.151709640 -0.017962402  0.073224885  0.007248852 -0.059664763
 [46]  0.037783277  0.038823965  0.161470729 -0.179626864  0.131804371
 [51]  0.039045164 -0.071527338 -0.082821305  0.283252174 -0.050018129
 [56] -0.049972825  0.382464110 -0.114410050  0.131217788  0.073337570
 [61] -0.419031884  0.117032481  0.237404163  0.327069324  0.197514040
 [66] -0.044825538 -0.231017042 -0.049033847  0.115380709 -0.016177525
 [71] -0.188627018  0.239415803  0.164936647 -0.134893016  0.210045092
 [76]  0.078997903 -0.059288067 -0.095773767 -0.014352642 -0.021859380
 [81]  0.032263480  0.208606714 -0.194475156  0.170936793  0.039464036
 [86]  0.152222783  0.125673791 -0.279572823 -0.238202988  0.009454677
 [91]  0.038897407  0.135902537  0.004260414  0.016042281  0.007988412
 [96] -0.057517438  0.009354561 -0.092872115  0.151053051  0.018541257
[101]  0.177946193 -0.130242326 -0.002022613  0.075591748 -0.167803571
[106]  0.281922475 -0.102090469  0.049525944 -0.304511374  0.115287143
[111]  0.052893529  0.040074356  0.059430097  0.003110540  0.087848827
[116]  0.120158214  0.114870970  0.035248465 -0.198225555 -0.041449795
[121] -0.115405597  0.041184250 -0.093637344  0.337598436  0.070389484
[126]  0.027266709  0.019997978  0.013929170 -0.269623556 -0.057563066
[131]  0.167582526  0.007022035 -0.155559625  0.096966697  0.028969767
[136]  0.162711843 -0.021640570  0.073121674 -0.169044558 -0.008995987
[141]  0.125083533  0.286837453  0.519140206 -0.048850843  0.135212971
[146]  0.114361315  0.281789392 -0.122161788  0.392355871  0.278956188
[151] -0.048613816 -0.236558481  0.083883959 -0.253969962 -0.400463929
[156] -0.440189830 -0.130182321  0.293622928 -0.311817491 -0.095862680
[161] -0.148808865  0.114115223 -0.386715266  0.261689137 -0.264862138
[166]  0.019964923  0.137178424
> postscript(file="/var/www/html/rcomp/tmp/2vlgn1261926007.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/3utv81261926007.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/4dkeb1261926007.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/548cy1261926007.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/6j6uk1261926007.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/730c41261926007.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/8wfbz1261926007.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/9kko31261926007.tab") 
> 
> try(system("convert tmp/1t9vj1261926007.ps tmp/1t9vj1261926007.png",intern=TRUE))
character(0)
> try(system("convert tmp/2vlgn1261926007.ps tmp/2vlgn1261926007.png",intern=TRUE))
character(0)
> try(system("convert tmp/3utv81261926007.ps tmp/3utv81261926007.png",intern=TRUE))
character(0)
> try(system("convert tmp/4dkeb1261926007.ps tmp/4dkeb1261926007.png",intern=TRUE))
character(0)
> try(system("convert tmp/548cy1261926007.ps tmp/548cy1261926007.png",intern=TRUE))
character(0)
> try(system("convert tmp/6j6uk1261926007.ps tmp/6j6uk1261926007.png",intern=TRUE))
character(0)
> try(system("convert tmp/730c41261926007.ps tmp/730c41261926007.png",intern=TRUE))
character(0)
> 
> 
> proc.time()
   user  system elapsed 
 20.161   1.415  21.044