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.

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(99.4,97.5,94.6,92.6,92.5,89.8,88.8,87.4,85.2,83.1,84.7,84.8,85.8,86.3,89,89,89.3,91.9,94.9,94.4,96.8,96.9,98,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,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,145.4,142.4,137.7,137,137.1,139.3,139.6,140.4,142.3,148.3)
> par10 = 'FALSE'
> par9 = '0'
> par8 = '0'
> par7 = '0'
> par6 = '1'
> par5 = '12'
> par4 = '0'
> par3 = '1'
> par2 = '1'
> par1 = '40'
> #'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
      0.4733
s.e.  0.1989

sigma^2 estimated as 2.531:  log likelihood = -35.91,  aic = 75.82
> (forecast <- predict(arima.out,fx))
$pred
Time Series:
Start = 21 
End = 60 
Frequency = 1 
 [1] 94.16336 94.05137 93.99837 93.97328 93.96141 93.95579 93.95313 93.95188
 [9] 93.95128 93.95100 93.95087 93.95080 93.95077 93.95076 93.95075 93.95075
[17] 93.95075 93.95075 93.95075 93.95075 93.95075 93.95075 93.95075 93.95075
[25] 93.95075 93.95075 93.95075 93.95075 93.95075 93.95075 93.95075 93.95075
[33] 93.95075 93.95075 93.95075 93.95075 93.95075 93.95075 93.95075 93.95075

$se
Time Series:
Start = 21 
End = 60 
Frequency = 1 
 [1]  1.590939  2.832822  3.913582  4.852483  5.678151  6.415638  7.084244
 [8]  7.698282  8.268295  8.802120  9.305661  9.783451 10.239029 10.675209
[15] 11.094265 11.498065 11.888159 12.265855 12.632263 12.988339 13.334910
[22] 13.672700 14.002343 14.324402 14.639378 14.947718 15.249825 15.546062
[29] 15.836759 16.122215 16.402704 16.678477 16.949763 17.216776 17.479710
[36] 17.738747 17.994055 18.245792 18.494102 18.739122

> (lb <- forecast$pred - 1.96 * forecast$se)
Time Series:
Start = 21 
End = 60 
Frequency = 1 
 [1] 91.04512 88.49904 86.32775 84.46242 82.83224 81.38114 80.06802 78.86324
 [9] 77.74542 76.69884 75.71177 74.77524 73.88228 73.02735 72.20599 71.41454
[17] 70.64995 69.90967 69.19151 68.49360 67.81432 67.15225 66.50615 65.87492
[25] 65.25757 64.65322 64.06109 63.48046 62.91070 62.35120 61.80145 61.26093
[33] 60.72921 60.20587 59.69051 59.18280 58.68240 58.18899 57.70231 57.22207
> (ub <- forecast$pred + 1.96 * forecast$se)
Time Series:
Start = 21 
End = 60 
Frequency = 1 
 [1]  97.2816  99.6037 101.6690 103.4842 105.0906 106.5304 107.8383 109.0405
 [9] 110.1571 111.2032 112.1900 113.1264 114.0193 114.8742 115.6955 116.4870
[17] 117.2515 117.9918 118.7100 119.4079 120.0872 120.7492 121.3953 122.0266
[25] 122.6439 123.2483 123.8404 124.4210 124.9908 125.5503 126.1000 126.6406
[33] 127.1723 127.6956 128.2110 128.7187 129.2191 129.7125 130.1992 130.6794
> 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] 99.40000 97.50000 94.60000 92.60000 92.50000 89.80000 88.80000 87.40000
 [9] 85.20000 83.10000 84.70000 84.80000 85.80000 86.30000 89.00000 89.00000
[17] 89.30000 91.90000 94.90000 94.40000 94.16336 94.05137 93.99837 93.97328
[25] 93.96141 93.95579 93.95313 93.95188 93.95128 93.95100 93.95087 93.95080
[33] 93.95077 93.95076 93.95075 93.95075 93.95075 93.95075 93.95075 93.95075
[41] 93.95075 93.95075 93.95075 93.95075 93.95075 93.95075 93.95075 93.95075
[49] 93.95075 93.95075 93.95075 93.95075 93.95075 93.95075 93.95075 93.95075
[57] 93.95075 93.95075 93.95075 93.95075
> (perc.se <- (ub-forecast$pred)/1.96/forecast$pred)
Time Series:
Start = 21 
End = 60 
Frequency = 1 
 [1] 0.01689552 0.03011994 0.04163457 0.05163684 0.06043067 0.06828358
 [7] 0.07540189 0.08193856 0.08800620 0.09368841 0.09904817 0.10413377
[13] 0.10898291 0.11362557 0.11808597 0.12238396 0.12653608 0.13055622
[19] 0.13445623 0.13824626 0.14193512 0.14553051 0.14903919 0.15246714
[25] 0.15581971 0.15910164 0.16231723 0.16547034 0.16856448 0.17160284
[31] 0.17458833 0.17752362 0.18041116 0.18325321 0.18605185 0.18880900
[37] 0.19152648 0.19420593 0.19684891 0.19945687
> postscript(file="/var/www/html/freestat/rcomp/tmp/1j1821229190303.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/freestat/rcomp/tmp/249pq1229190303.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 
> 
> #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/3sd6f1229190303.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/freestat/rcomp/tmp/4virf1229190303.tab") 
> 
> system("convert tmp/1j1821229190303.ps tmp/1j1821229190303.png")
> system("convert tmp/249pq1229190303.ps tmp/249pq1229190303.png")
> 
> 
> proc.time()
   user  system elapsed 
  0.994   0.449   1.122