R version 2.6.1 (2007-11-26)
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.

  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(96.9,98.0,97.9,100.9,103.9,103.1,102.5,104.3,102.6,101.7,102.8,105.4,110.9,113.5,116.3,124.0,128.8,133.5,132.6,128.4,127.3,126.7,123.3,123.2,124.4,128.2,128.7,135.7,139.0,145.4,142.4,137.7,137.0,137.1,139.3,139.6,140.4,142.3,148.3,157.7,161.6,161.7,171.8,185.1,176.7,184.4,183.0,180.9,187.0,189.9,193.8,194.5,198.7,204.7,213.2,214.7,211.0,213.2,206.3,210.8)
> par10 = 'FALSE'
> par9 = '1'
> par8 = '2'
> par7 = '2'
> par6 = '3'
> par5 = '12'
> par4 = '0'
> par3 = '1'
> par2 = '-0.6'
> par1 = '36'
> #'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    ma2    sar1     sar2    sma1
      -0.9091  -0.1811  0.4077  1.8007  1.000  1.2745  -0.9995  0.2167
s.e.   0.1030   0.1031  0.0782  0.0070  0.007  0.0117   0.0001  0.7134

sigma^2 estimated as 2.17e-10:  log likelihood = 139.19,  aic = -260.38
> (forecast <- predict(arima.out,fx))
$pred
Time Series:
Start = 25 
End = 60 
Frequency = 1 
 [1] 0.05481175 0.05420198 0.05308097 0.05145528 0.05095600 0.04919943
 [7] 0.04925716 0.05122704 0.05097586 0.05084422 0.05239421 0.05335148
[13] 0.05409249 0.05413369 0.05355346 0.05365512 0.05426758 0.05318150
[19] 0.05303939 0.05451183 0.05391076 0.05358807 0.05466322 0.05585620
[25] 0.05766008 0.05832210 0.05870303 0.06045756 0.06173720 0.06210871
[31] 0.06186989 0.06177759 0.06126258 0.06098287 0.06080391 0.06136759

$se
Time Series:
Start = 25 
End = 60 
Frequency = 1 
 [1] 2.087918e-05 3.809170e-05 4.586200e-05 5.681825e-05 6.665940e-05
 [6] 7.330091e-05 8.196502e-05 8.820879e-05 9.444704e-05 1.010286e-04
[11] 1.063002e-04 1.123291e-04 1.292023e-04 1.521459e-04 1.697047e-04
[16] 1.884542e-04 2.071884e-04 2.219557e-04 2.383592e-04 2.525268e-04
[21] 2.658162e-04 2.796231e-04 2.917167e-04 3.042407e-04 3.206588e-04
[26] 3.391401e-04 3.564906e-04 3.735634e-04 3.907857e-04 4.064027e-04
[31] 4.222024e-04 4.372641e-04 4.515591e-04 4.659376e-04 4.795320e-04
[36] 4.929643e-04

> (lb <- forecast$pred - 1.96 * forecast$se)
Time Series:
Start = 25 
End = 60 
Frequency = 1 
 [1] 0.05477082 0.05412732 0.05299108 0.05134392 0.05082535 0.04905576
 [7] 0.04909651 0.05105415 0.05079074 0.05064621 0.05218586 0.05313131
[13] 0.05383925 0.05383549 0.05322084 0.05328575 0.05386149 0.05274646
[19] 0.05257220 0.05401688 0.05338976 0.05304001 0.05409145 0.05525989
[25] 0.05703158 0.05765738 0.05800431 0.05972537 0.06097126 0.06131216
[31] 0.06104237 0.06092055 0.06037752 0.06006963 0.05986403 0.06040138
> (ub <- forecast$pred + 1.96 * forecast$se)
Time Series:
Start = 25 
End = 60 
Frequency = 1 
 [1] 0.05485267 0.05427664 0.05317086 0.05156665 0.05108665 0.04934309
 [7] 0.04941781 0.05139993 0.05116098 0.05104224 0.05260255 0.05357164
[13] 0.05434572 0.05443190 0.05388608 0.05402449 0.05467366 0.05361653
[19] 0.05350657 0.05500678 0.05443176 0.05413613 0.05523498 0.05645251
[25] 0.05828857 0.05898681 0.05940175 0.06118974 0.06250314 0.06290526
[31] 0.06269741 0.06263463 0.06214763 0.06189611 0.06174380 0.06233380
> 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)
+ }
> (actandfor <- c(x[1:nx], forecast$pred))
 [1]  96.9000  98.0000  97.9000 100.9000 103.9000 103.1000 102.5000 104.3000
 [9] 102.6000 101.7000 102.8000 105.4000 110.9000 113.5000 116.3000 124.0000
[17] 128.8000 133.5000 132.6000 128.4000 127.3000 126.7000 123.3000 123.2000
[25] 126.4379 128.8175 133.3834 140.4807 142.7823 151.3793 151.0837 141.5255
[33] 142.6896 143.3059 136.3101 132.2582 129.2523 129.0884 131.4279 131.0131
[41] 128.5581 132.9635 133.5578 127.5994 129.9793 131.2864 127.0110 122.5221
[49] 116.2005 114.0105 112.7801 107.3781 103.6944 102.6627 103.3240 103.5814
[57] 105.0367 105.8409 106.3606 104.7373
> (perc.se <- (ub-forecast$pred)/1.96/forecast$pred)
Time Series:
Start = 25 
End = 60 
Frequency = 1 
 [1] -0.0006342442 -0.0011691413 -0.0014367566 -0.0018350795 -0.0021728622
 [6] -0.0024734884 -0.0027613583 -0.0028570034 -0.0030730806 -0.0032945876
[11] -0.0033635823 -0.0034898802 -0.0039562005 -0.0046500880 -0.0052380646
[16] -0.0058005915 -0.0063002622 -0.0068808030 -0.0074029785 -0.0076284116
[21] -0.0081131311 -0.0085795573 -0.0087718851 -0.0089505383 -0.0091357237
[26] -0.0095463325 -0.0099629787 -0.0101343622 -0.0103778057 -0.0107220744
[31] -0.0111738330 -0.0115821631 -0.0120522716 -0.0124844206 -0.0128783418
[36] -0.0131125447
> postscript(file="/var/www/html/rcomp/tmp/1o71r1197122204.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/2xolj1197122204.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]
Warning message:
In NextMethod("[<-") :
  number of items to replace is not a multiple of replacement length
> 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/31myu1197122205.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/48fke1197122205.tab") 
> 
> system("convert tmp/1o71r1197122204.ps tmp/1o71r1197122204.png")
> system("convert tmp/2xolj1197122204.ps tmp/2xolj1197122204.png")
> 
> 
> proc.time()
   user  system elapsed 
  6.870   1.329   7.769