R version 2.13.0 (2011-04-13)
Copyright (C) 2011 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(348542
+ ,335658
+ ,330664
+ ,326814
+ ,322900
+ ,322310
+ ,385164
+ ,404861
+ ,412136
+ ,411057
+ ,410040
+ ,414980
+ ,413626
+ ,411062
+ ,408352
+ ,409780
+ ,411318
+ ,415555
+ ,479481
+ ,497826
+ ,501638
+ ,497990
+ ,499287
+ ,506247
+ ,510401
+ ,508642
+ ,501805
+ ,495476
+ ,490336
+ ,490042
+ ,553155
+ ,569999
+ ,573170
+ ,571687
+ ,575453
+ ,580177
+ ,579849
+ ,574346
+ ,563325
+ ,555604
+ ,545544
+ ,545109
+ ,605181
+ ,627856
+ ,631421
+ ,625671
+ ,613577
+ ,606463
+ ,601676
+ ,589121
+ ,573559
+ ,558487
+ ,552148
+ ,545720
+ ,606569
+ ,636067
+ ,630704
+ ,623275
+ ,617771
+ ,605401
+ ,619393
+ ,596019
+ ,569977
+ ,546213
+ ,528492
+ ,505944
+ ,554910
+ ,567831
+ ,564021
+ ,552800
+ ,541102
+ ,542378
+ ,540380
+ ,521219
+ ,504652
+ ,490626
+ ,481686
+ ,477930
+ ,522605
+ ,531432
+ ,532355
+ ,539954
+ ,524987
+ ,533307
+ ,530541
+ ,508392
+ ,495208
+ ,482223
+ ,470495
+ ,466106
+ ,515037
+ ,517752
+ ,515565
+ ,510727
+ ,499725
+ ,498369
+ ,493756
+ ,476141
+ ,458458
+ ,443182
+ ,429597
+ ,424476
+ ,476257
+ ,480555
+ ,469762
+ ,459820
+ ,451028
+ ,450065
+ ,444385
+ ,428846
+ ,421020
+ ,399778
+ ,389005
+ ,384018
+ ,431933
+ ,445844
+ ,431464
+ ,423263
+ ,415881
+ ,416208
+ ,413491
+ ,399153
+ ,385939
+ ,373917
+ ,364635
+ ,364696
+ ,418358
+ ,428212
+ ,423730
+ ,420677
+ ,417428
+ ,423245
+ ,423113
+ ,418873
+ ,405733
+ ,397812
+ ,389918
+ ,391116
+ ,443814
+ ,460373
+ ,455422
+ ,456288
+ ,452233
+ ,459256
+ ,461146
+ ,451391
+ ,443101
+ ,438810
+ ,430457
+ ,435721
+ ,488280
+ ,505814
+ ,502338
+ ,500910
+ ,501434
+ ,515476
+ ,520862
+ ,519517
+ ,511805
+ ,508607
+ ,505327
+ ,511435
+ ,570158
+ ,591665
+ ,593572
+ ,586346
+ ,586063
+ ,591504
+ ,594033
+ ,585597
+ ,572450
+ ,562917
+ ,554675
+ ,553997
+ ,601310
+ ,622255
+ ,616735
+ ,606480
+ ,595079
+ ,598588
+ ,599917
+ ,591573
+ ,575489
+ ,567223
+ ,555338
+ ,555252
+ ,608249
+ ,630859
+ ,628632
+ ,624435
+ ,609670
+ ,615830
+ ,621170
+ ,604212
+ ,584348
+ ,573717
+ ,555234
+ ,544897
+ ,598866
+ ,620081
+ ,607699
+ ,589960
+ ,578665
+ ,580166
+ ,579457
+ ,571560
+ ,560460
+ ,551397
+ ,536763
+ ,540562
+ ,588184
+ ,607049
+ ,598968
+ ,577644
+ ,562640
+ ,565867
+ ,561274
+ ,554144
+ ,539900
+ ,526271
+ ,511841
+ ,505282
+ ,554083
+ ,584225
+ ,568858
+ ,539516
+ ,521612
+ ,525562
+ ,526519
+ ,515713
+ ,503454
+ ,489301
+ ,479020
+ ,475102
+ ,523682
+ ,551528
+ ,531626
+ ,511037
+ ,492417
+ ,492188
+ ,492865
+ ,480961
+ ,461935
+ ,456608
+ ,441977
+ ,439148
+ ,488180
+ ,520564
+ ,501492
+ ,485025
+ ,464196
+ ,460170
+ ,467037
+ ,460070
+ ,447988
+ ,442867
+ ,436087
+ ,431328
+ ,484015
+ ,509673
+ ,512927
+ ,502831
+ ,470984
+ ,471067
+ ,476049
+ ,474605
+ ,470439
+ ,461251
+ ,454724
+ ,455626
+ ,516847
+ ,525192
+ ,522975
+ ,518585
+ ,509239
+ ,512238
+ ,519164
+ ,517009
+ ,509933
+ ,509127
+ ,500857
+ ,506971
+ ,569323
+ ,579714
+ ,577992
+ ,565464
+ ,547344
+ ,554788
+ ,562325
+ ,560854
+ ,555332
+ ,543599
+ ,536662
+ ,542722
+ ,593530
+ ,610763
+ ,612613
+ ,611324
+ ,594167
+ ,595454
+ ,590865
+ ,589379
+ ,584428
+ ,573100
+ ,567456
+ ,569028
+ ,620735
+ ,628884
+ ,628232
+ ,612117
+ ,595404
+ ,597141)
> par10 = 'FALSE'
> par9 = '1'
> par8 = '1'
> par7 = '1'
> par6 = '1'
> par5 = '12'
> par4 = '1'
> par3 = '1'
> par2 = '1.2'
> par1 = '24'
> 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      ma1    sar1     sma1
      0.9287  -0.7948  0.4056  -0.8739
s.e.  0.0504   0.0807  0.0912   0.0626

sigma^2 estimated as 8.548e+09:  log likelihood = -3539.09,  aic = 7088.18
> (forecast <- predict(arima.out,par1))
$pred
Time Series:
Start = 289 
End = 312 
Frequency = 1 
 [1] 7898169 7816179 7664477 7591417 7445341 7488037 8459177 8719880 8654212
[10] 8455390 8188373 8275846 8341551 8237350 8068608 7968302 7815795 7832060
[19] 8769770 9062201 8979702 8784699 8531876 8602172

$se
Time Series:
Start = 289 
End = 312 
Frequency = 1 
 [1]   92470.54  139797.26  181874.33  221835.71  260692.33  298857.59
 [7]  336516.10  373750.74  410595.43  447059.80  483141.47  518832.64
[13]  573250.14  627590.00  681747.89  735632.82  789168.76  842292.61
[19]  894952.39  947105.67  998718.22 1049762.85 1100218.41 1150068.91

> (lb <- forecast$pred - 1.96 * forecast$se)
Time Series:
Start = 289 
End = 312 
Frequency = 1 
 [1] 7716926 7542177 7308003 7156619 6934384 6902276 7799605 7987329 7849445
[10] 7579153 7241416 7258934 7217980 7007273 6732382 6526462 6269024 6181166
[19] 7015663 7205874 7022214 6727164 6375448 6348037
> (ub <- forecast$pred + 1.96 * forecast$se)
Time Series:
Start = 289 
End = 312 
Frequency = 1 
 [1]  8079411  8090182  8020951  8026215  7956298  8073798  9118748  9452431
 [9]  9458979  9331627  9135330  9292758  9465121  9467426  9404834  9410142
[17]  9362566  9482953 10523877 10918528 10937189 10842235 10688304 10856307
> 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] 348542.0 335658.0 330664.0 326814.0 322900.0 322310.0 385164.0 404861.0
  [9] 412136.0 411057.0 410040.0 414980.0 413626.0 411062.0 408352.0 409780.0
 [17] 411318.0 415555.0 479481.0 497826.0 501638.0 497990.0 499287.0 506247.0
 [25] 510401.0 508642.0 501805.0 495476.0 490336.0 490042.0 553155.0 569999.0
 [33] 573170.0 571687.0 575453.0 580177.0 579849.0 574346.0 563325.0 555604.0
 [41] 545544.0 545109.0 605181.0 627856.0 631421.0 625671.0 613577.0 606463.0
 [49] 601676.0 589121.0 573559.0 558487.0 552148.0 545720.0 606569.0 636067.0
 [57] 630704.0 623275.0 617771.0 605401.0 619393.0 596019.0 569977.0 546213.0
 [65] 528492.0 505944.0 554910.0 567831.0 564021.0 552800.0 541102.0 542378.0
 [73] 540380.0 521219.0 504652.0 490626.0 481686.0 477930.0 522605.0 531432.0
 [81] 532355.0 539954.0 524987.0 533307.0 530541.0 508392.0 495208.0 482223.0
 [89] 470495.0 466106.0 515037.0 517752.0 515565.0 510727.0 499725.0 498369.0
 [97] 493756.0 476141.0 458458.0 443182.0 429597.0 424476.0 476257.0 480555.0
[105] 469762.0 459820.0 451028.0 450065.0 444385.0 428846.0 421020.0 399778.0
[113] 389005.0 384018.0 431933.0 445844.0 431464.0 423263.0 415881.0 416208.0
[121] 413491.0 399153.0 385939.0 373917.0 364635.0 364696.0 418358.0 428212.0
[129] 423730.0 420677.0 417428.0 423245.0 423113.0 418873.0 405733.0 397812.0
[137] 389918.0 391116.0 443814.0 460373.0 455422.0 456288.0 452233.0 459256.0
[145] 461146.0 451391.0 443101.0 438810.0 430457.0 435721.0 488280.0 505814.0
[153] 502338.0 500910.0 501434.0 515476.0 520862.0 519517.0 511805.0 508607.0
[161] 505327.0 511435.0 570158.0 591665.0 593572.0 586346.0 586063.0 591504.0
[169] 594033.0 585597.0 572450.0 562917.0 554675.0 553997.0 601310.0 622255.0
[177] 616735.0 606480.0 595079.0 598588.0 599917.0 591573.0 575489.0 567223.0
[185] 555338.0 555252.0 608249.0 630859.0 628632.0 624435.0 609670.0 615830.0
[193] 621170.0 604212.0 584348.0 573717.0 555234.0 544897.0 598866.0 620081.0
[201] 607699.0 589960.0 578665.0 580166.0 579457.0 571560.0 560460.0 551397.0
[209] 536763.0 540562.0 588184.0 607049.0 598968.0 577644.0 562640.0 565867.0
[217] 561274.0 554144.0 539900.0 526271.0 511841.0 505282.0 554083.0 584225.0
[225] 568858.0 539516.0 521612.0 525562.0 526519.0 515713.0 503454.0 489301.0
[233] 479020.0 475102.0 523682.0 551528.0 531626.0 511037.0 492417.0 492188.0
[241] 492865.0 480961.0 461935.0 456608.0 441977.0 439148.0 488180.0 520564.0
[249] 501492.0 485025.0 464196.0 460170.0 467037.0 460070.0 447988.0 442867.0
[257] 436087.0 431328.0 484015.0 509673.0 512927.0 502831.0 470984.0 471067.0
[265] 476049.0 474605.0 470439.0 461251.0 454724.0 455626.0 516847.0 525192.0
[273] 522975.0 518585.0 509239.0 512238.0 519164.0 517009.0 509933.0 509127.0
[281] 500857.0 506971.0 569323.0 579714.0 577992.0 565464.0 547344.0 554788.0
[289] 559678.5 554832.8 545844.3 541504.9 532807.7 535352.7 592616.0 607797.3
[297] 603980.5 592395.0 576763.8 581893.7 585741.0 579637.2 569725.3 563817.0
[305] 554810.0 555772.0 610693.8 627617.1 622852.1 611560.0 596857.2 600952.4
> (perc.se <- (ub-forecast$pred)/1.96/forecast$pred)
Time Series:
Start = 289 
End = 312 
Frequency = 1 
 [1] 0.009738045 0.014861729 0.019699306 0.024237876 0.029015888 0.033048966
 [7] 0.032941858 0.035476022 0.039241317 0.043694656 0.048715672 0.051732686
[13] 0.056657131 0.062742693 0.069497325 0.075847623 0.082851455 0.088161967
[19] 0.083723576 0.085713136 0.091127181 0.097796285 0.105394322 0.109197430
> postscript(file="/var/wessaorg/rcomp/tmp/1juac1322906587.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/wessaorg/rcomp/tmp/2xn3s1322906587.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/wessaorg/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/wessaorg/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/wessaorg/rcomp/tmp/3mntp1322906587.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/wessaorg/rcomp/tmp/4vp9y1322906587.tab") 
> 
> try(system("convert tmp/1juac1322906587.ps tmp/1juac1322906587.png",intern=TRUE))
character(0)
> try(system("convert tmp/2xn3s1322906587.ps tmp/2xn3s1322906587.png",intern=TRUE))
character(0)
> 
> 
> proc.time()
   user  system elapsed 
  2.491   0.138   2.661