R version 2.12.0 (2010-10-15)
Copyright (C) 2010 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i486-pc-linux-gnu (32-bit)

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(655362,873127,1107897,1555964,1671159,1493308,2957796,2638691,1305669,1280496,921900,867888,652586,913831,1108544,1555827,1699283,1509458,3268975,2425016,1312703,1365498,934453,775019,651142,843192,1146766,1652601,1465906,1652734,2922334,2702805,1458956,1410363,1019279,936574,708917,885295,1099663,1576220,1487870,1488635,2882530,2677026,1404398,1344370,936865,872705,628151,953712,1160384,1400618,1661511,1495347,2918786,2775677,1407026,1370199,964526,850851,683118,847224,1073256,1514326,1503734,1507712,2865698,2788128,1391596,1366378,946295,859626)
> par10 = 'FALSE'
> par9 = '1'
> par8 = '1'
> par7 = '0'
> par6 = '0'
> par5 = '1'
> par4 = '1'
> par3 = '0'
> par2 = '-0.2'
> par1 = '24'
> #'GNU S' R Code compiled by R2WASP v. 1.0.44 ()
> #Author: Prof. Dr. P. Wessa
> #To cite this work: Wessa P., (2009), ARIMA Forecasting (v1.0.5) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_arimaforecasting.wasp/
> #Source of accompanying publication: 
> #Technical description: 
> 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.5549  0.7010
s.e.   0.3135  0.2568

sigma^2 estimated as 1.489e-05:  log likelihood = 194.48,  aic = -382.97
> (forecast <- predict(arima.out,par1))
$pred
Time Series:
Start = 49 
End = 72 
Frequency = 1 
 [1] 0.06434009 0.06461613 0.06446295 0.06454795 0.06450079 0.06452696
 [7] 0.06451244 0.06452049 0.06451602 0.06451850 0.06451713 0.06451789
[13] 0.06451747 0.06451770 0.06451757 0.06451764 0.06451760 0.06451763
[19] 0.06451761 0.06451762 0.06451762 0.06451762 0.06451762 0.06451762

$se
Time Series:
Start = 49 
End = 72 
Frequency = 1 
 [1] 0.003858243 0.005868492 0.007164081 0.008346581 0.009337618 0.010255103
 [7] 0.011085692 0.011864099 0.012591426 0.013280585 0.013934853 0.014560202
[13] 0.015159536 0.015736189 0.016292375 0.016830229 0.017351398 0.017857374
[19] 0.018349397 0.018828571 0.019295847 0.019752073 0.020197996 0.020634285

> (lb <- forecast$pred - 1.96 * forecast$se)
Time Series:
Start = 49 
End = 72 
Frequency = 1 
 [1] 0.05677794 0.05311389 0.05042136 0.04818865 0.04619905 0.04442696
 [7] 0.04278448 0.04126686 0.03983683 0.03848856 0.03720482 0.03597990
[13] 0.03480478 0.03367477 0.03258452 0.03153040 0.03050886 0.02951717
[19] 0.02855280 0.02761362 0.02669776 0.02580356 0.02492955 0.02407442
> (ub <- forecast$pred + 1.96 * forecast$se)
Time Series:
Start = 49 
End = 72 
Frequency = 1 
 [1] 0.07190225 0.07611837 0.07850455 0.08090725 0.08280252 0.08462696
 [7] 0.08624039 0.08777413 0.08919522 0.09054845 0.09182944 0.09305589
[13] 0.09423016 0.09536063 0.09645063 0.09750489 0.09852634 0.09951808
[19] 0.10048243 0.10142162 0.10233748 0.10323168 0.10410569 0.10496082
> 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]  655362.0  873127.0 1107897.0 1555964.0 1671159.0 1493308.0 2957796.0
 [8] 2638691.0 1305669.0 1280496.0  921900.0  867888.0  652586.0  913831.0
[15] 1108544.0 1555827.0 1699283.0 1509458.0 3268975.0 2425016.0 1312703.0
[22] 1365498.0  934453.0  775019.0  651142.0  843192.0 1146766.0 1652601.0
[29] 1465906.0 1652734.0 2922334.0 2702805.0 1458956.0 1410363.0 1019279.0
[36]  936574.0  708917.0  885295.0 1099663.0 1576220.0 1487870.0 1488635.0
[43] 2882530.0 2677026.0 1404398.0 1344370.0  936865.0  872705.0  906967.1
[50]  887759.5  898357.0  892457.7  895725.5  893910.4  894917.1  894358.3
[57]  894668.3  894496.3  894591.8  894538.8  894568.2  894551.9  894560.9
[64]  894555.9  894558.7  894557.1  894558.0  894557.5  894557.8  894557.6
[71]  894557.7  894557.7
> (perc.se <- (ub-forecast$pred)/1.96/forecast$pred)
Time Series:
Start = 49 
End = 72 
Frequency = 1 
 [1]  0.4431502  0.8493868  1.2324949  1.6898449  2.1962680  2.7875411
 [7]  3.4665606  4.2565769  5.1738485  6.2429578  7.4904082  8.9487964
[13] 10.6567879 12.6606831 15.0164398 17.7914056 21.0674818 24.9441442
[19] 29.5431028 35.0135444 41.5393776 49.3481338 58.7227947 70.0169492
> postscript(file="/var/www/rcomp/tmp/1d0g61323202357.ps",horizontal=F,onefile=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.mape1 <- array(0, dim=fx)
> perf.se <- array(0, dim=fx)
> perf.mse <- array(0, dim=fx)
> perf.mse1 <- 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.se[i] = (x[nx+i] - forecast$pred[i])^2
+ 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[1] = abs(perf.pe[1])
> perf.mse[1] = abs(perf.se[1])
> for (i in 2:fx) {
+ perf.mape[i] = perf.mape[i-1] + abs(perf.pe[i])
+ perf.mape1[i] = perf.mape[i] / i
+ perf.mse[i] = perf.mse[i-1] + perf.se[i]
+ perf.mse1[i] = perf.mse[i] / i
+ }
> perf.rmse = sqrt(perf.mse1)
> postscript(file="/var/www/rcomp/tmp/2cw8u1323202357.ps",horizontal=F,onefile=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:par1] <- x[(nx+1):lx]
> lines(dum, lty=1)
> lines(ub,lty=3)
> lines(lb,lty=3)
> dev.off()
null device 
          1 
> 
> #Note: the /var/www/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/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/rcomp/tmp/3huva1323202357.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.mape1[i],4))
+ a<-table.element(a,round(perf.se[i],4))
+ a<-table.element(a,round(perf.mse1[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/rcomp/tmp/40inq1323202357.tab") 
> 
> try(system("convert tmp/1d0g61323202357.ps tmp/1d0g61323202357.png",intern=TRUE))
character(0)
> try(system("convert tmp/2cw8u1323202357.ps tmp/2cw8u1323202357.png",intern=TRUE))
character(0)
> 
> 
> proc.time()
   user  system elapsed 
  0.680   0.040   0.698