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(-12.7,-2.4,7.1,-3.9,9.5,5,-16.1,-10.8,7,13.6,8.1,-8.1,4.9,-0.8,4.3,4,1.5,5.4,-11.3,-16.4,-2,8.9,-7.2,-18,1.3,6.3,-6,2.8,2,5.1,-7.6,-18.6,5.8,20.3,0.7,-11.2,-5.7,-0.1,3.4,3.3,-1.2,4.2,-8.8,-25.3,8.5,14.5,-3.1,-10.4,-2.9,0.3,22.6,15.4,9,29.1,2.8,-3.8,27.7,28.9,26.5,19.8,13.2,14.1,34.1,30,21.8,32.1,5.3,3,17.1,26.3,38.1,19.5,38,35.5,78.6,62.2,76.9,104.9,32.2,42.5,64.3,74.9,75.4,43,58.7,55.4,76.6,63.3,78.9,82.7)
> par10 = 'FALSE'
> par9 = '0'
> par8 = '1'
> par7 = '1'
> par6 = '3'
> par5 = '12'
> par4 = '1'
> 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      ma1     sar1
      0.0836  0.2425  0.4332  -0.5842  -0.2914
s.e.  0.2051  0.1420  0.1294   0.1942   0.1557

sigma^2 estimated as 74.54:  log likelihood = -233.27,  aic = 478.55
> (forecast <- predict(arima.out,fx))
$pred
Time Series:
Start = 79 
End = 90 
Frequency = 1 
 [1]  71.42445  80.59235 105.60243 113.08904 127.72840 115.85889 129.27632
 [8] 131.76729 170.42404 159.73691 170.22867 194.78390

$se
Time Series:
Start = 79 
End = 90 
Frequency = 1 
 [1]  8.633862  9.650529 11.386786 14.440116 16.348421 18.662798 21.195507
 [8] 23.407366 25.760120 28.109274 30.348220 32.609833

> (lb <- forecast$pred - 1.96 * forecast$se)
Time Series:
Start = 79 
End = 90 
Frequency = 1 
 [1]  54.50208  61.67731  83.28433  84.78641  95.68549  79.27980  87.73312
 [8]  85.88886 119.93421 104.64273 110.74616 130.86863
> (ub <- forecast$pred + 1.96 * forecast$se)
Time Series:
Start = 79 
End = 90 
Frequency = 1 
 [1]  88.34682  99.50739 127.92053 141.39167 159.77130 152.43797 170.81951
 [8] 177.64573 220.91388 214.83109 229.71118 258.69917
> 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] -12.70000  -2.40000   7.10000  -3.90000   9.50000   5.00000 -16.10000
 [8] -10.80000   7.00000  13.60000   8.10000  -8.10000   4.90000  -0.80000
[15]   4.30000   4.00000   1.50000   5.40000 -11.30000 -16.40000  -2.00000
[22]   8.90000  -7.20000 -18.00000   1.30000   6.30000  -6.00000   2.80000
[29]   2.00000   5.10000  -7.60000 -18.60000   5.80000  20.30000   0.70000
[36] -11.20000  -5.70000  -0.10000   3.40000   3.30000  -1.20000   4.20000
[43]  -8.80000 -25.30000   8.50000  14.50000  -3.10000 -10.40000  -2.90000
[50]   0.30000  22.60000  15.40000   9.00000  29.10000   2.80000  -3.80000
[57]  27.70000  28.90000  26.50000  19.80000  13.20000  14.10000  34.10000
[64]  30.00000  21.80000  32.10000   5.30000   3.00000  17.10000  26.30000
[71]  38.10000  19.50000  38.00000  35.50000  78.60000  62.20000  76.90000
[78] 104.90000  71.42445  80.59235 105.60243 113.08904 127.72840 115.85889
[85] 129.27632 131.76729 170.42404 159.73691 170.22867 194.78390
> (perc.se <- (ub-forecast$pred)/1.96/forecast$pred)
Time Series:
Start = 79 
End = 90 
Frequency = 1 
 [1] 0.1208811 0.1197450 0.1078269 0.1276880 0.1279936 0.1610821 0.1639551
 [8] 0.1776417 0.1511531 0.1759723 0.1782791 0.1674154
> postscript(file="/var/www/html/freestat/rcomp/tmp/15jas1201180873.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/2thbz1201180873.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/34j651201180873.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/4sj5u1201180873.tab") 
> 
> system("convert tmp/15jas1201180873.ps tmp/15jas1201180873.png")
> system("convert tmp/2thbz1201180873.ps tmp/2thbz1201180873.png")
> 
> 
> proc.time()
   user  system elapsed 
  1.863   0.455   1.927