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(5140,4749,3635,4305,5805,4260,3869,7325,9280,6222,3272,7598,1345,1900,1480,1472,3823,4454,3357,5393,8329,4152,4042,7747,1451,911,-406,1387,2150,1577,2642,4273,8064,3243,1112,2280,505,744,-1369,-531,1041,2076,577,5080,6584,3761,294,5020,1141,3805,2127,2531,3682,3263,2798,5936,10568,5296,1870,4390,3707,5201,3748,5282,5349,6249,5517,8640,15767,8850,5582,6496,3255,6189,6452,5099,6833,7046,7739,10142,16054,7721,6182,6490,3704,6235,4655,5072,3640,5147,5703,11889,15603,9589)
> par10 = 'FALSE'
> par9 = '1'
> par8 = '0'
> par7 = '1'
> par6 = '3'
> par5 = '12'
> par4 = '1'
> par3 = '1'
> 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     sma1
      -0.0848  -0.0440  -0.2759  -0.4372  -0.4250
s.e.   0.2494   0.1608   0.1296   0.2457   0.1233

sigma^2 estimated as 1720267:  log likelihood = -594.74,  aic = 1201.48
> (forecast <- predict(arima.out,par1))
$pred
Time Series:
Start = 83 
End = 94 
Frequency = 1 
 [1]  5790.706  7524.458  5437.615  7367.770  6889.782  6446.630  7814.106
 [8]  8134.622  8287.277 10986.032 16786.319  9479.559

$se
Time Series:
Start = 83 
End = 94 
Frequency = 1 
 [1] 1311.608 1453.773 1583.374 1610.729 1690.321 1765.197 1861.505 1933.937
 [9] 2004.571 2067.788 2133.973 2197.766

> (lb <- forecast$pred - 1.96 * forecast$se)
Time Series:
Start = 83 
End = 94 
Frequency = 1 
 [1]  3219.955  4675.062  2334.202  4210.742  3576.752  2986.843  4165.557
 [8]  4344.104  4358.318  6933.167 12603.731  5171.938
> (ub <- forecast$pred + 1.96 * forecast$se)
Time Series:
Start = 83 
End = 94 
Frequency = 1 
 [1]  8361.457 10373.854  8541.028 10524.798 10202.812  9906.417 11462.655
 [8] 11925.139 12216.237 15038.897 20968.906 13787.180
> 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]  5140.000  4749.000  3635.000  4305.000  5805.000  4260.000  3869.000
 [8]  7325.000  9280.000  6222.000  3272.000  7598.000  1345.000  1900.000
[15]  1480.000  1472.000  3823.000  4454.000  3357.000  5393.000  8329.000
[22]  4152.000  4042.000  7747.000  1451.000   911.000  -406.000  1387.000
[29]  2150.000  1577.000  2642.000  4273.000  8064.000  3243.000  1112.000
[36]  2280.000   505.000   744.000 -1369.000  -531.000  1041.000  2076.000
[43]   577.000  5080.000  6584.000  3761.000   294.000  5020.000  1141.000
[50]  3805.000  2127.000  2531.000  3682.000  3263.000  2798.000  5936.000
[57] 10568.000  5296.000  1870.000  4390.000  3707.000  5201.000  3748.000
[64]  5282.000  5349.000  6249.000  5517.000  8640.000 15767.000  8850.000
[71]  5582.000  6496.000  3255.000  6189.000  6452.000  5099.000  6833.000
[78]  7046.000  7739.000 10142.000 16054.000  7721.000  5790.706  7524.458
[85]  5437.615  7367.770  6889.782  6446.630  7814.106  8134.622  8287.277
[92] 10986.032 16786.319  9479.559
> (perc.se <- (ub-forecast$pred)/1.96/forecast$pred)
Time Series:
Start = 83 
End = 94 
Frequency = 1 
 [1] 0.2265022 0.1932064 0.2911891 0.2186182 0.2453374 0.2738171 0.2382236
 [8] 0.2377415 0.2418854 0.1882197 0.1271257 0.2318426
> postscript(file="/var/www/html/freestat/rcomp/tmp/1echc1293361928.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/2dnzj1293361929.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/3j6wd1293361929.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/4n6dj1293361929.tab") 
> 
> try(system("convert tmp/1echc1293361928.ps tmp/1echc1293361928.png",intern=TRUE))
character(0)
> try(system("convert tmp/2dnzj1293361929.ps tmp/2dnzj1293361929.png",intern=TRUE))
character(0)
> 
> 
> proc.time()
   user  system elapsed 
  1.379   0.492   1.462