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(1145.11
+ ,1176.86
+ ,1206.41
+ ,1192.72
+ ,1214.82
+ ,1199.07
+ ,1157.47
+ ,1100.1
+ ,1095.63
+ ,1105.63
+ ,1137.79
+ ,1124.72
+ ,1152.6
+ ,1211.85
+ ,1239.62
+ ,1244.13
+ ,1198.42
+ ,1227.99
+ ,1304.92
+ ,1340.26
+ ,1307.32
+ ,1356.51
+ ,1383.29
+ ,1437.87
+ ,1494.56
+ ,1521.42
+ ,1498.76
+ ,1488.75
+ ,1524.62
+ ,1439.27
+ ,1423.11
+ ,1466.85
+ ,1425.83
+ ,1363.45
+ ,1389.18
+ ,1395.89
+ ,1368.43
+ ,1349.03
+ ,1299.88
+ ,1365.41
+ ,1451.04
+ ,1433.75
+ ,1464.65
+ ,1475.57
+ ,1471.16
+ ,1429.12
+ ,1452.46
+ ,1538.09
+ ,1631.59
+ ,1665.5
+ ,1690.6
+ ,1711.74
+ ,1734.1
+ ,1748.09
+ ,1703.45
+ ,1745.74
+ ,1751.01
+ ,1795.65
+ ,1852.13
+ ,1877.1
+ ,1989.31
+ ,2097.76
+ ,2154.87
+ ,2152.18
+ ,2250.27
+ ,2346.9
+ ,2525.56
+ ,2409.36
+ ,2394.36
+ ,2401.33
+ ,2354.32
+ ,2450.41
+ ,2504.67
+ ,2661.39
+ ,2880.4
+ ,3064.42
+ ,3141.12
+ ,3327.7
+ ,3564.95
+ ,3403.13
+ ,3149.9
+ ,3006.84
+ ,3230.66
+ ,3361.13
+ ,3484.74
+ ,3411.13
+ ,3288.18
+ ,3280.37
+ ,3173.95
+ ,3165.26
+ ,3092.71
+ ,3053.05
+ ,3181.96
+ ,2999.93
+ ,3249.57
+ ,3210.52
+ ,3030.29
+ ,2803.47
+ ,2767.63
+ ,2882.6
+ ,2863.36
+ ,2897.06
+ ,3012.61
+ ,3142.95
+ ,3032.93
+ ,3045.78
+ ,3110.52
+ ,3013.24
+ ,2987.1
+ ,2995.55
+ ,2833.18
+ ,2848.96
+ ,2794.83
+ ,2845.26
+ ,2915.02
+ ,2892.63
+ ,2604.42
+ ,2641.65
+ ,2659.81
+ ,2638.53
+ ,2720.25
+ ,2745.88
+ ,2735.7
+ ,2811.7
+ ,2799.43
+ ,2555.28
+ ,2304.98
+ ,2214.95
+ ,2065.81
+ ,1940.49
+ ,2042
+ ,1995.37
+ ,1946.81
+ ,1765.9
+ ,1635.25
+ ,1833.42
+ ,1910.43
+ ,1959.67
+ ,1969.6
+ ,2061.41
+ ,2093.48
+ ,2120.88
+ ,2174.56
+ ,2196.72
+ ,2350.44
+ ,2440.25
+ ,2408.64
+ ,2472.81
+ ,2407.6
+ ,2454.62
+ ,2448.05
+ ,2497.84
+ ,2645.64
+ ,2756.76
+ ,2849.27
+ ,2921.44
+ ,2981.85
+ ,3080.58
+ ,3106.22
+ ,3119.31
+ ,3061.26
+ ,3097.31
+ ,3161.69
+ ,3257.16
+ ,3277.01
+ ,3295.32
+ ,3363.99
+ ,3494.17
+ ,3667.03
+ ,3813.06
+ ,3917.96
+ ,3895.51
+ ,3801.06
+ ,3570.12
+ ,3701.61
+ ,3862.27
+ ,3970.1
+ ,4138.52
+ ,4199.75
+ ,4290.89
+ ,4443.91
+ ,4502.64
+ ,4356.98
+ ,4591.27
+ ,4696.96
+ ,4621.4
+ ,4562.84
+ ,4202.52
+ ,4296.49
+ ,4435.23
+ ,4105.18
+ ,4116.68
+ ,3844.49
+ ,3720.98
+ ,3674.4
+ ,3857.62
+ ,3801.06
+ ,3504.37
+ ,3032.6
+ ,3047.03
+ ,2962.34
+ ,2197.82
+ ,2014.45
+ ,1862.83
+ ,1905.41
+ ,1810.99
+ ,1670.07
+ ,1864.44
+ ,2052.02
+ ,2029.6
+ ,2070.83
+ ,2293.41
+ ,2443.27
+ ,2513.17
+ ,2466.92
+ ,2502.66)
> par10 = 'FALSE'
> par9 = '1'
> par8 = '2'
> par7 = '0'
> par6 = '3'
> par5 = '12'
> par4 = '0'
> 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: 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     sar1     sar2    sma1
      0.2901  -0.0609  0.1600  -0.3441  -0.1030  0.4440
s.e.  0.0696   0.0729  0.0765   0.3439   0.0993  0.3368

sigma^2 estimated as 13795:  log likelihood = -1255.98,  aic = 2525.96
> (forecast <- predict(arima.out,par1))
$pred
Time Series:
Start = 205 
End = 216 
Frequency = 1 
 [1] 1670.915 1600.046 1592.432 1544.229 1501.054 1458.237 1419.605 1477.944
 [9] 1461.667 1374.981 1407.535 1393.058

$se
Time Series:
Start = 205 
End = 216 
Frequency = 1 
 [1] 117.4505 191.7124 246.0646 300.0820 351.0239 396.5574 438.5897 477.8570
 [9] 514.4614 548.8277 581.3257 612.1784

> (lb <- forecast$pred - 1.96 * forecast$se)
Time Series:
Start = 205 
End = 216 
Frequency = 1 
 [1] 1440.7126 1224.2895 1110.1450  956.0683  813.0474  680.9846  559.9692
 [8]  541.3443  453.3230  299.2788  268.1364  193.1882
> (ub <- forecast$pred + 1.96 * forecast$se)
Time Series:
Start = 205 
End = 216 
Frequency = 1 
 [1] 1901.118 1975.802 2074.718 2132.390 2189.061 2235.490 2279.241 2414.544
 [9] 2470.012 2450.684 2546.933 2592.927
> 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] 1145.110 1176.860 1206.410 1192.720 1214.820 1199.070 1157.470 1100.100
  [9] 1095.630 1105.630 1137.790 1124.720 1152.600 1211.850 1239.620 1244.130
 [17] 1198.420 1227.990 1304.920 1340.260 1307.320 1356.510 1383.290 1437.870
 [25] 1494.560 1521.420 1498.760 1488.750 1524.620 1439.270 1423.110 1466.850
 [33] 1425.830 1363.450 1389.180 1395.890 1368.430 1349.030 1299.880 1365.410
 [41] 1451.040 1433.750 1464.650 1475.570 1471.160 1429.120 1452.460 1538.090
 [49] 1631.590 1665.500 1690.600 1711.740 1734.100 1748.090 1703.450 1745.740
 [57] 1751.010 1795.650 1852.130 1877.100 1989.310 2097.760 2154.870 2152.180
 [65] 2250.270 2346.900 2525.560 2409.360 2394.360 2401.330 2354.320 2450.410
 [73] 2504.670 2661.390 2880.400 3064.420 3141.120 3327.700 3564.950 3403.130
 [81] 3149.900 3006.840 3230.660 3361.130 3484.740 3411.130 3288.180 3280.370
 [89] 3173.950 3165.260 3092.710 3053.050 3181.960 2999.930 3249.570 3210.520
 [97] 3030.290 2803.470 2767.630 2882.600 2863.360 2897.060 3012.610 3142.950
[105] 3032.930 3045.780 3110.520 3013.240 2987.100 2995.550 2833.180 2848.960
[113] 2794.830 2845.260 2915.020 2892.630 2604.420 2641.650 2659.810 2638.530
[121] 2720.250 2745.880 2735.700 2811.700 2799.430 2555.280 2304.980 2214.950
[129] 2065.810 1940.490 2042.000 1995.370 1946.810 1765.900 1635.250 1833.420
[137] 1910.430 1959.670 1969.600 2061.410 2093.480 2120.880 2174.560 2196.720
[145] 2350.440 2440.250 2408.640 2472.810 2407.600 2454.620 2448.050 2497.840
[153] 2645.640 2756.760 2849.270 2921.440 2981.850 3080.580 3106.220 3119.310
[161] 3061.260 3097.310 3161.690 3257.160 3277.010 3295.320 3363.990 3494.170
[169] 3667.030 3813.060 3917.960 3895.510 3801.060 3570.120 3701.610 3862.270
[177] 3970.100 4138.520 4199.750 4290.890 4443.910 4502.640 4356.980 4591.270
[185] 4696.960 4621.400 4562.840 4202.520 4296.490 4435.230 4105.180 4116.680
[193] 3844.490 3720.980 3674.400 3857.620 3801.060 3504.370 3032.600 3047.030
[201] 2962.340 2197.820 2014.450 1862.830 1670.915 1600.046 1592.432 1544.229
[209] 1501.054 1458.237 1419.605 1477.944 1461.667 1374.981 1407.535 1393.058
> (perc.se <- (ub-forecast$pred)/1.96/forecast$pred)
Time Series:
Start = 205 
End = 216 
Frequency = 1 
 [1] 0.0702911 0.1198168 0.1545213 0.1943248 0.2338516 0.2719430 0.3089519
 [8] 0.3233255 0.3519689 0.3991529 0.4130098 0.4394494
> postscript(file="/var/www/html/rcomp/tmp/1a0xq1292348219.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/2oadz1292348219.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/3o2rd1292348219.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/4938j1292348219.tab") 
> 
> try(system("convert tmp/1a0xq1292348219.ps tmp/1a0xq1292348219.png",intern=TRUE))
character(0)
> try(system("convert tmp/2oadz1292348219.ps tmp/2oadz1292348219.png",intern=TRUE))
character(0)
> 
> 
> proc.time()
   user  system elapsed 
  2.006   0.492   4.456