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(6.4,7.7,9.2,8.6,7.4,8.6,6.2,6,6.6,5.1,4.7,5,3.6,1.9,-0.1,-5.7,-5.6,-6.4,-7.7,-8,-11.9,-15.4,-15.5,-13.4,-10.9,-10.8,-7.3,-6.5,-5.1,-5.3,-6.8,-8.4,-8.4,-9.7,-8.8,-9.6,-11.5,-11,-14.9,-16.2,-14.4,-17.3,-15.7,-12.6,-9.4,-8.1,-5.4,-4.6,-4.9,-4,-3.1,-1.3,0,-0.4,3,0.4,1.2,0.6,-1.3,-3.2,-1.8,-3.6,-4.2,-6.9,-8,-7.5,-8.2,-7.6,-3.7,-1.7,-0.7,0.2,0.6,2.2,3.3,5.3,5.5,6.3,7.7,6.5,5.5,6.9,5.7,6.9,6.1,4.8,3.7,5.8,6.8,8.5,7.2,5,4.7,2.3,2.4,0.1,1.9,1.7,2,-1.9,0.5,-1.3,-3.3,-2.8,-8,-13.9,-21.9,-28.8,-27.6,-31.4,-31.8,-29.4,-27.6,-23.6,-22.8,-18.2,-17.8,-14.2,-8.8,-7.9,-7,-7,-3.6,-2.4,-4.9,-7.7,-6.5,-5.1,-3.4,-2.8,0.8)
> par10 = 'FALSE'
> par9 = '0'
> par8 = '0'
> par7 = '1'
> 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: 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:
          ar1     ar2     ar3     ma1
      -0.2797  0.3919  0.3059  0.5080
s.e.   0.2574  0.1121  0.1028  0.2565

sigma^2 estimated as 3.987:  log likelihood = -249.26,  aic = 508.52
> (forecast <- predict(arima.out,par1))
$pred
Time Series:
Start = 120 
End = 131 
Frequency = 1 
 [1] -6.9345399 -4.2390219 -2.6102290 -1.4388959 -0.3037490  0.3359694
 [7]  0.9601571  1.3834711  1.7053471  1.9721266  2.1531239  2.3054966

$se
Time Series:
Start = 120 
End = 131 
Frequency = 1 
 [1]  1.996701  3.162573  4.433800  5.783614  6.998272  8.221501  9.364894
 [8] 10.457082 11.500550 12.488073 13.432841 14.332641

> (lb <- forecast$pred - 1.96 * forecast$se)
Time Series:
Start = 120 
End = 131 
Frequency = 1 
 [1] -10.84807 -10.43767 -11.30048 -12.77478 -14.02036 -15.77817 -17.39504
 [8] -19.11241 -20.83573 -22.50450 -24.17525 -25.78648
> (ub <- forecast$pred + 1.96 * forecast$se)
Time Series:
Start = 120 
End = 131 
Frequency = 1 
 [1] -3.021005  1.959622  6.080020  9.896987 13.412863 16.450111 19.315349
 [8] 21.879351 24.246426 26.448749 28.481493 30.397472
> 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]   6.4000000   7.7000000   9.2000000   8.6000000   7.4000000   8.6000000
  [7]   6.2000000   6.0000000   6.6000000   5.1000000   4.7000000   5.0000000
 [13]   3.6000000   1.9000000  -0.1000000  -5.7000000  -5.6000000  -6.4000000
 [19]  -7.7000000  -8.0000000 -11.9000000 -15.4000000 -15.5000000 -13.4000000
 [25] -10.9000000 -10.8000000  -7.3000000  -6.5000000  -5.1000000  -5.3000000
 [31]  -6.8000000  -8.4000000  -8.4000000  -9.7000000  -8.8000000  -9.6000000
 [37] -11.5000000 -11.0000000 -14.9000000 -16.2000000 -14.4000000 -17.3000000
 [43] -15.7000000 -12.6000000  -9.4000000  -8.1000000  -5.4000000  -4.6000000
 [49]  -4.9000000  -4.0000000  -3.1000000  -1.3000000   0.0000000  -0.4000000
 [55]   3.0000000   0.4000000   1.2000000   0.6000000  -1.3000000  -3.2000000
 [61]  -1.8000000  -3.6000000  -4.2000000  -6.9000000  -8.0000000  -7.5000000
 [67]  -8.2000000  -7.6000000  -3.7000000  -1.7000000  -0.7000000   0.2000000
 [73]   0.6000000   2.2000000   3.3000000   5.3000000   5.5000000   6.3000000
 [79]   7.7000000   6.5000000   5.5000000   6.9000000   5.7000000   6.9000000
 [85]   6.1000000   4.8000000   3.7000000   5.8000000   6.8000000   8.5000000
 [91]   7.2000000   5.0000000   4.7000000   2.3000000   2.4000000   0.1000000
 [97]   1.9000000   1.7000000   2.0000000  -1.9000000   0.5000000  -1.3000000
[103]  -3.3000000  -2.8000000  -8.0000000 -13.9000000 -21.9000000 -28.8000000
[109] -27.6000000 -31.4000000 -31.8000000 -29.4000000 -27.6000000 -23.6000000
[115] -22.8000000 -18.2000000 -17.8000000 -14.2000000  -8.8000000  -6.9345399
[121]  -4.2390219  -2.6102290  -1.4388959  -0.3037490   0.3359694   0.9601571
[127]   1.3834711   1.7053471   1.9721266   2.1531239   2.3054966
> (perc.se <- (ub-forecast$pred)/1.96/forecast$pred)
Time Series:
Start = 120 
End = 131 
Frequency = 1 
 [1]  -0.2879357  -0.7460621  -1.6986251  -4.0194802 -23.0396553  24.4709802
 [7]   9.7535023   7.5585834   6.7438180   6.3322874   6.2387685   6.2167260
> postscript(file="/var/www/html/freestat/rcomp/tmp/1cmi81293202081.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/2jnxk1293202081.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/376bw1293202081.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/4gsz81293202081.tab") 
> 
> try(system("convert tmp/1cmi81293202081.ps tmp/1cmi81293202081.png",intern=TRUE))
character(0)
> try(system("convert tmp/2jnxk1293202081.ps tmp/2jnxk1293202081.png",intern=TRUE))
character(0)
> 
> 
> proc.time()
   user  system elapsed 
  0.962   0.474   1.039