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(364,351,380,319,322,386,221,187,343,342,365,313,356,337,389,326,343,357,220,218,391,425,332,298,360,336,325,393,301,426,265,210,429,440,357,431,442,422,544,420,396,482,261,211,448,468,464,425,415,433,531,457,380,481,302,216,509,417,370)
> par20 = ''
> par19 = ''
> par18 = ''
> par17 = ''
> par16 = ''
> par15 = ''
> par14 = ''
> par13 = ''
> par12 = ''
> par11 = ''
> par10 = 'TRUE'
> par9 = '0'
> par8 = '2'
> par7 = '0'
> par6 = '3'
> par5 = '6'
> par4 = '1'
> par3 = '0'
> par2 = '-0.7'
> par1 = '24'
> ylab = ''
> xlab = ''
> main = ''
> #'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#output/
> #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     sar1     sar2
      -0.0053  0.0347  0.0731  -0.9826  -0.0311
s.e.   0.1920  0.1915  0.2063   0.2230   0.2315

sigma^2 estimated as 1.801e-06:  log likelihood = 143.49,  aic = -274.97
> (forecast <- predict(arima.out,par1))
$pred
Time Series:
Start = 36 
End = 59 
Frequency = 1 
 [1] 0.01831307 0.01645323 0.01728132 0.01731791 0.01522068 0.01832879
 [7] 0.01462927 0.01993916 0.02336506 0.01451138 0.01416720 0.01643465
[13] 0.01812825 0.01662816 0.01758646 0.01717713 0.01516783 0.01823378
[19] 0.01480482 0.01977304 0.02307514 0.01464512 0.01421740 0.01652492

$se
Time Series:
Start = 36 
End = 59 
Frequency = 1 
 [1] 0.001341988 0.001342007 0.001342817 0.001346357 0.001346357 0.001346374
 [7] 0.001346718 0.001346718 0.001346719 0.001346721 0.001346721 0.001346721
[13] 0.001856205 0.001856217 0.001856748 0.001859069 0.001859069 0.001859081
[19] 0.001861404 0.001861404 0.001861408 0.001861420 0.001861420 0.001861420

> (lb <- forecast$pred - 1.96 * forecast$se)
Time Series:
Start = 36 
End = 59 
Frequency = 1 
 [1] 0.01568278 0.01382289 0.01464940 0.01467905 0.01258182 0.01568989
 [7] 0.01198970 0.01729960 0.02072549 0.01187180 0.01152763 0.01379508
[13] 0.01449009 0.01298997 0.01394723 0.01353335 0.01152405 0.01458998
[19] 0.01115647 0.01612469 0.01942678 0.01099674 0.01056902 0.01287653
> (ub <- forecast$pred + 1.96 * forecast$se)
Time Series:
Start = 36 
End = 59 
Frequency = 1 
 [1] 0.02094337 0.01908356 0.01991324 0.01995676 0.01785954 0.02096768
 [7] 0.01726884 0.02257873 0.02600463 0.01715095 0.01680677 0.01907423
[13] 0.02176641 0.02026634 0.02122569 0.02082090 0.01881161 0.02187758
[19] 0.01845317 0.02342140 0.02672350 0.01829350 0.01786579 0.02017330
> 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] 364.0000 351.0000 380.0000 319.0000 322.0000 386.0000 221.0000 187.0000
 [9] 343.0000 342.0000 365.0000 313.0000 356.0000 337.0000 389.0000 326.0000
[17] 343.0000 357.0000 220.0000 218.0000 391.0000 425.0000 332.0000 298.0000
[25] 360.0000 336.0000 325.0000 393.0000 301.0000 426.0000 265.0000 210.0000
[33] 429.0000 440.0000 357.0000 303.2283 353.3562 329.4179 328.4242 394.9325
[41] 302.8570 417.9367 268.5283 214.1010 422.7957 437.5451 353.9269 307.6544
[49] 348.0577 321.2832 332.2762 396.8997 305.1137 410.8751 271.7570 217.9542
[57] 417.2906 435.3396 351.1683
> (perc.se <- (ub-forecast$pred)/1.96/forecast$pred)
Time Series:
Start = 36 
End = 59 
Frequency = 1 
 [1] 0.12650566 0.14415842 0.13582805 0.13591386 0.15948516 0.12687314
 [7] 0.16773842 0.11474600 0.09530174 0.16947483 0.17475414 0.14498589
[13] 0.19241881 0.21579835 0.20033392 0.20703769 0.24523325 0.19134971
[19] 0.25412724 0.17259008 0.14220553 0.25804100 0.26912935 0.21844064
> postscript(file="/var/www/html/freestat/rcomp/tmp/1t7gl1293446602.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/20qwx1293446602.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/3p9tr1293446602.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/4ar9e1293446602.tab") 
> 
> try(system("convert tmp/1t7gl1293446602.ps tmp/1t7gl1293446602.png",intern=TRUE))
character(0)
> try(system("convert tmp/20qwx1293446602.ps tmp/20qwx1293446602.png",intern=TRUE))
character(0)
> 
> 
> proc.time()
   user  system elapsed 
  1.123   0.489   1.221