R version 2.8.0 (2008-10-20)
Copyright (C) 2008 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(3.4,3.4,3.4,-3.5,-3.5,-3.5,-8.3,-8.3,-8.3,-16.7,-16.7,-16.7,-11.6,-11.6,-11.6,-8.4,-8.4,-8.4,-8.6,-8.6,-8.6,0.6,0.6,0.6,-1.5,-1.5,-1.5,9.3,9.3,9.3,2.0,2.0,2.0,-5.5,-5.5,-5.5,4.0,4.0,4.0,-0.5,-0.5,-0.5,10.9,10.9,10.9,19.4,19.4,19.4,13.9,13.9,13.9,10.6,10.6,10.6,4.8,4.8,4.8,4.7,4.7,4.7,-3.9,-3.9,-3.9,-0.2,-0.2,-0.2)
> par10 = 'FALSE'
> par9 = '1'
> par8 = '0'
> par7 = '0'
> par6 = '0'
> par5 = '12'
> par4 = '0'
> par3 = '1'
> par2 = '1'
> par1 = '15'
> #'GNU S' R Code compiled by R2WASP v. 1.0.44 ()
> #Author: Prof. Dr. P. Wessa
> #To cite this work: Wessa P., (2009), ARIMA Forecasting (v1.0.5) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_arimaforecasting.wasp/
> #Source of accompanying publication: 
> #Technical description: 
> 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:
         sma1
      -0.9994
s.e.   1.4253

sigma^2 estimated as 7.777:  log likelihood = -132.05,  aic = 268.1
> (forecast <- predict(arima.out,par1))
$pred
Time Series:
Start = 52 
End = 66 
Frequency = 1 
 [1] 11.120004 11.120004 11.120004  7.420003  7.420003  7.420003  3.120005
 [8]  3.120005  3.120005  1.640006  1.640006  1.640006  1.640006  1.640006
[15]  1.640006

$se
Time Series:
Start = 52 
End = 66 
Frequency = 1 
 [1]  3.054058  4.319090  5.289783  6.108115  6.829081  7.480883  8.080277
 [8]  8.638180  9.162173  9.657778 10.116366 10.555048 10.555048 10.555048
[15] 10.555048

> (lb <- forecast$pred - 1.96 * forecast$se)
Time Series:
Start = 52 
End = 66 
Frequency = 1 
 [1]   5.134051   2.654588   0.752029  -4.551903  -5.964995  -7.242528
 [7] -12.717339 -13.810827 -14.837855 -17.289239 -18.188071 -19.047888
[13] -19.047888 -19.047888 -19.047889
> (ub <- forecast$pred + 1.96 * forecast$se)
Time Series:
Start = 52 
End = 66 
Frequency = 1 
 [1] 17.10596 19.58542 21.48798 19.39191 20.80500 22.08253 18.95735 20.05084
 [9] 21.07786 20.56925 21.46808 22.32790 22.32790 22.32790 22.32790
> 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]   3.400000   3.400000   3.400000  -3.500000  -3.500000  -3.500000
 [7]  -8.300000  -8.300000  -8.300000 -16.700000 -16.700000 -16.700000
[13] -11.600000 -11.600000 -11.600000  -8.400000  -8.400000  -8.400000
[19]  -8.600000  -8.600000  -8.600000   0.600000   0.600000   0.600000
[25]  -1.500000  -1.500000  -1.500000   9.300000   9.300000   9.300000
[31]   2.000000   2.000000   2.000000  -5.500000  -5.500000  -5.500000
[37]   4.000000   4.000000   4.000000  -0.500000  -0.500000  -0.500000
[43]  10.900000  10.900000  10.900000  19.400000  19.400000  19.400000
[49]  13.900000  13.900000  13.900000  11.120004  11.120004  11.120004
[55]   7.420003   7.420003   7.420003   3.120005   3.120005   3.120005
[61]   1.640006   1.640006   1.640006   1.640006   1.640006   1.640006
> (perc.se <- (ub-forecast$pred)/1.96/forecast$pred)
Time Series:
Start = 52 
End = 66 
Frequency = 1 
 [1] 0.2746454 0.3884072 0.4756998 0.8231958 0.9203609 1.0082048 2.5898286
 [8] 2.7686432 2.9365896 5.8888666 6.1684922 6.4359802 6.4359803 6.4359804
[15] 6.4359805
> postscript(file="/var/www/html/freestat/rcomp/tmp/1o9u71293200127.ps",horizontal=F,onefile=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/freestat/rcomp/tmp/2daaj1293200127.ps",horizontal=F,onefile=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/freestat/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> 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/3kc7v1293200127.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/freestat/rcomp/tmp/45u5j1293200127.tab") 
> 
> try(system("convert tmp/1o9u71293200127.ps tmp/1o9u71293200127.png",intern=TRUE))
character(0)
> try(system("convert tmp/2daaj1293200127.ps tmp/2daaj1293200127.png",intern=TRUE))
character(0)
> 
> 
> proc.time()
   user  system elapsed 
  0.954   0.484   1.037