R version 2.15.2 (2012-10-26) -- "Trick or Treat"
Copyright (C) 2012 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i686-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(8.9634
+ ,8.9522
+ ,8.8682
+ ,8.7331
+ ,8.3188
+ ,8.3462
+ ,8.3087
+ ,8.3836
+ ,8.8412
+ ,9.5001
+ ,10.1883
+ ,10.2931
+ ,10.1945
+ ,10.3014
+ ,10.0675
+ ,9.6715
+ ,9.5019
+ ,9.4597
+ ,9.4362
+ ,9.5919
+ ,10.0167
+ ,10.322
+ ,11.166
+ ,11.5454
+ ,11.3712
+ ,11.0723
+ ,10.813
+ ,10.3016
+ ,10.4227
+ ,10.3162
+ ,10.4519
+ ,10.8567
+ ,11.2716
+ ,11.4341
+ ,12.1273
+ ,11.9814
+ ,11.8352
+ ,11.9847
+ ,11.545
+ ,11.5285
+ ,11.5539
+ ,11.622
+ ,11.6578
+ ,11.6767
+ ,11.8752
+ ,13.2643
+ ,14.2297
+ ,14.308
+ ,13.7915
+ ,13.7633
+ ,13.9775
+ ,13.6478
+ ,13.2247
+ ,13.0971
+ ,13.1039
+ ,13.206
+ ,13.7901
+ ,14.6457
+ ,15.5764
+ ,15.6102
+ ,15.8855
+ ,16.0137
+ ,15.6186
+ ,15.384
+ ,15.2751
+ ,15.0912
+ ,14.9222
+ ,15.6231
+ ,16.6737
+ ,17.6805
+ ,19.1919
+ ,19.1711
+ ,18.5658
+ ,18.1285
+ ,16.791
+ ,16.9468
+ ,17.3164
+ ,17.1816
+ ,16.7627
+ ,17.239
+ ,17.8838
+ ,18.9038
+ ,20.0274
+ ,20.0087
+ ,19.6366
+ ,19.8163
+ ,18.8602
+ ,17.9206
+ ,17.6889
+ ,17.84
+ ,17.678
+ ,17.7258
+ ,18.5865
+ ,19.9804
+ ,21.1584
+ ,21.2921
+ ,20.9445
+ ,20.5731
+ ,19.3274
+ ,17.7866
+ ,17.7483
+ ,17.5648
+ ,17.4763
+ ,17.7264
+ ,18.5736
+ ,19.9236
+ ,21.3286
+ ,20.7249
+ ,20.3334
+ ,19.7658
+ ,18.7569
+ ,17.6963
+ ,17.7978
+ ,18.1771
+ ,18.3738
+ ,18.1996
+ ,18.8443
+ ,20.1001
+ ,21.2458
+ ,20.8381
+ ,20.1967
+ ,19.8159
+ ,18.5784
+ ,19.21
+ ,19.3419
+ ,19.12
+ ,19.1563
+ ,18.9783
+ ,20.2913
+ ,22.5439
+ ,23.2821
+ ,22.6191
+ ,22.1599
+ ,21.2766
+ ,19.0846
+ ,18.9096
+ ,18.8095
+ ,20.1164
+ ,20.7762
+ ,20.9044
+ ,22.0026
+ ,23.6401
+ ,25.04
+ ,24.7185
+ ,24.1752
+ ,24.1382
+ ,22.3949
+ ,21.3743
+ ,21.4911
+ ,21.2187
+ ,21.2137
+ ,21.6735
+ ,22.5096
+ ,24.3097
+ ,25.7989
+ ,25.4376
+ ,23.878
+ ,23.6966
+ ,23.3544
+ ,21.1993
+ ,22.0431
+ ,22.0203
+ ,21.886
+ ,21.9771
+ ,23.0759
+ ,24.9859
+ ,26.2614
+ ,26.1127
+ ,25.6296
+ ,25.2926
+ ,22.8146
+ ,22.2974
+ ,22.8868
+ ,22.4612
+ ,22.3165
+ ,22.7319
+ ,23.2692
+ ,24.9432
+ ,27.8272
+ ,27.4059
+ ,26.6232
+ ,26.8779
+ ,25.105
+ ,23.601
+ ,23.5374
+ ,23.5248
+ ,22.9465
+ ,23.6633
+ ,25.5932
+ ,27.7683
+ ,29.4691
+ ,28.3472
+ ,28.3879
+ ,27.9696
+ ,26.0075
+ ,24.2533
+ ,24.4999
+ ,23.8988
+ ,23.6683
+ ,23.9427
+ ,26.0155
+ ,28.9529
+ ,30.302
+ ,29.874
+ ,28.2257
+ ,28.0811
+ ,26.3398
+ ,25.4847
+ ,25.4823
+ ,24.9697
+ ,25.2282
+ ,25.9257
+ ,28.7818
+ ,27.9552
+ ,33.3475
+ ,32.7834
+ ,31.6586
+ ,31.6613
+ ,29.1839
+ ,28.8825
+ ,27.6334
+ ,27.7511
+ ,27.3792
+ ,27.7748
+ ,31.4329
+ ,33.2735
+ ,35.0962
+ ,34.9537
+ ,31.8307
+ ,30.9984
+ ,28.629
+ ,26.4379
+ ,25.4408
+ ,24.6681
+ ,24.0994
+ ,24.6043
+ ,27.2492
+ ,29.5511
+ ,29.8522
+ ,31.6989
+ ,29.6357
+ ,30.5197
+ ,32.7823
+ ,24.9942
+ ,23.5187
+ ,24.0249
+ ,24.5692
+ ,24.402
+ ,26.7089
+ ,31.6874
+ ,32.8801
+ ,32.7906
+ ,30.8785
+ ,30.3024
+ ,28.3679
+ ,25.6578
+ ,25.1598
+ ,24.6143
+ ,24.528
+ ,25.2905
+ ,30.0016
+ ,34.2728
+ ,34.4408
+ ,34.1907
+ ,33.6636
+ ,33.9073
+ ,30.2175
+ ,28.5274
+ ,25.9505
+ ,26.2398
+ ,26.2819
+ ,26.7362
+ ,28.8395
+ ,31.0951
+ ,33.7015
+ ,33.8091
+ ,32.1126
+ ,32
+ ,29.122
+ ,26.8124
+ ,25.4654
+ ,23.8331
+ ,24.714
+ ,28.3288
+ ,29.6391
+ ,32.4542
+ ,33.5657
+ ,33.1856
+ ,33.297
+ ,33.51
+ ,31.3789
+ ,29.4555
+ ,27.2699
+ ,27.2586
+ ,27.8591
+ ,29.6362
+ ,30.9587
+ ,31.8633
+ ,33.8188
+ ,33.7531
+ ,33.6103
+ ,32.9052
+ ,29.5005
+ ,27.3634
+ ,27.2298
+ ,26.5211
+ ,26.5228
+ ,27.2991
+ ,29.1726
+ ,30.297
+ ,32.5287
+ ,32.487
+ ,32.4197
+ ,30.854
+ ,28.6995
+ ,27.7881
+ ,26.5609
+ ,25.9431
+ ,25.5578
+ ,27.1275
+ ,30.2556
+ ,34.0976
+ ,34.5614
+ ,34.2948
+ ,33.3418
+ ,31.8187
+ ,29.0818
+ ,27.3444
+ ,26.6233
+ ,26.1869
+ ,26.2953
+ ,28.7043
+ ,32.0653
+ ,34.5401
+ ,34.6636
+ ,34.2557
+ ,32.0526
+ ,30.6892
+ ,28.012
+ ,26.1528
+ ,23.2276
+ ,24.244
+ ,24.8141
+ ,27.8632
+ ,29.6233
+ ,32.4245
+ ,33.3417
+ ,33.0442
+ ,32.0526
+ ,30.2182
+ ,28.9292
+ ,26.8221
+ ,26.1032
+ ,25.9792
+ ,27.1443
+ ,29.4993
+ ,31.656
+ ,33.3665
+ ,35.0521
+ ,34.4076
+ ,33.069
+ ,31.5816
+ ,30.0695
+ ,29.0035
+ ,28.6813
+ ,28.359
+ ,30.0447
+ ,31.5073
+ ,34.16
+ ,35.57
+ ,36.42
+ ,35.12
+ ,33.14
+ ,30.29
+ ,28.2
+ ,26.5
+ ,25.47
+ ,24.96
+ ,25.6
+ ,27.76
+ ,30.13
+ ,32.35
+ ,32.8
+ ,32.54
+ ,29.78
+ ,28.79
+ ,26.8
+ ,25.41
+ ,24.34
+ ,24.39
+ ,25
+ ,26.27
+ ,27.88
+ ,29.35
+ ,29.83
+ ,29.46)
> par9 = '1'
> par8 = '2'
> par7 = '1'
> par6 = '3'
> par5 = '12'
> par4 = '1'
> par3 = '1'
> par2 = '-0.5'
> 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
> 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)))
Error in optim(init[mask], armafn, method = optim.method, hessian = TRUE,  : 
  non-finite finite-difference value [3]
Calls: arimaSelect -> arima -> optim
Execution halted