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(296.95,296.84,287.54,287.81,283.99,275.79,269.52,278.35,283.43,289.46,282.30,293.55,304.78,300.99,315.29,316.21,331.79,329.38,317.27,317.98,340.28,339.21,336.71,340.11,347.72,328.68,303.05,299.83,320.04,317.94,303.31,308.85,319.19,314.52,312.39,315.77,320.23,309.45,296.54,297.28,301.39,306.68,305.91,314.76,323.34,341.58,330.12,318.16,317.84,325.39,327.56,329.77,333.29,346.10,358.00,344.82,313.30,301.26,306.38,319.31)
> par10 = 'FALSE'
> par9 = '0'
> par8 = '2'
> par7 = '1'
> par6 = '0'
> par5 = '4'
> par4 = '0'
> 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:
         ma1    sar1    sar2
      0.2008  0.2593  0.3645
s.e.  0.1654  0.1427  0.1570

sigma^2 estimated as 69.93:  log likelihood = -124.95,  aic = 257.9
> (forecast <- predict(arima.out,par1))
$pred
Time Series:
Start = 37 
End = 60 
Frequency = 1 
 [1] 326.0216 324.0450 318.1599 321.0558 327.4836 325.2687 322.9660 324.9491
 [9] 330.3529 329.0580 326.3156 327.8855 331.6300 330.4868 328.9362 330.0662
[17] 333.0070 332.2386 330.8368 331.7021 333.8297 333.2137 332.2849 332.9212

$se
Time Series:
Start = 37 
End = 60 
Frequency = 1 
 [1]  8.362584 13.067666 16.480139 19.298344 22.836754 26.104240 29.005953
 [8] 31.642682 35.574290 39.419470 42.921546 46.158678 49.805856 53.326503
[15] 56.628690 59.748650 63.271891 66.720570 69.999548 73.131655 76.444443
[22] 79.680799 82.790741 85.788016

> (lb <- forecast$pred - 1.96 * forecast$se)
Time Series:
Start = 37 
End = 60 
Frequency = 1 
 [1] 309.6310 298.4324 285.8588 283.2311 282.7235 274.1044 266.1143 262.9294
 [9] 260.6273 251.7958 242.1894 237.4145 234.0105 225.9668 217.9440 212.9588
[17] 208.9941 201.4662 193.6377 188.3641 183.9986 177.0393 170.0151 164.7767
> (ub <- forecast$pred + 1.96 * forecast$se)
Time Series:
Start = 37 
End = 60 
Frequency = 1 
 [1] 342.4123 349.6576 350.4609 358.8806 372.2436 376.4330 379.8177 386.9687
 [9] 400.0785 406.3201 410.4418 418.3565 429.2494 435.0067 439.9284 447.1736
[17] 457.0199 463.0109 468.0359 475.0401 483.6608 489.3880 494.5548 501.0658
> 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] 296.9500 296.8400 287.5400 287.8100 283.9900 275.7900 269.5200 278.3500
 [9] 283.4300 289.4600 282.3000 293.5500 304.7800 300.9900 315.2900 316.2100
[17] 331.7900 329.3800 317.2700 317.9800 340.2800 339.2100 336.7100 340.1100
[25] 347.7200 328.6800 303.0500 299.8300 320.0400 317.9400 303.3100 308.8500
[33] 319.1900 314.5200 312.3900 315.7700 326.0216 324.0450 318.1599 321.0558
[41] 327.4836 325.2687 322.9660 324.9491 330.3529 329.0580 326.3156 327.8855
[49] 331.6300 330.4868 328.9362 330.0662 333.0070 332.2386 330.8368 331.7021
[57] 333.8297 333.2137 332.2849 332.9212
> (perc.se <- (ub-forecast$pred)/1.96/forecast$pred)
Time Series:
Start = 37 
End = 60 
Frequency = 1 
 [1] 0.02565040 0.04032670 0.05179830 0.06010900 0.06973404 0.08025439
 [7] 0.08981117 0.09737735 0.10768572 0.11979490 0.13153384 0.14077681
[13] 0.15018504 0.16135745 0.17215707 0.18102020 0.19000166 0.20082128
[19] 0.21158332 0.22047390 0.22899235 0.23912824 0.24915587 0.25768261
> postscript(file="/var/www/html/freestat/rcomp/tmp/1bd4x1291995839.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/27nkp1291995839.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/3w6hi1291995839.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/4hof61291995839.tab") 
> 
> try(system("convert tmp/1bd4x1291995839.ps tmp/1bd4x1291995839.png",intern=TRUE))
character(0)
> try(system("convert tmp/27nkp1291995839.ps tmp/27nkp1291995839.png",intern=TRUE))
character(0)
> 
> 
> proc.time()
   user  system elapsed 
  0.946   0.431   1.015