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(117
+ ,116
+ ,166
+ ,180
+ ,202
+ ,290
+ ,298
+ ,441
+ ,388
+ ,260
+ ,175
+ ,105
+ ,137
+ ,142
+ ,176
+ ,231
+ ,240
+ ,316
+ ,363
+ ,537
+ ,487
+ ,324
+ ,185
+ ,133
+ ,169
+ ,157
+ ,206
+ ,244
+ ,243
+ ,393
+ ,405
+ ,579
+ ,525
+ ,373
+ ,198
+ ,148
+ ,201
+ ,177
+ ,222
+ ,275
+ ,290
+ ,402
+ ,534
+ ,614
+ ,578
+ ,419
+ ,203
+ ,173
+ ,229
+ ,192
+ ,294
+ ,310
+ ,365
+ ,509
+ ,537
+ ,655
+ ,643
+ ,444
+ ,259
+ ,229
+ ,276
+ ,245
+ ,324
+ ,323
+ ,349
+ ,480
+ ,530
+ ,676
+ ,670
+ ,476
+ ,281
+ ,240
+ ,259
+ ,237
+ ,400
+ ,367
+ ,497
+ ,593
+ ,696
+ ,969
+ ,878
+ ,581
+ ,373
+ ,232
+ ,358
+ ,318
+ ,410
+ ,480
+ ,604
+ ,713
+ ,844
+ ,1134
+ ,1013
+ ,755
+ ,371
+ ,280
+ ,417
+ ,417
+ ,514
+ ,548
+ ,583
+ ,839
+ ,924
+ ,1179
+ ,1109
+ ,896
+ ,452
+ ,337
+ ,484
+ ,524
+ ,575
+ ,622
+ ,664
+ ,926
+ ,1028
+ ,1361
+ ,1304
+ ,937
+ ,505
+ ,427
+ ,580
+ ,483
+ ,625
+ ,695
+ ,729
+ ,1099
+ ,1090
+ ,1393
+ ,1261
+ ,988
+ ,525
+ ,416
+ ,516
+ ,454
+ ,629
+ ,755
+ ,706
+ ,951
+ ,1099
+ ,1444
+ ,1316
+ ,1066
+ ,585
+ ,430
+ ,669
+ ,598
+ ,714
+ ,835
+ ,912
+ ,1031
+ ,1210
+ ,1581
+ ,1416
+ ,1120
+ ,652
+ ,505
+ ,741
+ ,675
+ ,782
+ ,956
+ ,996
+ ,1259
+ ,1389
+ ,1868
+ ,1609
+ ,1385
+ ,735
+ ,577
+ ,815
+ ,798
+ ,940
+ ,1007
+ ,1094
+ ,1413
+ ,1552
+ ,2038
+ ,1762
+ ,1411
+ ,805
+ ,729
+ ,912
+ ,753
+ ,989
+ ,1137
+ ,1256
+ ,1554
+ ,1629
+ ,2024
+ ,1900
+ ,1563
+ ,905
+ ,766
+ ,952
+ ,915
+ ,1197
+ ,1242
+ ,1197
+ ,1522
+ ,1591
+ ,2128
+ ,1962
+ ,1653
+ ,987
+ ,877
+ ,990
+ ,880
+ ,1258
+ ,1240
+ ,1312
+ ,1713
+ ,1683
+ ,2220
+ ,1996
+ ,1628
+ ,1119
+ ,890
+ ,1118
+ ,1164
+ ,1364
+ ,1412
+ ,1721
+ ,1752
+ ,1794
+ ,2434
+ ,2390
+ ,1929
+ ,1352
+ ,1060
+ ,1435
+ ,1196
+ ,1478
+ ,1648
+ ,1812
+ ,2118
+ ,2211
+ ,2826
+ ,2534
+ ,2290
+ ,1367
+ ,1105
+ ,1463
+ ,1299
+ ,1576
+ ,1850
+ ,1929
+ ,2367
+ ,2508
+ ,3073
+ ,2922
+ ,2377
+ ,1627
+ ,1259
+ ,1547
+ ,1436
+ ,1905
+ ,2079
+ ,1994
+ ,2501
+ ,2569
+ ,3467
+ ,2885
+ ,2211
+ ,1597
+ ,1141
+ ,1533
+ ,1546
+ ,1967
+ ,2171
+ ,2021
+ ,2753
+ ,2626
+ ,3532
+ ,3096
+ ,2639
+ ,1653
+ ,1425
+ ,1802
+ ,1674
+ ,1970
+ ,2092
+ ,2280
+ ,2715
+ ,2971
+ ,3937
+ ,3110
+ ,2662
+ ,1728
+ ,1609
+ ,1922
+ ,1863
+ ,1945
+ ,2365
+ ,2275
+ ,2962
+ ,2930
+ ,4062
+ ,3445
+ ,2943
+ ,1879
+ ,1694
+ ,2147
+ ,1999
+ ,2266
+ ,2562
+ ,2583
+ ,2965
+ ,3142
+ ,4115
+ ,3654
+ ,2992
+ ,2031
+ ,1699
+ ,2313
+ ,1970
+ ,2382
+ ,2830
+ ,2614
+ ,3321
+ ,3418
+ ,4468
+ ,3657
+ ,3250
+ ,2174
+ ,2014
+ ,2118
+ ,2227
+ ,2563
+ ,2817
+ ,2680
+ ,3337
+ ,3559
+ ,4608
+ ,3930
+ ,3133
+ ,2042
+ ,1999
+ ,2679
+ ,2425
+ ,2693
+ ,2760
+ ,2941
+ ,3611
+ ,3779
+ ,4945
+ ,4034
+ ,2906
+ ,2132
+ ,1932
+ ,2268
+ ,2178
+ ,2317
+ ,2552
+ ,2582
+ ,2886
+ ,3283
+ ,4125
+ ,3536
+ ,2568
+ ,1802
+ ,1598
+ ,2013
+ ,1872
+ ,2227
+ ,2497
+ ,2530
+ ,3119
+ ,3411
+ ,4511
+ ,3528
+ ,2833
+ ,1760
+ ,1517
+ ,1968
+ ,1809
+ ,2104
+ ,2391
+ ,2691
+ ,3023
+ ,3188
+ ,4057
+ ,3476)
> par10 = 'FALSE'
> par9 = '1'
> par8 = '1'
> par7 = '1'
> par6 = '0'
> par5 = '12'
> par4 = '1'
> par3 = '1'
> par2 = '0.2'
> 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:
          ma1    sar1     sma1
      -0.6864  0.1583  -0.8198
s.e.   0.0409  0.0733   0.0506

sigma^2 estimated as 0.002676:  log likelihood = 543.75,  aic = -1079.49
> (forecast <- predict(arima.out,par1))
$pred
Time Series:
Start = 370 
End = 381 
Frequency = 1 
 [1] 4.910342 4.549504 4.440579 4.637536 4.578586 4.717035 4.814781 4.823350
 [9] 5.016038 5.084160 5.361312 5.166519

$se
Time Series:
Start = 370 
End = 381 
Frequency = 1 
 [1] 0.05172753 0.05421201 0.05658750 0.05886722 0.06106188 0.06318035
 [7] 0.06523006 0.06721730 0.06914744 0.07102516 0.07285449 0.07463901

> (lb <- forecast$pred - 1.96 * forecast$se)
Time Series:
Start = 370 
End = 381 
Frequency = 1 
 [1] 4.808956 4.443248 4.329667 4.522156 4.458905 4.593201 4.686930 4.691604
 [9] 4.880509 4.944950 5.218517 5.020227
> (ub <- forecast$pred + 1.96 * forecast$se)
Time Series:
Start = 370 
End = 381 
Frequency = 1 
 [1] 5.011728 4.655759 4.551490 4.752915 4.698268 4.840868 4.942632 4.955096
 [9] 5.151567 5.223369 5.504107 5.312812
> 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]  117.000  116.000  166.000  180.000  202.000  290.000  298.000  441.000
  [9]  388.000  260.000  175.000  105.000  137.000  142.000  176.000  231.000
 [17]  240.000  316.000  363.000  537.000  487.000  324.000  185.000  133.000
 [25]  169.000  157.000  206.000  244.000  243.000  393.000  405.000  579.000
 [33]  525.000  373.000  198.000  148.000  201.000  177.000  222.000  275.000
 [41]  290.000  402.000  534.000  614.000  578.000  419.000  203.000  173.000
 [49]  229.000  192.000  294.000  310.000  365.000  509.000  537.000  655.000
 [57]  643.000  444.000  259.000  229.000  276.000  245.000  324.000  323.000
 [65]  349.000  480.000  530.000  676.000  670.000  476.000  281.000  240.000
 [73]  259.000  237.000  400.000  367.000  497.000  593.000  696.000  969.000
 [81]  878.000  581.000  373.000  232.000  358.000  318.000  410.000  480.000
 [89]  604.000  713.000  844.000 1134.000 1013.000  755.000  371.000  280.000
 [97]  417.000  417.000  514.000  548.000  583.000  839.000  924.000 1179.000
[105] 1109.000  896.000  452.000  337.000  484.000  524.000  575.000  622.000
[113]  664.000  926.000 1028.000 1361.000 1304.000  937.000  505.000  427.000
[121]  580.000  483.000  625.000  695.000  729.000 1099.000 1090.000 1393.000
[129] 1261.000  988.000  525.000  416.000  516.000  454.000  629.000  755.000
[137]  706.000  951.000 1099.000 1444.000 1316.000 1066.000  585.000  430.000
[145]  669.000  598.000  714.000  835.000  912.000 1031.000 1210.000 1581.000
[153] 1416.000 1120.000  652.000  505.000  741.000  675.000  782.000  956.000
[161]  996.000 1259.000 1389.000 1868.000 1609.000 1385.000  735.000  577.000
[169]  815.000  798.000  940.000 1007.000 1094.000 1413.000 1552.000 2038.000
[177] 1762.000 1411.000  805.000  729.000  912.000  753.000  989.000 1137.000
[185] 1256.000 1554.000 1629.000 2024.000 1900.000 1563.000  905.000  766.000
[193]  952.000  915.000 1197.000 1242.000 1197.000 1522.000 1591.000 2128.000
[201] 1962.000 1653.000  987.000  877.000  990.000  880.000 1258.000 1240.000
[209] 1312.000 1713.000 1683.000 2220.000 1996.000 1628.000 1119.000  890.000
[217] 1118.000 1164.000 1364.000 1412.000 1721.000 1752.000 1794.000 2434.000
[225] 2390.000 1929.000 1352.000 1060.000 1435.000 1196.000 1478.000 1648.000
[233] 1812.000 2118.000 2211.000 2826.000 2534.000 2290.000 1367.000 1105.000
[241] 1463.000 1299.000 1576.000 1850.000 1929.000 2367.000 2508.000 3073.000
[249] 2922.000 2377.000 1627.000 1259.000 1547.000 1436.000 1905.000 2079.000
[257] 1994.000 2501.000 2569.000 3467.000 2885.000 2211.000 1597.000 1141.000
[265] 1533.000 1546.000 1967.000 2171.000 2021.000 2753.000 2626.000 3532.000
[273] 3096.000 2639.000 1653.000 1425.000 1802.000 1674.000 1970.000 2092.000
[281] 2280.000 2715.000 2971.000 3937.000 3110.000 2662.000 1728.000 1609.000
[289] 1922.000 1863.000 1945.000 2365.000 2275.000 2962.000 2930.000 4062.000
[297] 3445.000 2943.000 1879.000 1694.000 2147.000 1999.000 2266.000 2562.000
[305] 2583.000 2965.000 3142.000 4115.000 3654.000 2992.000 2031.000 1699.000
[313] 2313.000 1970.000 2382.000 2830.000 2614.000 3321.000 3418.000 4468.000
[321] 3657.000 3250.000 2174.000 2014.000 2118.000 2227.000 2563.000 2817.000
[329] 2680.000 3337.000 3559.000 4608.000 3930.000 3133.000 2042.000 1999.000
[337] 2679.000 2425.000 2693.000 2760.000 2941.000 3611.000 3779.000 4945.000
[345] 4034.000 2906.000 2132.000 1932.000 2268.000 2178.000 2317.000 2552.000
[353] 2582.000 2886.000 3283.000 4125.000 3536.000 2568.000 1802.000 1598.000
[361] 2013.000 1872.000 2227.000 2497.000 2530.000 3119.000 3411.000 4511.000
[369] 3528.000 2854.689 1949.037 1726.624 2145.044 2012.135 2335.314 2587.514
[377] 2610.622 3175.440 3397.003 4429.508 3681.208
> (perc.se <- (ub-forecast$pred)/1.96/forecast$pred)
Time Series:
Start = 370 
End = 381 
Frequency = 1 
 [1] 0.05489248 0.06242893 0.06697972 0.06670588 0.07026039 0.07058021
 [7] 0.07143367 0.07359091 0.07275303 0.07378073 0.07166165 0.07644147
> postscript(file="/var/wessaorg/rcomp/tmp/1zup91323791723.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/2lxbc1323791723.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/3j3bv1323791723.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/4t2d41323791723.tab") 
> 
> try(system("convert tmp/1zup91323791723.ps tmp/1zup91323791723.png",intern=TRUE))
character(0)
> try(system("convert tmp/2lxbc1323791723.ps tmp/2lxbc1323791723.png",intern=TRUE))
character(0)
> 
> 
> proc.time()
   user  system elapsed 
  1.599   0.166   2.751