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.

  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(186448,190530,194207,190855,200779,204428,207617,212071,214239,215883,223484,221529,225247,226699,231406,232324,237192,236727,240698,240688,245283,243556,247826,245798,250479,249216,251896,247616,249994,246552,248771,247551,249745,245742,249019,245841,248771,244723,246878,246014,248496,244351,248016,246509,249426,247840,251035,250161,254278,250801,253985,249174,251287,247947,249992,243805,255812,250417,253033,248705,253950,251484,251093,245996,252721,248019,250464,245571,252690,250183,253639,254436,265280,268705,270643,271480)
> par10 = 'FALSE'
> par9 = '1'
> par8 = '0'
> par7 = '0'
> par6 = '3'
> par5 = '4'
> par4 = '1'
> par3 = '1'
> par2 = '2.0'
> 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:
          ar1     ar2     ar3     sma1
      -0.1069  0.0438  0.4563  -0.5052
s.e.   0.1472  0.1452  0.1396   0.1654

sigma^2 estimated as 7.034e+17:  log likelihood = -1033.44,  aic = 2076.88
> (forecast <- predict(arima.out,par1))
$pred
Time Series:
Start = 53 
End = 76 
Frequency = 1 
 [1] 63824360227 62211235263 62984152586 61561804593 63220488043 61278736239
 [7] 62538939022 61014517820 62555508539 60844226222 62028060117 60468196930
[13] 62114799164 60355826793 61533213916 60020143455 61639699756 59882727294
[19] 61080068273 59552611034 61175492260 59426639446 60616693256 59091888047

$se
Time Series:
Start = 53 
End = 76 
Frequency = 1 
 [1]   838714985  1124533513  1377424944  1806807246  2351812557  2793271903
 [7]  3265907252  3741699016  4358665016  4919971051  5455039393  6016201708
[13]  6703768088  7325805696  7937481472  8571879968  9310747120  9999363215
[19] 10677748984 11373586019 12167593718 12915946187 13653849567 14409009199

> (lb <- forecast$pred - 1.96 * forecast$se)
Time Series:
Start = 53 
End = 76 
Frequency = 1 
 [1] 62180478856 60007149578 60284399696 58020462391 58610935430 55803923309
 [7] 56137760808 53680787749 54012525107 51201082961 51336182906 48676441582
[13] 48975413712 45997247629 45975750231 43219258718 43390635401 40283975393
[19] 40151680264 37260382437 37327008574 34111384919 33855148106 30850230018
> (ub <- forecast$pred + 1.96 * forecast$se)
Time Series:
Start = 53 
End = 76 
Frequency = 1 
 [1] 65468241599 64415320948 65683905475 65103146796 67830040655 66753549169
 [7] 68940117236 68348247891 71098491971 70487369483 72719937327 72259952279
[13] 75254184616 74714405957 77090677602 76821028192 79888764111 79481479195
[19] 82008456281 81844839631 85023975947 84741893973 87378238407 87333546076
> 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] 186448.0 190530.0 194207.0 190855.0 200779.0 204428.0 207617.0 212071.0
 [9] 214239.0 215883.0 223484.0 221529.0 225247.0 226699.0 231406.0 232324.0
[17] 237192.0 236727.0 240698.0 240688.0 245283.0 243556.0 247826.0 245798.0
[25] 250479.0 249216.0 251896.0 247616.0 249994.0 246552.0 248771.0 247551.0
[33] 249745.0 245742.0 249019.0 245841.0 248771.0 244723.0 246878.0 246014.0
[41] 248496.0 244351.0 248016.0 246509.0 249426.0 247840.0 251035.0 250161.0
[49] 254278.0 250801.0 253985.0 249174.0 252634.8 249421.8 250966.4 248116.5
[57] 251436.8 247545.4 250077.9 247011.2 250111.0 246666.2 249054.3 245902.8
[65] 249228.4 245674.2 248058.9 244990.1 248273.4 244709.5 247143.8 244034.0
[73] 247336.8 243775.8 246204.6 243088.2
> (perc.se <- (ub-forecast$pred)/1.96/forecast$pred)
Time Series:
Start = 53 
End = 76 
Frequency = 1 
 [1] 0.006528721 0.008959361 0.010819964 0.014469562 0.018272863 0.022304007
 [7] 0.025474994 0.029792524 0.033723828 0.038944535 0.042225056 0.047532659
[13] 0.051376023 0.057453576 0.060866858 0.067008079 0.070635926 0.077591224
[19] 0.080981033 0.087917091 0.091282436 0.099055581 0.102357070 0.110051302
> postscript(file="/var/www/html/freestat/rcomp/tmp/1kurn1292433136.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/html/freestat/rcomp/tmp/2rvog1292433136.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/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/3gela1292433136.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/html/freestat/rcomp/tmp/4jf2y1292433136.tab") 
> 
> try(system("convert tmp/1kurn1292433136.ps tmp/1kurn1292433136.png",intern=TRUE))
character(0)
> try(system("convert tmp/2rvog1292433136.ps tmp/2rvog1292433136.png",intern=TRUE))
character(0)
> 
> 
> proc.time()
   user  system elapsed 
  0.966   0.440   1.034