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(46,62,66,59,58,61,41,27,58,70,49,59,44,36,72,45,56,54,53,35,61,52,47,51,52,63,74,45,51,64,36,30,55,64,39,40,63,45,59,55,40,64,27,28,45,57,45,69,60,56,58,50,51,53,37,22,55,70,62,58,39,49,58,47,42,62,39,40,72,70,54,65)
> par10 = 'FALSE'
> par9 = '1'
> par8 = '2'
> par7 = '2'
> par6 = '3'
> par5 = '12'
> par4 = '1'
> par3 = '0'
> 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     ma2    sar1    sar2     sma1
      0.6445  -0.0002  -0.2142  -0.6380  0.1048  0.0237  0.0142  -1.0000
s.e.  0.8287   0.6459   0.1537   0.8449  0.6265  0.2341  0.2776   0.6878

sigma^2 estimated as 62.22:  log likelihood = -176.65,  aic = 371.3
> (forecast <- predict(arima.out,par1))
$pred
Time Series:
Start = 61 
End = 72 
Frequency = 1 
 [1] 54.32558 49.48768 62.93191 49.04644 50.52331 59.32354 39.18743 28.74906
 [9] 54.89650 62.73908 48.13557 55.57387

$se
Time Series:
Start = 61 
End = 72 
Frequency = 1 
 [1] 8.633671 8.633854 8.682969 8.768128 8.804330 8.833061 8.835212 8.834560
 [9] 8.835446 8.834574 8.833840 8.833860

> (lb <- forecast$pred - 1.96 * forecast$se)
Time Series:
Start = 61 
End = 72 
Frequency = 1 
 [1] 37.40359 32.56533 45.91329 31.86091 33.26682 42.01074 21.87042 11.43333
 [9] 37.57902 45.42332 30.82124 38.25951
> (ub <- forecast$pred + 1.96 * forecast$se)
Time Series:
Start = 61 
End = 72 
Frequency = 1 
 [1] 71.24758 66.41004 79.95053 66.23197 67.77979 76.63634 56.50445 46.06480
 [9] 72.21397 80.05485 65.44989 72.88824
> 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] 46.00000 62.00000 66.00000 59.00000 58.00000 61.00000 41.00000 27.00000
 [9] 58.00000 70.00000 49.00000 59.00000 44.00000 36.00000 72.00000 45.00000
[17] 56.00000 54.00000 53.00000 35.00000 61.00000 52.00000 47.00000 51.00000
[25] 52.00000 63.00000 74.00000 45.00000 51.00000 64.00000 36.00000 30.00000
[33] 55.00000 64.00000 39.00000 40.00000 63.00000 45.00000 59.00000 55.00000
[41] 40.00000 64.00000 27.00000 28.00000 45.00000 57.00000 45.00000 69.00000
[49] 60.00000 56.00000 58.00000 50.00000 51.00000 53.00000 37.00000 22.00000
[57] 55.00000 70.00000 62.00000 58.00000 54.32558 49.48768 62.93191 49.04644
[65] 50.52331 59.32354 39.18743 28.74906 54.89650 62.73908 48.13557 55.57387
> (perc.se <- (ub-forecast$pred)/1.96/forecast$pred)
Time Series:
Start = 61 
End = 72 
Frequency = 1 
 [1] 0.1589246 0.1744647 0.1379740 0.1787720 0.1742627 0.1488964 0.2254604
 [8] 0.3072990 0.1609473 0.1408145 0.1835200 0.1589571
> postscript(file="/var/www/html/freestat/rcomp/tmp/10z6a1291396864.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/2f9m11291396864.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/3ej0y1291396864.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/4z1hl1291396864.tab") 
> 
> try(system("convert tmp/10z6a1291396864.ps tmp/10z6a1291396864.png",intern=TRUE))
character(0)
> try(system("convert tmp/2f9m11291396864.ps tmp/2f9m11291396864.png",intern=TRUE))
character(0)
> 
> 
> proc.time()
   user  system elapsed 
  6.313   1.059   7.033