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(51.220
+ ,50.487
+ ,49.415
+ ,49.398
+ ,48.196
+ ,47.348
+ ,49.331
+ ,49.644
+ ,49.588
+ ,49.567
+ ,49.010
+ ,49.563
+ ,49.741
+ ,49.487
+ ,48.278
+ ,47.478
+ ,46.985
+ ,45.216
+ ,46.581
+ ,49.266
+ ,48.121
+ ,46.412
+ ,46.285
+ ,46.824
+ ,46.949
+ ,45.355
+ ,44.924
+ ,45.059
+ ,44.202
+ ,44.149
+ ,46.151
+ ,47.703
+ ,48.436
+ ,47.089
+ ,47.492
+ ,49.295
+ ,49.127
+ ,50.041
+ ,48.857
+ ,48.428
+ ,48.788
+ ,48.820
+ ,50.743
+ ,52.590
+ ,51.959
+ ,53.451
+ ,55.674
+ ,56.120
+ ,55.685
+ ,56.714
+ ,54.882
+ ,55.173
+ ,53.574
+ ,53.954
+ ,58.055
+ ,61.062
+ ,58.353
+ ,59.693
+ ,58.833
+ ,60.417
+ ,61.696
+ ,62.515
+ ,62.687
+ ,61.794
+ ,63.014
+ ,63.134
+ ,68.057
+ ,67.327
+ ,68.310
+ ,69.780
+ ,69.944
+ ,69.881
+ ,71.397
+ ,70.631
+ ,70.452
+ ,69.862
+ ,69.114
+ ,69.358
+ ,71.133
+ ,73.128
+ ,73.528
+ ,73.677
+ ,72.273
+ ,71.962
+ ,73.654
+ ,73.305
+ ,73.355
+ ,73.346
+ ,72.881
+ ,72.424
+ ,74.540
+ ,74.847
+ ,75.904
+ ,76.870
+ ,76.370
+ ,77.631
+ ,78.335
+ ,77.926
+ ,77.236
+ ,76.755
+ ,74.710
+ ,73.486
+ ,76.034
+ ,76.389
+ ,77.767
+ ,78.124
+ ,76.696
+ ,77.375
+ ,77.431
+ ,77.347
+ ,77.013
+ ,76.666
+ ,75.225
+ ,75.579
+ ,77.100
+ ,78.592
+ ,79.502
+ ,78.528
+ ,77.775
+ ,77.271
+ ,78.738
+ ,77.885
+ ,76.896
+ ,75.813
+ ,74.958
+ ,75.340
+ ,77.187
+ ,78.602
+ ,81.653
+ ,78.125
+ ,76.092
+ ,74.870
+ ,75.615
+ ,74.776
+ ,72.528
+ ,71.894
+ ,71.641
+ ,71.145
+ ,73.320
+ ,72.186
+ ,72.854
+ ,74.243
+ ,74.628
+ ,72.368
+ ,75.361
+ ,72.746
+ ,70.536
+ ,69.410
+ ,66.219
+ ,66.739
+ ,67.626
+ ,70.602
+ ,71.758
+ ,71.786
+ ,69.641
+ ,68.055
+ ,70.148
+ ,69.390
+ ,68.562
+ ,68.622
+ ,68.120
+ ,68.308
+ ,70.421
+ ,69.766
+ ,72.157
+ ,72.928
+ ,75.340
+ ,74.812
+ ,74.593
+ ,76.003
+ ,75.112
+ ,75.452
+ ,75.634
+ ,75.653
+ ,78.645
+ ,73.100
+ ,79.699
+ ,82.848
+ ,81.834
+ ,81.736
+ ,82.267
+ ,84.120
+ ,83.819
+ ,82.734
+ ,81.842
+ ,81.735
+ ,83.227
+ ,81.934
+ ,89.521
+ ,88.827
+ ,85.874
+ ,85.211
+ ,87.130
+ ,88.620
+ ,89.563
+ ,89.056
+ ,88.542
+ ,89.504
+ ,89.428
+ ,86.040
+ ,96.240
+ ,94.423
+ ,93.028
+ ,92.285
+ ,91.685
+ ,94.260
+ ,93.858
+ ,92.437
+ ,92.980
+ ,92.099
+ ,92.803
+ ,88.551
+ ,98.334
+ ,98.329
+ ,96.455
+ ,97.109
+ ,97.687
+ ,98.512
+ ,98.673
+ ,96.028
+ ,98.014
+ ,95.580
+ ,97.838
+ ,97.760
+ ,99.913
+ ,97.588
+ ,93.942
+ ,93.656
+ ,93.365
+ ,92.881
+ ,93.120
+ ,91.063
+ ,90.930
+ ,91.946
+ ,94.624
+ ,95.484
+ ,95.862
+ ,95.530
+ ,94.574
+ ,94.677
+ ,93.845
+ ,91.533
+ ,91.214
+ ,90.922
+ ,89.563
+ ,89.945
+ ,91.850
+ ,92.505
+ ,92.437
+ ,93.876)
> par10 = 'FALSE'
> par9 = '1'
> par8 = '1'
> par7 = '1'
> par6 = '0'
> par5 = '12'
> par4 = '1'
> par3 = '1'
> par2 = '1'
> 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:
          ma1    sar1     sma1
      -0.2467  0.3790  -0.8552
s.e.   0.0693  0.1073   0.0758

sigma^2 estimated as 2.087:  log likelihood = -406.17,  aic = 820.34
> (forecast <- predict(arima.out,fx))
$pred
Time Series:
Start = 239 
End = 250 
Frequency = 1 
 [1] 93.07137 92.78879 93.01773 93.03093 92.91728 91.40981 91.22994 91.56944
 [9] 93.79277 93.74481 95.98295 95.64901

$se
Time Series:
Start = 239 
End = 250 
Frequency = 1 
 [1] 1.445001 1.809134 2.111371 2.375431 2.612940 2.830590 3.032660 3.222082
 [9] 3.400970 3.570907 3.733117 3.888566

> (lb <- forecast$pred - 1.96 * forecast$se)
Time Series:
Start = 239 
End = 250 
Frequency = 1 
 [1] 90.23916 89.24289 88.87945 88.37508 87.79591 85.86186 85.28593 85.25416
 [9] 87.12687 86.74583 88.66604 88.02742
> (ub <- forecast$pred + 1.96 * forecast$se)
Time Series:
Start = 239 
End = 250 
Frequency = 1 
 [1]  95.90357  96.33470  97.15602  97.68677  98.03864  96.95777  97.17396
 [8]  97.88472 100.45867 100.74378 103.29986 103.27060
> 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] 51.22000 50.48700 49.41500 49.39800 48.19600 47.34800 49.33100 49.64400
  [9] 49.58800 49.56700 49.01000 49.56300 49.74100 49.48700 48.27800 47.47800
 [17] 46.98500 45.21600 46.58100 49.26600 48.12100 46.41200 46.28500 46.82400
 [25] 46.94900 45.35500 44.92400 45.05900 44.20200 44.14900 46.15100 47.70300
 [33] 48.43600 47.08900 47.49200 49.29500 49.12700 50.04100 48.85700 48.42800
 [41] 48.78800 48.82000 50.74300 52.59000 51.95900 53.45100 55.67400 56.12000
 [49] 55.68500 56.71400 54.88200 55.17300 53.57400 53.95400 58.05500 61.06200
 [57] 58.35300 59.69300 58.83300 60.41700 61.69600 62.51500 62.68700 61.79400
 [65] 63.01400 63.13400 68.05700 67.32700 68.31000 69.78000 69.94400 69.88100
 [73] 71.39700 70.63100 70.45200 69.86200 69.11400 69.35800 71.13300 73.12800
 [81] 73.52800 73.67700 72.27300 71.96200 73.65400 73.30500 73.35500 73.34600
 [89] 72.88100 72.42400 74.54000 74.84700 75.90400 76.87000 76.37000 77.63100
 [97] 78.33500 77.92600 77.23600 76.75500 74.71000 73.48600 76.03400 76.38900
[105] 77.76700 78.12400 76.69600 77.37500 77.43100 77.34700 77.01300 76.66600
[113] 75.22500 75.57900 77.10000 78.59200 79.50200 78.52800 77.77500 77.27100
[121] 78.73800 77.88500 76.89600 75.81300 74.95800 75.34000 77.18700 78.60200
[129] 81.65300 78.12500 76.09200 74.87000 75.61500 74.77600 72.52800 71.89400
[137] 71.64100 71.14500 73.32000 72.18600 72.85400 74.24300 74.62800 72.36800
[145] 75.36100 72.74600 70.53600 69.41000 66.21900 66.73900 67.62600 70.60200
[153] 71.75800 71.78600 69.64100 68.05500 70.14800 69.39000 68.56200 68.62200
[161] 68.12000 68.30800 70.42100 69.76600 72.15700 72.92800 75.34000 74.81200
[169] 74.59300 76.00300 75.11200 75.45200 75.63400 75.65300 78.64500 73.10000
[177] 79.69900 82.84800 81.83400 81.73600 82.26700 84.12000 83.81900 82.73400
[185] 81.84200 81.73500 83.22700 81.93400 89.52100 88.82700 85.87400 85.21100
[193] 87.13000 88.62000 89.56300 89.05600 88.54200 89.50400 89.42800 86.04000
[201] 96.24000 94.42300 93.02800 92.28500 91.68500 94.26000 93.85800 92.43700
[209] 92.98000 92.09900 92.80300 88.55100 98.33400 98.32900 96.45500 97.10900
[217] 97.68700 98.51200 98.67300 96.02800 98.01400 95.58000 97.83800 97.76000
[225] 99.91300 97.58800 93.94200 93.65600 93.36500 92.88100 93.12000 91.06300
[233] 90.93000 91.94600 94.62400 95.48400 95.86200 95.53000 93.07137 92.78879
[241] 93.01773 93.03093 92.91728 91.40981 91.22994 91.56944 93.79277 93.74481
[249] 95.98295 95.64901
> (perc.se <- (ub-forecast$pred)/1.96/forecast$pred)
Time Series:
Start = 239 
End = 250 
Frequency = 1 
 [1] 0.01552573 0.01949733 0.02269858 0.02553377 0.02812114 0.03096593
 [7] 0.03324193 0.03518730 0.03626047 0.03809179 0.03889355 0.04065454
> postscript(file="/var/www/html/rcomp/tmp/1ajip1229089522.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/23zi71229089522.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 
> 
> #Note: the /var/www/html/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> 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/37pfl1229089523.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/4yfa81229089523.tab") 
> 
> system("convert tmp/1ajip1229089522.ps tmp/1ajip1229089522.png")
> system("convert tmp/23zi71229089522.ps tmp/23zi71229089522.png")
> 
> 
> proc.time()
   user  system elapsed 
  1.447   0.321   3.502