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(15136
+ ,16733
+ ,20016
+ ,17708
+ ,18019
+ ,19227
+ ,22893
+ ,23739
+ ,21133
+ ,22591
+ ,26786
+ ,29740
+ ,15028
+ ,17977
+ ,20008
+ ,21354
+ ,19498
+ ,22125
+ ,25817
+ ,28779
+ ,20960
+ ,22254
+ ,27392
+ ,29945
+ ,16933
+ ,17892
+ ,20533
+ ,23569
+ ,22417
+ ,22084
+ ,26580
+ ,27454
+ ,24081
+ ,23451
+ ,28991
+ ,31386
+ ,16896
+ ,20045
+ ,23471
+ ,21747
+ ,25621
+ ,23859
+ ,25500
+ ,30998
+ ,24475
+ ,23145
+ ,29701
+ ,34365
+ ,17556
+ ,22077
+ ,25702
+ ,22214
+ ,26886
+ ,23191
+ ,27831
+ ,35406
+ ,23195
+ ,25110
+ ,30009
+ ,36242
+ ,18450
+ ,21845
+ ,26488
+ ,22394
+ ,28057
+ ,25451
+ ,24872
+ ,33424
+ ,24052
+ ,28449
+ ,33533
+ ,37351
+ ,19969
+ ,21701
+ ,26249
+ ,24493
+ ,24603
+ ,26485
+ ,30723
+ ,34569
+ ,26689
+ ,26157
+ ,32064
+ ,38870
+ ,21337
+ ,19419
+ ,23166
+ ,28286
+ ,24570
+ ,24001
+ ,33151
+ ,24878
+ ,26804
+ ,28967
+ ,33311
+ ,40226
+ ,20504
+ ,23060
+ ,23562
+ ,27562
+ ,23940
+ ,24584
+ ,34303
+ ,25517
+ ,23494
+ ,29095
+ ,32903
+ ,34379
+ ,16991
+ ,21109
+ ,23740
+ ,25552
+ ,21752
+ ,20294
+ ,29009
+ ,25500
+ ,24166
+ ,26960
+ ,31222
+ ,38641
+ ,14672
+ ,17543
+ ,25453
+ ,32683
+ ,22449
+ ,22316
+ ,27595
+ ,25451
+ ,25421
+ ,25288
+ ,32568
+ ,35110
+ ,16052
+ ,22146
+ ,21198
+ ,19543
+ ,22084
+ ,23816
+ ,29961
+ ,26773
+ ,26635
+ ,26972
+ ,30207
+ ,38687
+ ,16974
+ ,21697
+ ,24179
+ ,23757
+ ,25013
+ ,24019
+ ,30345
+ ,24488
+ ,25156
+ ,25650
+ ,30923
+ ,37240
+ ,17466
+ ,19463
+ ,24352
+ ,26805
+ ,25236
+ ,24735
+ ,29356
+ ,31234
+ ,22724
+ ,28496
+ ,32857
+ ,37198
+ ,13652
+ ,22784
+ ,23565
+ ,26323
+ ,23779
+ ,27549
+ ,29660
+ ,23356)
> par10 = 'FALSE'
> par9 = '1'
> par8 = '0'
> par7 = '1'
> par6 = '3'
> par5 = '12'
> par4 = '1'
> par3 = '0'
> 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     sma1
      0.8458  -0.0618  0.2003  -0.8214  -0.6357
s.e.  0.1012   0.1095  0.0896   0.0682   0.0978

sigma^2 estimated as 5130319:  log likelihood = -1393.05,  aic = 2798.11
> (forecast <- predict(arima.out,fx))
$pred
Time Series:
Start = 165 
End = 176 
Frequency = 1 
 [1] 26126.11 26559.14 32447.64 38425.53 18037.59 21456.16 24870.93 26369.94
 [9] 25132.01 24857.08 30621.44 28749.42

$se
Time Series:
Start = 165 
End = 176 
Frequency = 1 
 [1] 2265.026 2265.696 2267.615 2297.841 2321.555 2333.756 2347.705 2362.921
 [9] 2377.008 2390.457 2403.686 2416.562

> (lb <- forecast$pred - 1.96 * forecast$se)
Time Series:
Start = 165 
End = 176 
Frequency = 1 
 [1] 21686.66 22118.37 28003.11 33921.76 13487.35 16882.00 20269.43 21738.62
 [9] 20473.07 20171.78 25910.22 24012.96
> (ub <- forecast$pred + 1.96 * forecast$se)
Time Series:
Start = 165 
End = 176 
Frequency = 1 
 [1] 30565.56 30999.90 36892.16 42929.30 22587.84 26030.33 29472.43 31001.27
 [9] 29790.94 29542.37 35332.67 33485.89
> 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] 15136.00 16733.00 20016.00 17708.00 18019.00 19227.00 22893.00 23739.00
  [9] 21133.00 22591.00 26786.00 29740.00 15028.00 17977.00 20008.00 21354.00
 [17] 19498.00 22125.00 25817.00 28779.00 20960.00 22254.00 27392.00 29945.00
 [25] 16933.00 17892.00 20533.00 23569.00 22417.00 22084.00 26580.00 27454.00
 [33] 24081.00 23451.00 28991.00 31386.00 16896.00 20045.00 23471.00 21747.00
 [41] 25621.00 23859.00 25500.00 30998.00 24475.00 23145.00 29701.00 34365.00
 [49] 17556.00 22077.00 25702.00 22214.00 26886.00 23191.00 27831.00 35406.00
 [57] 23195.00 25110.00 30009.00 36242.00 18450.00 21845.00 26488.00 22394.00
 [65] 28057.00 25451.00 24872.00 33424.00 24052.00 28449.00 33533.00 37351.00
 [73] 19969.00 21701.00 26249.00 24493.00 24603.00 26485.00 30723.00 34569.00
 [81] 26689.00 26157.00 32064.00 38870.00 21337.00 19419.00 23166.00 28286.00
 [89] 24570.00 24001.00 33151.00 24878.00 26804.00 28967.00 33311.00 40226.00
 [97] 20504.00 23060.00 23562.00 27562.00 23940.00 24584.00 34303.00 25517.00
[105] 23494.00 29095.00 32903.00 34379.00 16991.00 21109.00 23740.00 25552.00
[113] 21752.00 20294.00 29009.00 25500.00 24166.00 26960.00 31222.00 38641.00
[121] 14672.00 17543.00 25453.00 32683.00 22449.00 22316.00 27595.00 25451.00
[129] 25421.00 25288.00 32568.00 35110.00 16052.00 22146.00 21198.00 19543.00
[137] 22084.00 23816.00 29961.00 26773.00 26635.00 26972.00 30207.00 38687.00
[145] 16974.00 21697.00 24179.00 23757.00 25013.00 24019.00 30345.00 24488.00
[153] 25156.00 25650.00 30923.00 37240.00 17466.00 19463.00 24352.00 26805.00
[161] 25236.00 24735.00 29356.00 31234.00 26126.11 26559.14 32447.64 38425.53
[169] 18037.59 21456.16 24870.93 26369.94 25132.01 24857.08 30621.44 28749.42
> (perc.se <- (ub-forecast$pred)/1.96/forecast$pred)
Time Series:
Start = 165 
End = 176 
Frequency = 1 
 [1] 0.08669586 0.08530760 0.06988536 0.05979984 0.12870647 0.10876854
 [7] 0.09439556 0.08960661 0.09458092 0.09616807 0.07849683 0.08405601
> postscript(file="/var/www/html/rcomp/tmp/1rrwa1198788085.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/2dovh1198788085.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/3hjwk1198788085.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/4o8s81198788085.tab") 
> 
> system("convert tmp/1rrwa1198788085.ps tmp/1rrwa1198788085.png")
> system("convert tmp/2dovh1198788085.ps tmp/2dovh1198788085.png")
> 
> 
> proc.time()
   user  system elapsed 
  1.801   0.361   3.409