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(105.7,105.7,111.1,82.4,60,107.3,99.3,113.5,108.9,100.2,103.9,138.7,120.2,100.2,143.2,70.9,85.2,133,136.6,117.9,106.3,122.3,125.5,148.4,126.3,99.6,140.4,80.3,92.6,138.5,110.9,119.6,105,109,129.4,148.6,101.4,134.8,143.7,81.6,90.3,141.5,140.7,140.2,100.2,125.7,119.6,134.7,109,116.3,146.9,97.4,89.4,132.1,139.8,129,112.5,121.9,121.7,123.1)
> par10 = 'FALSE'
> par9 = '1'
> par8 = '0'
> par7 = '1'
> par6 = '1'
> par5 = '12'
> par4 = '1'
> par3 = '0'
> par2 = '0.5'
> par1 = '24'
> par1 <- as.numeric(par1) #cut off periods
> par1 <- 28
> par2 <- as.numeric(par2) #lambda
> 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) #p
> par6 <- 3
> par7 <- as.numeric(par7) #q
> par7 <- 3
> par8 <- as.numeric(par8) #P
> par9 <- as.numeric(par9) #Q
> if (par10 == 'TRUE') par10 <- TRUE
> if (par10 == 'FALSE') par10 <- FALSE
> if (par2 == 0) x <- log(x)
> if (par2 != 0) x <- x^par2
> lx <- length(x)
> first <- lx - 2*par1
> nx <- lx - par1
> nx1 <- nx + 1
> fx <- lx - nx
> if (fx < 1) {
+ fx <- par5
+ nx1 <- lx + fx - 1
+ first <- lx - 2*fx
+ }
> first <- 1
> if (fx < 3) fx <- round(lx/10,0)
> (arima.out <- arima(x[1:nx], order=c(par6,par3,par7), seasonal=list(order=c(par8,par4,par9), period=par5), include.mean=par10, method='ML'))

Call:
arima(x = x[1:nx], order = c(par6, par3, par7), seasonal = list(order = c(par8, 
    par4, par9), period = par5), include.mean = par10, method = "ML")

Coefficients:
          ar1     ar2     ar3     ma1     ma2      ma3     sma1
      -0.1553  0.5542  0.4992  0.4622  0.0088  -0.7909  -0.4732
s.e.   0.5709  0.2023  0.4706  0.4228  0.3137   0.2923   0.6223

sigma^2 estimated as 0.3075:  log likelihood = -19.58,  aic = 55.16
> (forecast <- predict(arima.out,par1))
$pred
Time Series:
Start = 33 
End = 60 
Frequency = 1 
 [1] 10.340209 11.066383 10.818864 12.294787 11.245004 10.187210 11.888259
 [8]  8.995467  9.352616 11.629040 10.900059 11.004716 10.496027 10.993506
[15] 10.927659 12.315286 11.265738 10.249664 11.900286  9.038565  9.383767
[22] 11.654093 10.934949 11.028734 10.524142 11.019869 10.951138 12.340287

$se
Time Series:
Start = 33 
End = 60 
Frequency = 1 
 [1] 0.5813790 0.6012056 0.6619126 0.6679624 0.7103611 0.7119383 0.7162758
 [8] 0.7285772 0.7300590 0.7376417 0.7415247 0.7448048 0.8449737 0.8575993
[15] 0.8851504 0.8848515 0.9039256 0.9081695 0.9135083 0.9222878 0.9256326
[22] 0.9317910 0.9361495 0.9399474 1.0238102 1.0363723 1.0611679 1.0611167

> (lb <- forecast$pred - 1.96 * forecast$se)
Time Series:
Start = 33 
End = 60 
Frequency = 1 
 [1]  9.200706  9.888020  9.521516 10.985581  9.852696  8.791811 10.484359
 [8]  7.567456  7.921701 10.183262  9.446671  9.544898  8.839879  9.312611
[15]  9.192764 10.580977  9.494044  8.469652 10.109810  7.230881  7.569527
[22]  9.827783  9.100096  9.186437  8.517474  8.988580  8.871248 10.260498
> (ub <- forecast$pred + 1.96 * forecast$se)
Time Series:
Start = 33 
End = 60 
Frequency = 1 
 [1] 11.47971 12.24475 12.11621 13.60399 12.63731 11.58261 13.29216 10.42348
 [9] 10.78353 13.07482 12.35345 12.46453 12.15218 12.67440 12.66255 14.04960
[17] 13.03743 12.02968 13.69076 10.84625 11.19801 13.48040 12.76980 12.87103
[25] 12.53081 13.05116 13.03103 14.42008
> if (par2 == 0) {
+ x <- exp(x)
+ forecast$pred <- exp(forecast$pred)
+ lb <- exp(lb)
+ ub <- exp(ub)
+ }
> if (par2 != 0) {
+ x <- x^(1/par2)
+ forecast$pred <- forecast$pred^(1/par2)
+ lb <- lb^(1/par2)
+ ub <- ub^(1/par2)
+ }
> if (par2 < 0) {
+ olb <- lb
+ lb <- ub
+ ub <- olb
+ }
> (actandfor <- c(x[1:nx], forecast$pred))
 [1] 105.70000 105.70000 111.10000  82.40000  60.00000 107.30000  99.30000
 [8] 113.50000 108.90000 100.20000 103.90000 138.70000 120.20000 100.20000
[15] 143.20000  70.90000  85.20000 133.00000 136.60000 117.90000 106.30000
[22] 122.30000 125.50000 148.40000 126.30000  99.60000 140.40000  80.30000
[29]  92.60000 138.50000 110.90000 119.60000 106.91993 122.46484 117.04782
[36] 151.16179 126.45011 103.77925 141.33071  80.91843  87.47143 135.23457
[43] 118.81130 121.10376 110.16659 120.85717 119.41373 151.66627 126.91686
[50] 105.05561 141.61681  81.69566  88.05508 135.81788 119.57311 121.63297
[57] 110.75756 121.43752 119.92741 152.28268
> (perc.se <- (ub-forecast$pred)/1.96/forecast$pred)
Time Series:
Start = 33 
End = 60 
Frequency = 1 
 [1] 0.1186462 0.1144392 0.1296993 0.1144430 0.1341641 0.1493436 0.1276165
 [8] 0.1748452 0.1680614 0.1347481 0.1451297 0.1443391 0.1737109 0.1679468
[15] 0.1748617 0.1538180 0.1730917 0.1925972 0.1650767 0.2244859 0.2163550
[22] 0.1724375 0.1855868 0.1846911 0.2131132 0.2054270 0.2122042 0.1864681
> postscript(file="/var/www/html/rcomp/tmp/1fbuy1260390868.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) 
> opar <- par(mar=c(4,4,2,2),las=1)
> ylim <- c( min(x[first:nx],lb), max(x[first:nx],ub))
> plot(x,ylim=ylim,type='n',xlim=c(first,lx))
> usr <- par('usr')
> rect(usr[1],usr[3],nx+1,usr[4],border=NA,col='lemonchiffon')
> rect(nx1,usr[3],usr[2],usr[4],border=NA,col='lavender')
> abline(h= (-3:3)*2 , col ='gray', lty =3)
> polygon( c(nx1:lx,lx:nx1), c(lb,rev(ub)), col = 'orange', lty=2,border=NA)
> lines(nx1:lx, lb , lty=2)
> lines(nx1:lx, ub , lty=2)
> lines(x, lwd=2)
> lines(nx1:lx, forecast$pred , lwd=2 , col ='white')
> box()
> par(opar)
> dev.off()
null device 
          1 
> prob.dec <- array(NA, dim=fx)
> prob.sdec <- array(NA, dim=fx)
> prob.ldec <- array(NA, dim=fx)
> prob.pval <- array(NA, dim=fx)
> perf.pe <- array(0, dim=fx)
> perf.mape <- array(0, dim=fx)
> perf.mape1 <- array(0, dim=fx)
> perf.se <- array(0, dim=fx)
> perf.mse <- array(0, dim=fx)
> perf.mse1 <- array(0, dim=fx)
> perf.rmse <- array(0, dim=fx)
> for (i in 1:fx) {
+ locSD <- (ub[i] - forecast$pred[i]) / 1.96
+ perf.pe[i] = (x[nx+i] - forecast$pred[i]) / forecast$pred[i]
+ perf.se[i] = (x[nx+i] - forecast$pred[i])^2
+ prob.dec[i] = pnorm((x[nx+i-1] - forecast$pred[i]) / locSD)
+ prob.sdec[i] = pnorm((x[nx+i-par5] - forecast$pred[i]) / locSD)
+ prob.ldec[i] = pnorm((x[nx] - forecast$pred[i]) / locSD)
+ prob.pval[i] = pnorm(abs(x[nx+i] - forecast$pred[i]) / locSD)
+ }
> perf.mape[1] = abs(perf.pe[1])
> perf.mse[1] = abs(perf.se[1])
> for (i in 2:fx) {
+ perf.mape[i] = perf.mape[i-1] + abs(perf.pe[i])
+ perf.mape1[i] = perf.mape[i] / i
+ perf.mse[i] = perf.mse[i-1] + perf.se[i]
+ perf.mse1[i] = perf.mse[i] / i
+ }
> perf.rmse = sqrt(perf.mse1)
> postscript(file="/var/www/html/rcomp/tmp/2qfmt1260390868.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) 
> plot(forecast$pred, pch=19, type='b',main='ARIMA Extrapolation Forecast', ylab='Forecast and 95% CI', xlab='time',ylim=c(min(lb),max(ub)))
> dum <- forecast$pred
> dum[1:par1] <- x[(nx+1):lx]
> lines(dum, lty=1)
> lines(ub,lty=3)
> lines(lb,lty=3)
> dev.off()
null device 
          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,'Univariate ARIMA Extrapolation Forecast',9,TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a,'time',1,header=TRUE)
> a<-table.element(a,'Y[t]',1,header=TRUE)
> a<-table.element(a,'F[t]',1,header=TRUE)
> a<-table.element(a,'95% LB',1,header=TRUE)
> a<-table.element(a,'95% UB',1,header=TRUE)
> a<-table.element(a,'p-value<br />(H0: Y[t] = F[t])',1,header=TRUE)
> a<-table.element(a,'P(F[t]>Y[t-1])',1,header=TRUE)
> a<-table.element(a,'P(F[t]>Y[t-s])',1,header=TRUE)
> mylab <- paste('P(F[t]>Y[',nx,sep='')
> mylab <- paste(mylab,'])',sep='')
> a<-table.element(a,mylab,1,header=TRUE)
> a<-table.row.end(a)
> for (i in (nx-par5):nx) {
+ a<-table.row.start(a)
+ a<-table.element(a,i,header=TRUE)
+ a<-table.element(a,x[i])
+ a<-table.element(a,'-')
+ a<-table.element(a,'-')
+ a<-table.element(a,'-')
+ a<-table.element(a,'-')
+ a<-table.element(a,'-')
+ a<-table.element(a,'-')
+ a<-table.element(a,'-')
+ a<-table.row.end(a)
+ }
> for (i in 1:fx) {
+ a<-table.row.start(a)
+ a<-table.element(a,nx+i,header=TRUE)
+ a<-table.element(a,round(x[nx+i],4))
+ a<-table.element(a,round(forecast$pred[i],4))
+ a<-table.element(a,round(lb[i],4))
+ a<-table.element(a,round(ub[i],4))
+ a<-table.element(a,round((1-prob.pval[i]),4))
+ a<-table.element(a,round((1-prob.dec[i]),4))
+ a<-table.element(a,round((1-prob.sdec[i]),4))
+ a<-table.element(a,round((1-prob.ldec[i]),4))
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/www/html/rcomp/tmp/3te3j1260390868.tab") 
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a,'Univariate ARIMA Extrapolation Forecast Performance',7,TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a,'time',1,header=TRUE)
> a<-table.element(a,'% S.E.',1,header=TRUE)
> a<-table.element(a,'PE',1,header=TRUE)
> a<-table.element(a,'MAPE',1,header=TRUE)
> a<-table.element(a,'Sq.E',1,header=TRUE)
> a<-table.element(a,'MSE',1,header=TRUE)
> a<-table.element(a,'RMSE',1,header=TRUE)
> a<-table.row.end(a)
> for (i in 1:fx) {
+ a<-table.row.start(a)
+ a<-table.element(a,nx+i,header=TRUE)
+ a<-table.element(a,round(perc.se[i],4))
+ a<-table.element(a,round(perf.pe[i],4))
+ a<-table.element(a,round(perf.mape1[i],4))
+ a<-table.element(a,round(perf.se[i],4))
+ a<-table.element(a,round(perf.mse1[i],4))
+ a<-table.element(a,round(perf.rmse[i],4))
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/www/html/rcomp/tmp/4gy061260390868.tab") 
> 
> system("convert tmp/1fbuy1260390868.ps tmp/1fbuy1260390868.png")
> system("convert tmp/2qfmt1260390868.ps tmp/2qfmt1260390868.png")
> 
> 
> proc.time()
   user  system elapsed 
  1.446   0.357   1.983