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(88.74,88.92,88.77,89.17,89.61,89.52,89.74,89.40,89.36,89.38,89.36,89.29,89.59,89.79,89.86,90.21,90.37,90.19,90.33,90.22,90.42,90.54,90.73,91.02,91.19,91.53,91.88,92.06,92.32,92.67,92.85,92.82,93.46,93.23,93.54,93.29,93.20,93.60,93.81,94.62,95.22,95.38,95.31,95.30,95.57,95.42,95.53,95.33,95.90,96.06,96.31,96.34,96.49,96.22,96.53,96.50,96.77,96.66,96.58,96.63,97.06,97.73,98.01,97.76,97.49,97.77,97.96,98.23,98.51,98.19,98.37,98.31,98.60,98.97,99.11,99.64,100.03,99.98,100.32,100.44,100.51,101.00,100.88,100.55,100.83,101.51,102.16,102.39,102.54,102.85,103.47,103.57,103.69,103.50,103.47,103.45,103.48,103.93,103.89,104.40,104.79,104.77,105.13,105.26,104.96)
> par10 = 'FALSE'
> par9 = '1'
> par8 = '1'
> par7 = '0'
> par6 = '0'
> par5 = '12'
> par4 = '0'
> par3 = '1'
> par2 = '0.0'
> 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:
        sar1     sma1
      0.9762  -0.8162
s.e.  0.0488   0.1907

sigma^2 estimated as 6.314e-06:  log likelihood = 414.57,  aic = -823.13
> (forecast <- predict(arima.out,fx))
$pred
Time Series:
Start = 94 
End = 105 
Frequency = 1 
 [1] 4.641369 4.641926 4.641059 4.643512 4.647108 4.649423 4.651744 4.653579
 [9] 4.654311 4.656717 4.656981 4.658788

$se
Time Series:
Start = 94 
End = 105 
Frequency = 1 
 [1] 0.002528512 0.003575856 0.004379512 0.005057025 0.005651538 0.006189205
 [7] 0.006683759 0.007144160 0.007576635 0.007985724 0.008374853 0.008746688

> (lb <- forecast$pred - 1.96 * forecast$se)
Time Series:
Start = 94 
End = 105 
Frequency = 1 
 [1] 4.636413 4.634917 4.632475 4.633600 4.636031 4.637293 4.638643 4.639576
 [9] 4.639460 4.641065 4.640566 4.641644
> (ub <- forecast$pred + 1.96 * forecast$se)
Time Series:
Start = 94 
End = 105 
Frequency = 1 
 [1] 4.646325 4.648935 4.649643 4.653424 4.658185 4.661554 4.664844 4.667581
 [9] 4.669161 4.672369 4.673395 4.675931
> 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]  88.7400  88.9200  88.7700  89.1700  89.6100  89.5200  89.7400  89.4000
  [9]  89.3600  89.3800  89.3600  89.2900  89.5900  89.7900  89.8600  90.2100
 [17]  90.3700  90.1900  90.3300  90.2200  90.4200  90.5400  90.7300  91.0200
 [25]  91.1900  91.5300  91.8800  92.0600  92.3200  92.6700  92.8500  92.8200
 [33]  93.4600  93.2300  93.5400  93.2900  93.2000  93.6000  93.8100  94.6200
 [41]  95.2200  95.3800  95.3100  95.3000  95.5700  95.4200  95.5300  95.3300
 [49]  95.9000  96.0600  96.3100  96.3400  96.4900  96.2200  96.5300  96.5000
 [57]  96.7700  96.6600  96.5800  96.6300  97.0600  97.7300  98.0100  97.7600
 [65]  97.4900  97.7700  97.9600  98.2300  98.5100  98.1900  98.3700  98.3100
 [73]  98.6000  98.9700  99.1100  99.6400 100.0300  99.9800 100.3200 100.4400
 [81] 100.5100 101.0000 100.8800 100.5500 100.8300 101.5100 102.1600 102.3900
 [89] 102.5400 102.8500 103.4700 103.5700 103.6900 103.6862 103.7440 103.6540
 [97] 103.9087 104.2830 104.5247 104.7675 104.9599 105.0368 105.2899 105.3176
[105] 105.5081
> (perc.se <- (ub-forecast$pred)/1.96/forecast$pred)
Time Series:
Start = 94 
End = 105 
Frequency = 1 
 [1] 0.002534788 0.003588417 0.004398362 0.005082170 0.005682955 0.006226897
 [7] 0.006727730 0.007194412 0.007633172 0.008048547 0.008443966 0.008822092
> postscript(file="/var/www/html/rcomp/tmp/1ywz61197317613.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/2x2mg1197317613.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/316wi1197317614.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/4aznz1197317614.tab") 
> 
> system("convert tmp/1ywz61197317613.ps tmp/1ywz61197317613.png")
> system("convert tmp/2x2mg1197317613.ps tmp/2x2mg1197317613.png")
> 
> 
> proc.time()
   user  system elapsed 
  1.842   0.515   1.931