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(2617.2,2506.13,2679.07,2589.73,2457.46,2517.3,2386.53,2453.37,2529.66,2475.14,2525.93,2480.93,2229.85,2169.14,2030.98,2071.37,1953.35,1748.74,1696.58,1900.09,1908.64,1881.46,2100.18,2672.2,3136,2994.38,3168.22,3751.41,3925.43,3719.52,3757.12,3722.23,4127.47,4162.5,4441.82,4325.29,4350.83,4384.47,4639.4,4697.86,4614.76,4471.65,4305.23,4433.57,4388.53,4140.3,4144.38,4070.78,3906.01,3795.91,3703.32,3675.8,3911.06,3912.28,3839.25,3744.63,3549.25,3394.14,3264.26,3328.8,3223.98,3228.01,3112.83,3051.67,3039.71,3125.67,3106.54)
> par10 = 'FALSE'
> par9 = '1'
> par8 = '2'
> par7 = '1'
> par6 = '3'
> par5 = '12'
> par4 = '1'
> par3 = '1'
> par2 = '1'
> par1 = '24'
> #'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     sar1     sar2     sma1
      0.8969  -0.2362  0.2125  -0.6133  -0.7280  -0.8484  -0.9785
s.e.  0.3890   0.2575  0.2067   0.3285   0.1797   0.1118   2.0333

sigma^2 estimated as 6715:  log likelihood = -204.47,  aic = 424.95
> (forecast <- predict(arima.out,par1))
$pred
Time Series:
Start = 44 
End = 67 
Frequency = 1 
 [1] 4350.980 4463.147 4382.045 4430.307 4317.829 4035.282 3918.616 3806.685
 [9] 3753.988 3601.255 3496.125 3423.484 3600.945 3591.810 3551.508 3721.705
[17] 4193.289 4510.104 4357.535 4446.723 4922.561 5040.895 4856.883 4890.878

$se
Time Series:
Start = 44 
End = 67 
Frequency = 1 
 [1] 110.8613 179.2101 227.8022 276.1855 327.5276 377.5520 419.9684 462.5282
 [9] 504.9288 546.7648 587.7735 627.8103 657.1265 683.6512 709.6485 733.9016
[17] 756.3769 777.0412 794.6930 811.3498 827.0786 841.9746 856.1030 869.5115

> (lb <- forecast$pred - 1.96 * forecast$se)
Time Series:
Start = 44 
End = 67 
Frequency = 1 
 [1] 4133.691 4111.895 3935.553 3888.984 3675.875 3295.280 3095.478 2900.130
 [9] 2764.327 2529.596 2344.089 2192.976 2312.977 2251.854 2160.597 2283.258
[17] 2710.790 2987.103 2799.937 2856.478 3301.487 3390.625 3178.921 3186.635
> (ub <- forecast$pred + 1.96 * forecast$se)
Time Series:
Start = 44 
End = 67 
Frequency = 1 
 [1] 4568.268 4814.398 4828.537 4971.631 4959.783 4775.284 4741.754 4713.240
 [9] 4743.648 4672.914 4648.161 4653.992 4888.912 4931.766 4942.419 5160.152
[17] 5675.788 6033.104 5915.134 6036.969 6543.635 6691.165 6534.844 6595.120
> 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] 2617.200 2506.130 2679.070 2589.730 2457.460 2517.300 2386.530 2453.370
 [9] 2529.660 2475.140 2525.930 2480.930 2229.850 2169.140 2030.980 2071.370
[17] 1953.350 1748.740 1696.580 1900.090 1908.640 1881.460 2100.180 2672.200
[25] 3136.000 2994.380 3168.220 3751.410 3925.430 3719.520 3757.120 3722.230
[33] 4127.470 4162.500 4441.820 4325.290 4350.830 4384.470 4639.400 4697.860
[41] 4614.760 4471.650 4305.230 4350.980 4463.147 4382.045 4430.307 4317.829
[49] 4035.282 3918.616 3806.685 3753.988 3601.255 3496.125 3423.484 3600.945
[57] 3591.810 3551.508 3721.705 4193.289 4510.104 4357.535 4446.723 4922.561
[65] 5040.895 4856.883 4890.878
> (perc.se <- (ub-forecast$pred)/1.96/forecast$pred)
Time Series:
Start = 44 
End = 67 
Frequency = 1 
 [1] 0.02547961 0.04015331 0.05198537 0.06234004 0.07585470 0.09356274
 [7] 0.10717262 0.12150419 0.13450466 0.15182618 0.16812140 0.18338343
[13] 0.18248726 0.19033611 0.19981612 0.19719498 0.18037798 0.17228898
[19] 0.18237212 0.18246015 0.16801795 0.16702878 0.17626594 0.17778230
> postscript(file="/var/www/html/freestat/rcomp/tmp/1z09q1292621800.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/2vspz1292621800.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/32bmb1292621800.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/45b3h1292621800.tab") 
> 
> try(system("convert tmp/1z09q1292621800.ps tmp/1z09q1292621800.png",intern=TRUE))
character(0)
> try(system("convert tmp/2vspz1292621800.ps tmp/2vspz1292621800.png",intern=TRUE))
character(0)
> 
> 
> proc.time()
   user  system elapsed 
  4.678   0.832   5.200