R version 2.6.1 (2007-11-26)
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(13.5,16.2,17.6,15.8,17.6,15.2,15.9,12.0,13.3,14.8,16.1,16.9,17.6,13.9,10.0,7.6,7.1,8.1,8.1,7.7,4.0,1.4,0.3,-1.0,-1.9,-1.5,-0.2,3.4,3.0,4.1,3.4,3.2,6.1,5.8,6.2,5.8,5.9,6.7,5.9,3.8,1.7,1.4,1.8,3.0,3.6,4.8,4.3,4.2,2.9,4.9,7.2,8.7,9.1,8.9,9.0,11.6,9.6,9.1,9.2,10.8,11.0,8.5,6.5,7.2,7.8,8.7,7.8,7.5,7.7,7.5,8.3,7.9,10.4,11.5,14.0,11.9,11.9,10.3,11.3,9.9,8.9,9.2,8.8,6.7,7.1,6.6,7.2,5.0,5.3,6.3)
> par10 = 'FALSE'
> par9 = '0'
> par8 = '0'
> par7 = '0'
> par6 = '3'
> par5 = '12'
> par4 = '0'
> par3 = '1'
> par2 = '1'
> par1 = '12'
> #'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!)
> par1 <- as.numeric(par1) #cut off periods
> 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
> par7 <- as.numeric(par7) #q
> 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
      0.2241  0.1562  -0.3319
s.e.  0.1089  0.1102   0.1105

sigma^2 estimated as 2.256:  log likelihood = -140.79,  aic = 289.57
> (forecast <- predict(arima.out,fx))
$pred
Time Series:
Start = 79 
End = 90 
Frequency = 1 
 [1] 10.63847 10.46433 11.00923 10.99178 11.13080 10.97837 10.97173 10.90029
 [9] 10.93383 10.93239 10.96102 10.95608

$se
Time Series:
Start = 79 
End = 90 
Frequency = 1 
 [1] 1.501947 2.373998 3.201896 3.659582 4.003948 4.251554 4.504411 4.753816
 [9] 5.015974 5.266227 5.505984 5.728613

> (lb <- forecast$pred - 1.96 * forecast$se)
Time Series:
Start = 79 
End = 90 
Frequency = 1 
 [1]  7.6946588  5.8112990  4.7335155  3.8189996  3.2830604  2.6453262
 [7]  2.1430853  1.5828095  1.1025239  0.6105853  0.1692917 -0.2720049
> (ub <- forecast$pred + 1.96 * forecast$se)
Time Series:
Start = 79 
End = 90 
Frequency = 1 
 [1] 13.58229 15.11737 17.28495 18.16456 18.97854 19.31142 19.80038 20.21777
 [9] 20.76514 21.25420 21.75275 22.18416
> 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] 13.50000 16.20000 17.60000 15.80000 17.60000 15.20000 15.90000 12.00000
 [9] 13.30000 14.80000 16.10000 16.90000 17.60000 13.90000 10.00000  7.60000
[17]  7.10000  8.10000  8.10000  7.70000  4.00000  1.40000  0.30000 -1.00000
[25] -1.90000 -1.50000 -0.20000  3.40000  3.00000  4.10000  3.40000  3.20000
[33]  6.10000  5.80000  6.20000  5.80000  5.90000  6.70000  5.90000  3.80000
[41]  1.70000  1.40000  1.80000  3.00000  3.60000  4.80000  4.30000  4.20000
[49]  2.90000  4.90000  7.20000  8.70000  9.10000  8.90000  9.00000 11.60000
[57]  9.60000  9.10000  9.20000 10.80000 11.00000  8.50000  6.50000  7.20000
[65]  7.80000  8.70000  7.80000  7.50000  7.70000  7.50000  8.30000  7.90000
[73] 10.40000 11.50000 14.00000 11.90000 11.90000 10.30000 10.63847 10.46433
[81] 11.00923 10.99178 11.13080 10.97837 10.97173 10.90029 10.93383 10.93239
[89] 10.96102 10.95608
> (perc.se <- (ub-forecast$pred)/1.96/forecast$pred)
Time Series:
Start = 79 
End = 90 
Frequency = 1 
 [1] 0.1411807 0.2268656 0.2908374 0.3329380 0.3597180 0.3872664 0.4105470
 [8] 0.4361184 0.4587571 0.4817087 0.5023240 0.5228708
> postscript(file="/var/www/html/freestat/rcomp/tmp/1v1u21201180408.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.se <- array(0, dim=fx)
> perf.mse <- 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.mape[i] = perf.mape[i] + abs(perf.pe[i])
+ perf.se[i] = (x[nx+i] - forecast$pred[i])^2
+ perf.mse[i] = perf.mse[i] + perf.se[i]
+ 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 = perf.mape / fx
> perf.mse = perf.mse / fx
> perf.rmse = sqrt(perf.mse)
> postscript(file="/var/www/html/freestat/rcomp/tmp/2tszh1201180408.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:12] <- x[(nx+1):lx]
> lines(dum, lty=1)
> lines(ub,lty=3)
> lines(lb,lty=3)
> dev.off()
null device 
          1 
> load(file='/var/www/html/freestat/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/freestat/rcomp/tmp/3qpfk1201180408.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.mape[i],4))
+ a<-table.element(a,round(perf.se[i],4))
+ a<-table.element(a,round(perf.mse[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/freestat/rcomp/tmp/422vk1201180408.tab") 
> 
> system("convert tmp/1v1u21201180408.ps tmp/1v1u21201180408.png")
> system("convert tmp/2tszh1201180408.ps tmp/2tszh1201180408.png")
> 
> 
> proc.time()
   user  system elapsed 
  1.206   0.441   1.259