R version 2.6.0 (2007-10-03)
Copyright (C) 2007 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.

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(7272.2,6680.1,8427.6,8752.8,7952.7,8694.3,7787,8474.2,9154.7,8557.2,7951.1,9156.7,7865.7,7337.4,9131.7,8814.6,8598.8,8439.6,7451.8,8016.2,9544.1,8270.7,8102.2,9369,7657.7,7816.6,9391.3,9445.4,9533.1,10068.7,8955.5,10423.9,11617.2,9391.1,10872,10230.4,9221,9428.6,10934.5,10986,11724.6,11180.9,11163.2,11240.9,12107.1,10762.3,11340.4,11266.8,9542.7,9227.7,10571.9,10774.4,10392.8,9920.2,9884.9,10174.5,11395.4,10760.2,10570.1,10536,9902.6,8889,10837.3,11624.1,10509,10984.9,10649.1,10855.7,11677.4,10760.2,10046.2,10772.8,9987.7,8638.7,11063.7,11855.7,10684.5,11337.4,10478,11123.9,12909.3,11339.9,10462.2,12733.5,10519.2,10414.9,12476.8,12384.6,12266.7,12919.9,11497.3,12142,13919.4,12656.8,12034.1,13199.7,10881.3,11301.2,13643.9,12517,13981.1,14275.7,13435,13565.7,16216.3,12970,14079.9,14235,12213.4,12581,14130.4,14210.8,14378.5,13142.8,13714.7,13621.9,15379.8,13306.3,14391.2,14909.9)
> par10 = 'FALSE'
> par9 = '1'
> par8 = '2'
> 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: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/
> #Source of accompanying publication: Office for Research, Development, and Education
> #Technical description: Write here your technical program description (don't use hard returns!)
> 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.3821  0.0714  0.3314  -0.4101  0.2143  -0.1860  -0.9999
s.e.   0.2795  0.2355  0.1351   0.2833  0.1447   0.1249   0.3798

sigma^2 estimated as 192830:  log likelihood = -726.52,  aic = 1469.03
> (forecast <- predict(arima.out,fx))
$pred
Time Series:
Start = 109 
End = 120 
Frequency = 1 
 [1] 12374.18 12670.19 14260.82 14150.18 14464.60 14480.34 13962.64 14381.57
 [9] 15906.44 14075.02 14427.12 14852.99

$se
Time Series:
Start = 109 
End = 120 
Frequency = 1 
 [1]  465.7956  475.6079  546.9383  639.0860  669.2594  737.8752  788.6777
 [8]  828.9809  881.0825  921.4760  962.0671 1003.2380

> (lb <- forecast$pred - 1.96 * forecast$se)
Time Series:
Start = 109 
End = 120 
Frequency = 1 
 [1] 11461.22 11737.99 13188.82 12897.57 13152.85 13034.10 12416.83 12756.77
 [9] 14179.52 12268.93 12541.47 12886.64
> (ub <- forecast$pred + 1.96 * forecast$se)
Time Series:
Start = 109 
End = 120 
Frequency = 1 
 [1] 13287.14 13602.38 15332.82 15402.79 15776.35 15926.57 15508.44 16006.37
 [9] 17633.36 15881.11 16312.77 16819.33
> 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]  7272.20  6680.10  8427.60  8752.80  7952.70  8694.30  7787.00  8474.20
  [9]  9154.70  8557.20  7951.10  9156.70  7865.70  7337.40  9131.70  8814.60
 [17]  8598.80  8439.60  7451.80  8016.20  9544.10  8270.70  8102.20  9369.00
 [25]  7657.70  7816.60  9391.30  9445.40  9533.10 10068.70  8955.50 10423.90
 [33] 11617.20  9391.10 10872.00 10230.40  9221.00  9428.60 10934.50 10986.00
 [41] 11724.60 11180.90 11163.20 11240.90 12107.10 10762.30 11340.40 11266.80
 [49]  9542.70  9227.70 10571.90 10774.40 10392.80  9920.20  9884.90 10174.50
 [57] 11395.40 10760.20 10570.10 10536.00  9902.60  8889.00 10837.30 11624.10
 [65] 10509.00 10984.90 10649.10 10855.70 11677.40 10760.20 10046.20 10772.80
 [73]  9987.70  8638.70 11063.70 11855.70 10684.50 11337.40 10478.00 11123.90
 [81] 12909.30 11339.90 10462.20 12733.50 10519.20 10414.90 12476.80 12384.60
 [89] 12266.70 12919.90 11497.30 12142.00 13919.40 12656.80 12034.10 13199.70
 [97] 10881.30 11301.20 13643.90 12517.00 13981.10 14275.70 13435.00 13565.70
[105] 16216.30 12970.00 14079.90 14235.00 12374.18 12670.19 14260.82 14150.18
[113] 14464.60 14480.34 13962.64 14381.57 15906.44 14075.02 14427.12 14852.99
> (perc.se <- (ub-forecast$pred)/1.96/forecast$pred)
Time Series:
Start = 109 
End = 120 
Frequency = 1 
 [1] 0.03764254 0.03753756 0.03835252 0.04516451 0.04626878 0.05095704
 [7] 0.05648487 0.05764188 0.05539157 0.06546889 0.06668461 0.06754453
> postscript(file="/var/www/html/rcomp/tmp/1pg4w1199354743.ps",horizontal=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.se <- array(0, dim=fx)
> perf.mse <- 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.mape[i] = perf.mape[i] + abs(perf.pe[i])
+ perf.se[i] = (x[nx+i] - forecast$pred[i])^2
+ perf.mse[i] = perf.mse[i] + perf.se[i]
+ 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 = perf.mape / fx
> perf.mse = perf.mse / fx
> perf.rmse = sqrt(perf.mse)
> postscript(file="/var/www/html/rcomp/tmp/2ypc71199354744.ps",horizontal=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:12] <- x[(nx+1):lx]
> lines(dum, lty=1)
> lines(ub,lty=3)
> lines(lb,lty=3)
> dev.off()
null device 
          1 
> load(file='/var/www/html/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/rcomp/tmp/3o78t1199354744.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.mape[i],4))
+ a<-table.element(a,round(perf.se[i],4))
+ a<-table.element(a,round(perf.mse[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/rcomp/tmp/4l9lc1199354744.tab") 
> 
> system("convert tmp/1pg4w1199354743.ps tmp/1pg4w1199354743.png")
> system("convert tmp/2ypc71199354744.ps tmp/2ypc71199354744.png")
> 
> 
> proc.time()
   user  system elapsed 
  4.550   0.512  13.374