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(12008,9169,8788,8417,8247,8197,8236,8253,7733,8366,8626,8863,10102,8463,9114,8563,8872,8301,8301,8278,7736,7973,8268,9476,11100,8962,9173,8738,8459,8078,8411,8291,7810,8616,8312,9692,9911,8915,9452,9112,8472,8230,8384,8625,8221,8649,8625,10443,10357,8586,8892,8329,8101,7922,8120,7838,7735,8406,8209,9451,10041,9411,10405,8467,8464,8102,7627,7513,7510,8291,8064,9383,9706,8579,9474,8318,8213,8059,9111,7708,7680,8014,8007,8718,9486,9113,9025,8476,7952,7759,7835,7600,7651,8319,8812,8630)
> par10 = 'FALSE'
> par9 = '1'
> par8 = '2'
> par7 = '1'
> par6 = '1'
> 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      ma1     sar1    sar2     sma1
      -0.0427  -0.3029  -0.2407  0.3319  -0.2672
s.e.   0.0058   0.1695      NaN     NaN   0.1835

sigma^2 estimated as 193827:  log likelihood = -446.14,  aic = 904.29
Warning message:
In sqrt(diag(x$var.coef)) : NaNs produced
> (forecast <- predict(arima.out,par1))
$pred
Time Series:
Start = 73 
End = 96 
Frequency = 1 
 [1] 9615.146 8305.422 8919.171 7641.888 7579.740 7331.305 7201.987 6917.471
 [9] 6927.564 7687.286 7447.901 8643.451 9155.062 8389.022 9322.617 7429.887
[17] 7456.659 7120.143 6684.223 6496.521 6536.653 7338.008 7091.648 8342.475

$se
Time Series:
Start = 73 
End = 96 
Frequency = 1 
 [1]  440.2600  526.1553  603.0207  671.0181  732.7371  789.6466  842.7216
 [8]  892.6465  939.9233  984.9334 1027.9747 1069.2848 1185.1013 1262.7971
[15] 1337.0285 1407.3074 1474.2415 1538.2658 1599.7298 1658.9181 1716.0661
[22] 1771.3714 1825.0015 1877.1000

> (lb <- forecast$pred - 1.96 * forecast$se)
Time Series:
Start = 73 
End = 96 
Frequency = 1 
 [1] 8752.236 7274.158 7737.250 6326.693 6143.575 5783.598 5550.253 5167.884
 [9] 5085.314 5756.816 5433.071 6547.652 6832.264 5913.940 6702.041 4671.564
[17] 4567.145 4105.142 3548.753 3245.042 3173.163 3866.120 3514.645 4663.359
> (ub <- forecast$pred + 1.96 * forecast$se)
Time Series:
Start = 73 
End = 96 
Frequency = 1 
 [1] 10478.056  9336.687 10101.092  8957.084  9015.904  8879.012  8853.721
 [8]  8667.059  8769.813  9617.755  9462.732 10739.249 11477.861 10864.105
[15] 11943.193 10188.209 10346.172 10135.144  9819.694  9748.001  9900.142
[22] 10809.896 10668.651 12021.591
> 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] 12008.000  9169.000  8788.000  8417.000  8247.000  8197.000  8236.000
 [8]  8253.000  7733.000  8366.000  8626.000  8863.000 10102.000  8463.000
[15]  9114.000  8563.000  8872.000  8301.000  8301.000  8278.000  7736.000
[22]  7973.000  8268.000  9476.000 11100.000  8962.000  9173.000  8738.000
[29]  8459.000  8078.000  8411.000  8291.000  7810.000  8616.000  8312.000
[36]  9692.000  9911.000  8915.000  9452.000  9112.000  8472.000  8230.000
[43]  8384.000  8625.000  8221.000  8649.000  8625.000 10443.000 10357.000
[50]  8586.000  8892.000  8329.000  8101.000  7922.000  8120.000  7838.000
[57]  7735.000  8406.000  8209.000  9451.000 10041.000  9411.000 10405.000
[64]  8467.000  8464.000  8102.000  7627.000  7513.000  7510.000  8291.000
[71]  8064.000  9383.000  9615.146  8305.422  8919.171  7641.888  7579.740
[78]  7331.305  7201.987  6917.471  6927.564  7687.286  7447.901  8643.451
[85]  9155.062  8389.022  9322.617  7429.887  7456.659  7120.143  6684.223
[92]  6496.521  6536.653  7338.008  7091.648  8342.475
> (perc.se <- (ub-forecast$pred)/1.96/forecast$pred)
Time Series:
Start = 73 
End = 96 
Frequency = 1 
 [1] 0.04578817 0.06335082 0.06760950 0.08780790 0.09667049 0.10770887
 [7] 0.11701238 0.12904231 0.13567877 0.12812499 0.13802206 0.12371041
[13] 0.12944764 0.15052971 0.14341772 0.18941169 0.19770806 0.21604424
[19] 0.23932920 0.25535484 0.26252980 0.24139676 0.25734519 0.22500516
> postscript(file="/var/www/html/freestat/rcomp/tmp/1bdmx1291738765.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/2zejr1291738765.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/36fyk1291738765.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/4sfwq1291738765.tab") 
> 
> try(system("convert tmp/1bdmx1291738765.ps tmp/1bdmx1291738765.png",intern=TRUE))
character(0)
> try(system("convert tmp/2zejr1291738765.ps tmp/2zejr1291738765.png",intern=TRUE))
character(0)
> 
> 
> proc.time()
   user  system elapsed 
  2.448   0.510   2.595