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(396
+ ,297
+ ,559
+ ,967
+ ,270
+ ,143
+ ,1562
+ ,109
+ ,371
+ ,656
+ ,511
+ ,655
+ ,465
+ ,525
+ ,885
+ ,497
+ ,1436
+ ,612
+ ,865
+ ,385
+ ,567
+ ,639
+ ,963
+ ,398
+ ,410
+ ,966
+ ,801
+ ,892
+ ,513
+ ,469
+ ,683
+ ,643
+ ,535
+ ,625
+ ,264
+ ,992
+ ,238
+ ,818
+ ,937
+ ,70
+ ,507
+ ,260
+ ,503
+ ,927
+ ,1269
+ ,537
+ ,910
+ ,532
+ ,345
+ ,918
+ ,1635
+ ,330
+ ,557
+ ,1178
+ ,740
+ ,452
+ ,218
+ ,764
+ ,255
+ ,454
+ ,866
+ ,574
+ ,1276
+ ,379
+ ,825
+ ,798
+ ,663
+ ,1069
+ ,921
+ ,858
+ ,711
+ ,503
+ ,382
+ ,464
+ ,717
+ ,690
+ ,462
+ ,657
+ ,385
+ ,577
+ ,619
+ ,479
+ ,817
+ ,752
+ ,430
+ ,451
+ ,537
+ ,519
+ ,1000
+ ,637
+ ,465
+ ,437
+ ,711
+ ,299
+ ,248
+ ,1162
+ ,714
+ ,905
+ ,649
+ ,512
+ ,472
+ ,905
+ ,786
+ ,489
+ ,479
+ ,617
+ ,925
+ ,351
+ ,1144
+ ,669
+ ,707
+ ,458
+ ,214
+ ,599
+ ,572
+ ,897
+ ,819
+ ,720
+ ,273
+ ,508
+ ,506
+ ,451
+ ,699
+ ,407
+ ,465
+ ,245
+ ,370
+ ,316
+ ,603
+ ,154
+ ,229
+ ,577
+ ,192
+ ,617
+ ,411
+ ,975
+ ,146
+ ,705
+ ,184
+ ,200
+ ,274
+ ,502
+ ,382
+ ,964
+ ,537
+ ,438
+ ,369
+ ,417
+ ,276
+ ,514
+ ,822
+ ,389
+ ,466
+ ,1255
+ ,694
+ ,1024
+ ,400
+ ,397
+ ,350
+ ,719
+ ,1277
+ ,356
+ ,457
+ ,1402
+ ,600
+ ,480
+ ,595
+ ,436
+ ,230
+ ,651
+ ,1367
+ ,564
+ ,716
+ ,747
+ ,467
+ ,671
+ ,861
+ ,319
+ ,612
+ ,433
+ ,434
+ ,503
+ ,85
+ ,564
+ ,824
+ ,74
+ ,259
+ ,69
+ ,535
+ ,239
+ ,438
+ ,459
+ ,426
+ ,288
+ ,498
+ ,454
+ ,376
+ ,225
+ ,555
+ ,252
+ ,208
+ ,130
+ ,481
+ ,389
+ ,565
+ ,173
+ ,278
+ ,609
+ ,422
+ ,445
+ ,387
+ ,339
+ ,181
+ ,245
+ ,384
+ ,212
+ ,399
+ ,229
+ ,224
+ ,203
+ ,333
+ ,384
+ ,636
+ ,185
+ ,93
+ ,581
+ ,248
+ ,304
+ ,344
+ ,407
+ ,170
+ ,312
+ ,507
+ ,224
+ ,340
+ ,168
+ ,443
+ ,204
+ ,367
+ ,210
+ ,335
+ ,364
+ ,178
+ ,206
+ ,279
+ ,387
+ ,490
+ ,238
+ ,343
+ ,232
+ ,530
+ ,291
+ ,67
+ ,397
+ ,467
+ ,178
+ ,175
+ ,299
+ ,154
+ ,106
+ ,189
+ ,194
+ ,135
+ ,201
+ ,207
+ ,280
+ ,260
+ ,227
+ ,239
+ ,333
+ ,428
+ ,230
+ ,292
+ ,350
+ ,186
+ ,326
+ ,155
+ ,75
+ ,361
+ ,261
+ ,299
+ ,300
+ ,450
+ ,183
+ ,238
+ ,165
+ ,234
+ ,176
+ ,329)
> par10 = 'FALSE'
> par9 = '1'
> par8 = '0'
> par7 = '2'
> par6 = '0'
> par5 = '1'
> par4 = '0'
> par3 = '0'
> par2 = '-0.1'
> par1 = '12'
> 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 <- 7 #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     ma2    sma1
      0.6030  0.9214  0.8384
s.e.  0.0274  0.0230  0.0252

sigma^2 estimated as 0.01887:  log likelihood = 151.52,  aic = -295.04
> (forecast <- predict(arima.out,par1))
$pred
Time Series:
Start = 278 
End = 289 
Frequency = 1 
 [1] 0.48363370 0.36673740 0.27285807 0.21104512 0.27369902 0.29197881
 [7] 0.29087845 0.19937649 0.04188702 0.00000000 0.00000000 0.00000000

$se
Time Series:
Start = 278 
End = 289 
Frequency = 1 
 [1] 0.1373599 0.1604025 0.2043190 0.2043189 0.2043189 0.2043189 0.2043189
 [8] 0.2345386 0.2446041 0.2666265 0.2666265 0.2666265

> (lb <- forecast$pred - 1.96 * forecast$se)
Time Series:
Start = 278 
End = 289 
Frequency = 1 
 [1]  0.21440838  0.05234843 -0.12760707 -0.18942002 -0.12676611 -0.10848633
 [7] -0.10958669 -0.26031925 -0.43753695 -0.52258796 -0.52258796 -0.52258796
> (ub <- forecast$pred + 1.96 * forecast$se)
Time Series:
Start = 278 
End = 289 
Frequency = 1 
 [1] 0.7528590 0.6811264 0.6733232 0.6115103 0.6741642 0.6924439 0.6913436
 [8] 0.6590722 0.5213110 0.5225880 0.5225880 0.5225880
> 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] 3.960000e+02 2.970000e+02 5.590000e+02 9.670000e+02 2.700000e+02
  [6] 1.430000e+02 1.562000e+03 1.090000e+02 3.710000e+02 6.560000e+02
 [11] 5.110000e+02 6.550000e+02 4.650000e+02 5.250000e+02 8.850000e+02
 [16] 4.970000e+02 1.436000e+03 6.120000e+02 8.650000e+02 3.850000e+02
 [21] 5.670000e+02 6.390000e+02 9.630000e+02 3.980000e+02 4.100000e+02
 [26] 9.660000e+02 8.010000e+02 8.920000e+02 5.130000e+02 4.690000e+02
 [31] 6.830000e+02 6.430000e+02 5.350000e+02 6.250000e+02 2.640000e+02
 [36] 9.920000e+02 2.380000e+02 8.180000e+02 9.370000e+02 7.000000e+01
 [41] 5.070000e+02 2.600000e+02 5.030000e+02 9.270000e+02 1.269000e+03
 [46] 5.370000e+02 9.100000e+02 5.320000e+02 3.450000e+02 9.180000e+02
 [51] 1.635000e+03 3.300000e+02 5.570000e+02 1.178000e+03 7.400000e+02
 [56] 4.520000e+02 2.180000e+02 7.640000e+02 2.550000e+02 4.540000e+02
 [61] 8.660000e+02 5.740000e+02 1.276000e+03 3.790000e+02 8.250000e+02
 [66] 7.980000e+02 6.630000e+02 1.069000e+03 9.210000e+02 8.580000e+02
 [71] 7.110000e+02 5.030000e+02 3.820000e+02 4.640000e+02 7.170000e+02
 [76] 6.900000e+02 4.620000e+02 6.570000e+02 3.850000e+02 5.770000e+02
 [81] 6.190000e+02 4.790000e+02 8.170000e+02 7.520000e+02 4.300000e+02
 [86] 4.510000e+02 5.370000e+02 5.190000e+02 1.000000e+03 6.370000e+02
 [91] 4.650000e+02 4.370000e+02 7.110000e+02 2.990000e+02 2.480000e+02
 [96] 1.162000e+03 7.140000e+02 9.050000e+02 6.490000e+02 5.120000e+02
[101] 4.720000e+02 9.050000e+02 7.860000e+02 4.890000e+02 4.790000e+02
[106] 6.170000e+02 9.250000e+02 3.510000e+02 1.144000e+03 6.690000e+02
[111] 7.070000e+02 4.580000e+02 2.140000e+02 5.990000e+02 5.720000e+02
[116] 8.970000e+02 8.190000e+02 7.200000e+02 2.730000e+02 5.080000e+02
[121] 5.060000e+02 4.510000e+02 6.990000e+02 4.070000e+02 4.650000e+02
[126] 2.450000e+02 3.700000e+02 3.160000e+02 6.030000e+02 1.540000e+02
[131] 2.290000e+02 5.770000e+02 1.920000e+02 6.170000e+02 4.110000e+02
[136] 9.750000e+02 1.460000e+02 7.050000e+02 1.840000e+02 2.000000e+02
[141] 2.740000e+02 5.020000e+02 3.820000e+02 9.640000e+02 5.370000e+02
[146] 4.380000e+02 3.690000e+02 4.170000e+02 2.760000e+02 5.140000e+02
[151] 8.220000e+02 3.890000e+02 4.660000e+02 1.255000e+03 6.940000e+02
[156] 1.024000e+03 4.000000e+02 3.970000e+02 3.500000e+02 7.190000e+02
[161] 1.277000e+03 3.560000e+02 4.570000e+02 1.402000e+03 6.000000e+02
[166] 4.800000e+02 5.950000e+02 4.360000e+02 2.300000e+02 6.510000e+02
[171] 1.367000e+03 5.640000e+02 7.160000e+02 7.470000e+02 4.670000e+02
[176] 6.710000e+02 8.610000e+02 3.190000e+02 6.120000e+02 4.330000e+02
[181] 4.340000e+02 5.030000e+02 8.500000e+01 5.640000e+02 8.240000e+02
[186] 7.400000e+01 2.590000e+02 6.900000e+01 5.350000e+02 2.390000e+02
[191] 4.380000e+02 4.590000e+02 4.260000e+02 2.880000e+02 4.980000e+02
[196] 4.540000e+02 3.760000e+02 2.250000e+02 5.550000e+02 2.520000e+02
[201] 2.080000e+02 1.300000e+02 4.810000e+02 3.890000e+02 5.650000e+02
[206] 1.730000e+02 2.780000e+02 6.090000e+02 4.220000e+02 4.450000e+02
[211] 3.870000e+02 3.390000e+02 1.810000e+02 2.450000e+02 3.840000e+02
[216] 2.120000e+02 3.990000e+02 2.290000e+02 2.240000e+02 2.030000e+02
[221] 3.330000e+02 3.840000e+02 6.360000e+02 1.850000e+02 9.300000e+01
[226] 5.810000e+02 2.480000e+02 3.040000e+02 3.440000e+02 4.070000e+02
[231] 1.700000e+02 3.120000e+02 5.070000e+02 2.240000e+02 3.400000e+02
[236] 1.680000e+02 4.430000e+02 2.040000e+02 3.670000e+02 2.100000e+02
[241] 3.350000e+02 3.640000e+02 1.780000e+02 2.060000e+02 2.790000e+02
[246] 3.870000e+02 4.900000e+02 2.380000e+02 3.430000e+02 2.320000e+02
[251] 5.300000e+02 2.910000e+02 6.700000e+01 3.970000e+02 4.670000e+02
[256] 1.780000e+02 1.750000e+02 2.990000e+02 1.540000e+02 1.060000e+02
[261] 1.890000e+02 1.940000e+02 1.350000e+02 2.010000e+02 2.070000e+02
[266] 2.800000e+02 2.600000e+02 2.270000e+02 2.390000e+02 3.330000e+02
[271] 4.280000e+02 2.300000e+02 2.920000e+02 3.500000e+02 1.860000e+02
[276] 3.260000e+02 1.550000e+02 1.428349e+03 2.272208e+04 4.371515e+05
[281] 5.704885e+06 4.239040e+05 2.220684e+05 2.306134e+05 1.007536e+07
[286] 6.014587e+13          Inf          Inf          Inf
> (perc.se <- (ub-forecast$pred)/1.96/forecast$pred)
Time Series:
Start = 278 
End = 289 
Frequency = 1 
 [1]  1.739260e+03  1.452978e+08  1.018937e+03  9.937629e-01  1.122658e+03
 [6]  1.017385e+04  8.856384e+03 -4.747700e-01 -5.102041e-01           NaN
[11]           NaN           NaN
> postscript(file="/var/wessaorg/rcomp/tmp/1vtrz1324396089.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/2sjka1324396089.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/38dml1324396089.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/4mjnm1324396089.tab") 
> 
> try(system("convert tmp/1vtrz1324396089.ps tmp/1vtrz1324396089.png",intern=TRUE))
character(0)
> try(system("convert tmp/2sjka1324396089.ps tmp/2sjka1324396089.png",intern=TRUE))
character(0)
> 
> 
> proc.time()
   user  system elapsed 
  0.948   0.163   1.105