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(9769,9321,9939,9336,10195,9464,10010,10213,9563,9890,9305,9391,9928,8686,9843,9627,10074,9503,10119,10000,9313,9866,9172,9241,9659,8904,9755,9080,9435,8971,10063,9793,9454,9759,8820,9403,9676,8642,9402,9610,9294,9448,10319,9548,9801,9596,8923,9746,9829,9125,9782,9441,9162,9915,10444,10209,9985,9842,9429,10132,9849,9172,10313,9819,9955,10048,10082,10541,10208,10233,9439,9963,10158,9225,10474,9757,10490,10281,10444,10640,10695,10786,9832,9747,10411,9511,10402,9701,10540,10112,10915,11183,10384,10834,9886,10216)
> par10 = 'FALSE'
> par9 = '1'
> par8 = '1'
> par7 = '1'
> par6 = '1'
> par5 = '12'
> par4 = '1'
> par3 = '0'
> par2 = '-0.8'
> 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      ma1    sar1     sma1
      0.9937  -0.8286  0.4595  -1.0000
s.e.  0.0253   0.0555  0.1503   0.1611

sigma^2 estimated as 1.815e-10:  log likelihood = 698.6,  aic = -1387.19
Warning messages:
1: In log(s2) : NaNs produced
2: In log(s2) : NaNs produced
> (forecast <- predict(arima.out,par1))
$pred
Time Series:
Start = 85 
End = 96 
Frequency = 1 
 [1] 0.0006139382 0.0006613806 0.0006042810 0.0006336525 0.0006067729
 [6] 0.0006171626 0.0005993466 0.0005962363 0.0006031321 0.0005971788
[11] 0.0006394027 0.0006322278

$se
Time Series:
Start = 85 
End = 96 
Frequency = 1 
 [1] 1.418817e-05 1.436303e-05 1.453150e-05 1.469383e-05 1.485024e-05
 [6] 1.500094e-05 1.514613e-05 1.528598e-05 1.542067e-05 1.555034e-05
[11] 1.567515e-05 1.579522e-05

> (lb <- forecast$pred - 1.96 * forecast$se)
Time Series:
Start = 85 
End = 96 
Frequency = 1 
 [1] 0.0005861294 0.0006332291 0.0005757993 0.0006048526 0.0005776664
 [6] 0.0005877607 0.0005696602 0.0005662757 0.0005729076 0.0005667001
[11] 0.0006086795 0.0006012692
> (ub <- forecast$pred + 1.96 * forecast$se)
Time Series:
Start = 85 
End = 96 
Frequency = 1 
 [1] 0.0006417470 0.0006895321 0.0006327628 0.0006624525 0.0006358794
 [6] 0.0006465644 0.0006290330 0.0006261968 0.0006333566 0.0006276574
[11] 0.0006701260 0.0006631864
> 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]  9769.000  9321.000  9939.000  9336.000 10195.000  9464.000 10010.000
 [8] 10213.000  9563.000  9890.000  9305.000  9391.000  9928.000  8686.000
[15]  9843.000  9627.000 10074.000  9503.000 10119.000 10000.000  9313.000
[22]  9866.000  9172.000  9241.000  9659.000  8904.000  9755.000  9080.000
[29]  9435.000  8971.000 10063.000  9793.000  9454.000  9759.000  8820.000
[36]  9403.000  9676.000  8642.000  9402.000  9610.000  9294.000  9448.000
[43] 10319.000  9548.000  9801.000  9596.000  8923.000  9746.000  9829.000
[50]  9125.000  9782.000  9441.000  9162.000  9915.000 10444.000 10209.000
[57]  9985.000  9842.000  9429.000 10132.000  9849.000  9172.000 10313.000
[64]  9819.000  9955.000 10048.000 10082.000 10541.000 10208.000 10233.000
[71]  9439.000  9963.000 10158.000  9225.000 10474.000  9757.000 10490.000
[78] 10281.000 10444.000 10640.000 10695.000 10786.000  9832.000  9747.000
[85] 10347.707  9428.350 10554.832  9946.860 10500.676 10280.175 10663.565
[92] 10733.146 10579.970 10711.975  9835.170  9974.887
> (perc.se <- (ub-forecast$pred)/1.96/forecast$pred)
Time Series:
Start = 85 
End = 96 
Frequency = 1 
 [1] 0.03043554 0.02850859 0.03173907 0.03054518 0.03233401 0.03209978
 [7] 0.03344875 0.03396263 0.03386464 0.03452776 0.03239149 0.03304593
> postscript(file="/var/www/html/freestat/rcomp/tmp/1yurg1291388311.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/2umpp1291388311.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/305401291388311.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/44o2o1291388311.tab") 
> 
> try(system("convert tmp/1yurg1291388311.ps tmp/1yurg1291388311.png",intern=TRUE))
character(0)
> try(system("convert tmp/2umpp1291388311.ps tmp/2umpp1291388311.png",intern=TRUE))
character(0)
> 
> 
> proc.time()
   user  system elapsed 
  1.882   0.529   2.057