R version 2.9.0 (2009-04-17)
Copyright (C) 2009 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(25.2609
+ ,16.8622
+ ,13.3181
+ ,12.5621
+ ,14.2754
+ ,12.2961
+ ,10.0871
+ ,13.5117
+ ,13.9921
+ ,13.6932
+ ,14.4211
+ ,15.3397
+ ,16.5182
+ ,15.2809
+ ,15.6204
+ ,15.5698
+ ,15.9458
+ ,16.4063
+ ,17.55
+ ,17.0353
+ ,16.0591
+ ,16.3643
+ ,14.6527
+ ,13.4664
+ ,13.3266
+ ,13.1823
+ ,12.113
+ ,13.354
+ ,13.4537
+ ,13.2715
+ ,13.1959
+ ,13.5542
+ ,12.124
+ ,10.967
+ ,10.9201
+ ,12.5971
+ ,14.3177
+ ,14.2471
+ ,16.0926
+ ,17.1334
+ ,16.5866
+ ,16.361
+ ,15.8494
+ ,15.5932
+ ,16.6387
+ ,16.8312
+ ,16.5044
+ ,16.5556
+ ,16.7469
+ ,15.9543
+ ,15.5888
+ ,14.3945
+ ,13.8889
+ ,12.9999
+ ,14.1022
+ ,19.6245
+ ,24.7658
+ ,25.9843
+ ,22.9635
+ ,19.6288
+ ,17.3363
+ ,13.311
+ ,14.6359
+ ,15.834
+ ,16.2415
+ ,15.9808
+ ,16.9726
+ ,16.8708
+ ,16.923
+ ,18.1138
+ ,16.7716
+ ,14.0299
+ ,13.822
+ ,14.2537
+ ,14.3985
+ ,15.2454
+ ,15.6683
+ ,16.1721
+ ,14.8679
+ ,14.1948
+ ,14.7056
+ ,15.3819
+ ,15.5001
+ ,14.7886
+ ,14.563
+ ,15.5528
+ ,15.9781
+ ,15.5139
+ ,15.3603
+ ,15.0512
+ ,14.7874
+ ,14.9624
+ ,13.9188
+ ,14.5146
+ ,13.7115
+ ,12.0738
+ ,12.5688
+ ,12.2547
+ ,11.8741
+ ,13.0261
+ ,13.8681
+ ,14.2137
+ ,14.4743
+ ,13.9764
+ ,13.1558
+ ,13.0991
+ ,13.7831
+ ,13.2546
+ ,13.3426
+ ,13.5011
+ ,12.8245
+ ,13.6596
+ ,13.8754
+ ,12.9011
+ ,11.871
+ ,12.3954
+ ,12.8179
+ ,12.1219
+ ,12.6176
+ ,13.6362
+ ,13.5422
+ ,13.362
+ ,14.5735
+ ,15.8357
+ ,14.9927
+ ,14.5078
+ ,15.2648
+ ,15.7163
+ ,17.7969
+ ,19.0408
+ ,17.8571
+ ,18.815
+ ,19.0961
+ ,17.6215
+ ,17.0163
+ ,15.8286
+ ,16.7818
+ ,15.8726
+ ,16.6621
+ ,17.5709
+ ,16.9914
+ ,18.0412
+ ,16.9764
+ ,15.7649
+ ,14.3928
+ ,13.5061
+ ,12.7433
+ ,13.017
+ ,13.0171
+ ,12.2412
+ ,11.8878
+ ,11.2511
+ ,11.8583
+ ,11.1202
+ ,10.185
+ ,8.7563
+ ,9.5267
+ ,9.4106
+ ,11.878
+ ,14.4228
+ ,14.896
+ ,15.6664
+ ,18.147
+ ,19.3069
+ ,21.6807
+ ,20.7934
+ ,23.4241
+ ,24.8273
+ ,24.9276
+ ,27.4256
+ ,28.1746
+ ,24.5615
+ ,30.2532
+ ,31.2514
+ ,30.4733
+ ,33.3047
+ ,37.2103
+ ,36.7711
+ ,37.7163
+ ,28.8488
+ ,27.4682
+ ,29.8793
+ ,28.0598
+ ,29.7733
+ ,32.6926
+ ,32.4803
+ ,29.4168
+ ,28.7054
+ ,28.7614
+ ,23.8075
+ ,21.6987
+ ,21.4691
+ ,22.5709
+ ,23.4546
+ ,27.8976
+ ,29.2965
+ ,28.1191
+ ,25.812
+ ,25.931
+ ,26.9925
+ ,28.9213
+ ,27.8898
+ ,24.2473
+ ,27.1056
+ ,28.2833
+ ,29.8076
+ ,27.1826
+ ,22.8764
+ ,21.938
+ ,23.3076
+ ,24.9572
+ ,26.4694
+ ,23.9297
+ ,24.7033
+ ,24.646
+ ,24.0496
+ ,24.2096
+ ,24.0717
+ ,26.6673
+ ,27.6457
+ ,30.8791
+ ,29.3278
+ ,30.7268
+ ,34.1204
+ ,35.0205
+ ,39.3565
+ ,34.4724
+ ,29.9762
+ ,33.6008
+ ,35.2464
+ ,40.4137
+ ,41.3922
+ ,39.4243
+ ,45.7259
+ ,48.2549
+ ,52.0461
+ ,52.1871
+ ,49.3474
+ ,47.8653
+ ,48.5179
+ ,52.4815
+ ,51.8171
+ ,52.5811
+ ,57.5617
+ ,55.7091
+ ,55.4378
+ ,58.7493
+ ,57.794
+ ,50.282
+ ,47.6976
+ ,46.7381
+ ,47.4282
+ ,42.2269
+ ,44.9066
+ ,47.2648
+ ,50.2325
+ ,50.2504
+ ,52.5685
+ ,55.2325
+ ,52.3674
+ ,55.1692
+ ,57.7252
+ ,62.8232
+ ,62.7599
+ ,62.4387
+ ,64.0862
+ ,66.1209
+ ,69.8474
+ ,80.1039
+ ,85.9319
+ ,85.2843
+ ,77.0383
+ ,69.9981
+ ,55.2039
+ ,43.1188
+ ,32.077
+ ,34.2974
+ ,34.5613
+ ,36.5235
+ ,39.0474
+ ,42.8033
+ ,49.5164
+ ,46.459
+ ,51.1313
+ ,46.9331
+ ,49.7654
+ ,52.0729
+ ,51.6425
+ ,53.9784
+ ,54.4891
+ ,59.0665
+ ,63.9929
+ ,61.6167
+ ,62.1816
+ ,58.9178
+ ,59.9151)
> par10 = 'FALSE'
> par9 = '0'
> par8 = '0'
> par7 = '0'
> par6 = '1'
> par5 = '12'
> par4 = '0'
> par3 = '1'
> par2 = '-0.4'
> par1 = '12'
> #'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
      0.2486
s.e.  0.0597

sigma^2 estimated as 0.0001112:  log likelihood = 886.6,  aic = -1769.21
> (forecast <- predict(arima.out,par1))
$pred
Time Series:
Start = 285 
End = 296 
Frequency = 1 
 [1] 0.2052511 0.2047505 0.2046260 0.2045951 0.2045874 0.2045855 0.2045850
 [8] 0.2045849 0.2045849 0.2045849 0.2045849 0.2045849

$se
Time Series:
Start = 285 
End = 296 
Frequency = 1 
 [1] 0.01054687 0.01687189 0.02181010 0.02590766 0.02945949 0.03263124
 [7] 0.03552188 0.03819459 0.04069219 0.04304513 0.04527596 0.04740191

> (lb <- forecast$pred - 1.96 * forecast$se)
Time Series:
Start = 285 
End = 296 
Frequency = 1 
 [1] 0.1845793 0.1716816 0.1618783 0.1538161 0.1468468 0.1406283 0.1349621
 [8] 0.1297235 0.1248282 0.1202164 0.1158440 0.1116771
> (ub <- forecast$pred + 1.96 * forecast$se)
Time Series:
Start = 285 
End = 296 
Frequency = 1 
 [1] 0.2259230 0.2378194 0.2473738 0.2553741 0.2623280 0.2685427 0.2742079
 [8] 0.2794463 0.2843416 0.2889533 0.2933257 0.2974926
> 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] 25.26090 16.86220 13.31810 12.56210 14.27540 12.29610 10.08710 13.51170
  [9] 13.99210 13.69320 14.42110 15.33970 16.51820 15.28090 15.62040 15.56980
 [17] 15.94580 16.40630 17.55000 17.03530 16.05910 16.36430 14.65270 13.46640
 [25] 13.32660 13.18230 12.11300 13.35400 13.45370 13.27150 13.19590 13.55420
 [33] 12.12400 10.96700 10.92010 12.59710 14.31770 14.24710 16.09260 17.13340
 [41] 16.58660 16.36100 15.84940 15.59320 16.63870 16.83120 16.50440 16.55560
 [49] 16.74690 15.95430 15.58880 14.39450 13.88890 12.99990 14.10220 19.62450
 [57] 24.76580 25.98430 22.96350 19.62880 17.33630 13.31100 14.63590 15.83400
 [65] 16.24150 15.98080 16.97260 16.87080 16.92300 18.11380 16.77160 14.02990
 [73] 13.82200 14.25370 14.39850 15.24540 15.66830 16.17210 14.86790 14.19480
 [81] 14.70560 15.38190 15.50010 14.78860 14.56300 15.55280 15.97810 15.51390
 [89] 15.36030 15.05120 14.78740 14.96240 13.91880 14.51460 13.71150 12.07380
 [97] 12.56880 12.25470 11.87410 13.02610 13.86810 14.21370 14.47430 13.97640
[105] 13.15580 13.09910 13.78310 13.25460 13.34260 13.50110 12.82450 13.65960
[113] 13.87540 12.90110 11.87100 12.39540 12.81790 12.12190 12.61760 13.63620
[121] 13.54220 13.36200 14.57350 15.83570 14.99270 14.50780 15.26480 15.71630
[129] 17.79690 19.04080 17.85710 18.81500 19.09610 17.62150 17.01630 15.82860
[137] 16.78180 15.87260 16.66210 17.57090 16.99140 18.04120 16.97640 15.76490
[145] 14.39280 13.50610 12.74330 13.01700 13.01710 12.24120 11.88780 11.25110
[153] 11.85830 11.12020 10.18500  8.75630  9.52670  9.41060 11.87800 14.42280
[161] 14.89600 15.66640 18.14700 19.30690 21.68070 20.79340 23.42410 24.82730
[169] 24.92760 27.42560 28.17460 24.56150 30.25320 31.25140 30.47330 33.30470
[177] 37.21030 36.77110 37.71630 28.84880 27.46820 29.87930 28.05980 29.77330
[185] 32.69260 32.48030 29.41680 28.70540 28.76140 23.80750 21.69870 21.46910
[193] 22.57090 23.45460 27.89760 29.29650 28.11910 25.81200 25.93100 26.99250
[201] 28.92130 27.88980 24.24730 27.10560 28.28330 29.80760 27.18260 22.87640
[209] 21.93800 23.30760 24.95720 26.46940 23.92970 24.70330 24.64600 24.04960
[217] 24.20960 24.07170 26.66730 27.64570 30.87910 29.32780 30.72680 34.12040
[225] 35.02050 39.35650 34.47240 29.97620 33.60080 35.24640 40.41370 41.39220
[233] 39.42430 45.72590 48.25490 52.04610 52.18710 49.34740 47.86530 48.51790
[241] 52.48150 51.81710 52.58110 57.56170 55.70910 55.43780 58.74930 57.79400
[249] 50.28200 47.69760 46.73810 47.42820 42.22690 44.90660 47.26480 50.23250
[257] 50.25040 52.56850 55.23250 52.36740 55.16920 57.72520 62.82320 62.75990
[265] 62.43870 64.08620 66.12090 69.84740 80.10390 85.93190 85.28430 77.03830
[273] 69.99810 55.20390 43.11880 32.07700 34.29740 34.56130 36.52350 39.04740
[281] 42.80330 49.51640 46.45900 51.13130 52.39456 52.71540 52.79560 52.81556
[289] 52.82053 52.82176 52.82207 52.82215 52.82217 52.82217 52.82217 52.82217
> (perc.se <- (ub-forecast$pred)/1.96/forecast$pred)
Time Series:
Start = 285 
End = 296 
Frequency = 1 
 [1] 0.1550699 0.2822917 0.4063855 0.5308602 0.6587028 0.7922111 0.9332354
 [8] 1.0834030 1.2442669 1.4174022 1.6044725 1.8072826
> postscript(file="/var/www/html/rcomp/tmp/19msz1292173878.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/rcomp/tmp/2gnps1292173878.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/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/3np4m1292173878.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/rcomp/tmp/4qpls1292173878.tab") 
> 
> try(system("convert tmp/19msz1292173878.ps tmp/19msz1292173878.png",intern=TRUE))
character(0)
> try(system("convert tmp/2gnps1292173878.ps tmp/2gnps1292173878.png",intern=TRUE))
character(0)
> 
> 
> proc.time()
   user  system elapsed 
  0.632   0.355   1.365