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(1149822
+ ,1086979
+ ,1276674
+ ,1522522
+ ,1742117
+ ,1737275
+ ,1979900
+ ,2061036
+ ,1867943
+ ,1707752
+ ,1298756
+ ,1281814
+ ,1281151
+ ,1164976
+ ,1454329
+ ,1645288
+ ,1817743
+ ,1895785
+ ,2236311
+ ,2295951
+ ,2087315
+ ,1980891
+ ,1465446
+ ,1445026
+ ,1488120
+ ,1338333
+ ,1715789
+ ,1806090
+ ,2083316
+ ,2092278
+ ,2430800
+ ,2424894
+ ,2299016
+ ,2130688
+ ,1652221
+ ,1608162
+ ,1647074
+ ,1479691
+ ,1884978
+ ,2007898
+ ,2208954
+ ,2217164
+ ,2534291
+ ,2560312
+ ,2429069
+ ,2315077
+ ,1799608
+ ,1772590
+ ,1744799
+ ,1659093
+ ,2099821
+ ,2135736
+ ,2427894
+ ,2468882
+ ,2703217
+ ,2766841
+ ,2655236
+ ,2550373
+ ,2052097
+ ,1998055
+ ,1920748
+ ,1876694
+ ,2380930
+ ,2467402
+ ,2770771
+ ,2781340
+ ,3143926
+ ,3172235
+ ,2952540
+ ,2920877
+ ,2384552
+ ,2248987
+ ,2208616
+ ,2178756
+ ,2632870
+ ,2706905
+ ,3029745
+ ,3015402
+ ,3391414
+ ,3507805
+ ,3177852
+ ,3142961
+ ,2545815
+ ,2414007
+ ,2372578
+ ,2332664
+ ,2825328
+ ,2901478
+ ,3263955
+ ,3226738
+ ,3610786
+ ,3709274
+ ,3467185
+ ,3449646
+ ,2802951
+ ,2462530
+ ,2490645
+ ,2561520
+ ,3067554
+ ,3226951
+ ,3546493
+ ,3492787
+ ,3952263
+ ,3932072
+ ,3720284
+ ,3651555
+ ,2914972
+ ,2713514
+ ,2703997
+ ,2591373
+ ,3163748
+ ,3355137
+ ,3613702
+ ,3686773
+ ,4098716
+ ,4063517
+ ,3551489
+ ,3226663
+ ,2656842
+ ,2597484
+ ,2572399
+ ,2596631
+ ,3165225
+ ,3303145
+ ,3698247
+ ,3668631
+ ,4130433
+ ,4131400
+ ,3864358
+ ,3721110
+ ,2892532
+ ,2843451
+ ,2747502
+ ,2668775
+ ,3018602
+ ,3013392
+ ,3393657
+ ,3544233
+ ,4075832
+ ,4032923
+ ,3734509
+ ,3761285
+ ,2970090
+ ,2847849
+ ,2741680
+ ,2830639
+ ,3257673
+ ,3480085
+ ,3843271
+ ,3796961
+ ,4337767
+ ,4243630
+ ,3927202
+ ,3915296
+ ,3087396
+ ,2963792
+ ,2955792
+ ,2829925
+ ,3281195
+ ,3548011
+ ,4059648
+ ,3941175
+ ,4528594
+ ,4433151
+ ,4145737
+ ,4077132
+ ,3198519
+ ,3078660
+ ,3028202
+ ,2858642
+ ,3398954
+ ,3808883
+ ,4175961
+ ,4227542
+ ,4744616
+ ,4608012
+ ,4295049
+ ,4201144
+ ,3353276
+ ,3286851
+ ,3169889
+ ,3051720
+ ,3695426
+ ,3905501
+ ,4296458
+ ,4246247
+ ,4921849
+ ,4821446
+ ,4425064
+ ,4379099
+ ,3472889
+ ,3359160
+ ,3200944
+ ,3153170
+ ,3741498
+ ,3918719
+ ,4403449
+ ,4400407
+ ,4847473
+ ,4716136
+ ,4297440
+ ,4272253
+ ,3271834
+ ,3168388
+ ,2911748
+ ,2720999
+ ,3199918
+ ,3672623
+ ,3892013
+ ,3850845
+ ,4532467
+ ,4484739
+ ,4014972
+ ,3983758
+ ,3158459
+ ,3100569
+ ,2935404
+ ,2855719
+ ,3465611
+ ,3006985
+ ,4095110
+ ,4104793
+ ,4730788
+ ,4642726
+ ,4246919
+ ,4308117)
> par9 = '1'
> par8 = '2'
> par7 = '1'
> par6 = '3'
> par5 = '12'
> par4 = '1'
> par3 = '1'
> par2 = '0.3'
> 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)))
Error in optim(init[mask], armafn, method = optim.method, hessian = TRUE,  : 
  non-finite finite-difference value [2]
Calls: arimaSelect -> arima -> optim
Execution halted