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(1657
+ ,1418
+ ,1501
+ ,1315
+ ,1621
+ ,2308
+ ,3554
+ ,3318
+ ,3252
+ ,2921
+ ,2133
+ ,2040
+ ,1858
+ ,1833
+ ,2094
+ ,2173
+ ,2366
+ ,2074
+ ,2522
+ ,1822
+ ,1952
+ ,2232
+ ,1755
+ ,1791
+ ,2075
+ ,1850
+ ,2137
+ ,2467
+ ,2154
+ ,2289
+ ,2628
+ ,2074
+ ,2798
+ ,2194
+ ,2442
+ ,2565
+ ,2063
+ ,2070
+ ,2539
+ ,1898
+ ,2139
+ ,2408
+ ,2725
+ ,2201
+ ,2311
+ ,2548
+ ,2276
+ ,2351
+ ,2280
+ ,2057
+ ,2479
+ ,2379
+ ,2295
+ ,2456
+ ,2546
+ ,2844
+ ,2260
+ ,2981
+ ,2678
+ ,3440
+ ,2842
+ ,2450
+ ,2669
+ ,2570
+ ,2540
+ ,2318
+ ,2930
+ ,2947
+ ,2799
+ ,2695
+ ,2498
+ ,2260
+ ,2160
+ ,2058
+ ,2533
+ ,2150
+ ,2172
+ ,2155
+ ,3016
+ ,2333
+ ,2355
+ ,2825
+ ,2214
+ ,2360
+ ,2299
+ ,1746
+ ,2069
+ ,2267
+ ,1878
+ ,2266
+ ,2282
+ ,2085
+ ,2277
+ ,2251
+ ,1828
+ ,1954
+ ,1851
+ ,1570
+ ,1852
+ ,2187
+ ,1855
+ ,2218
+ ,2253
+ ,2028
+ ,2169
+ ,1997
+ ,2034
+ ,1791
+ ,1627
+ ,1631
+ ,2319
+ ,1707
+ ,1747
+ ,2397
+ ,2059
+ ,2251
+ ,2558
+ ,2406
+ ,2049
+ ,2074
+ ,1734
+ ,1983
+ ,2121
+ ,1905
+ ,2126
+ ,2363
+ ,2173
+ ,2710
+ ,2137
+ ,2742
+ ,2419
+ ,2194
+ ,2660
+ ,2189
+ ,2310
+ ,2349
+ ,2540
+ ,2434
+ ,2916
+ ,2446
+ ,2375
+ ,3032
+ ,2218
+ ,1920
+ ,2039
+ ,1889
+ ,2014
+ ,2105
+ ,2153
+ ,2309
+ ,2955
+ ,2225
+ ,2160
+ ,2386
+ ,1653
+ ,1099
+ ,5010
+ ,2672
+ ,2729
+ ,2955
+ ,2409
+ ,3086
+ ,3384
+ ,2458
+ ,2913
+ ,2448
+ ,2215
+ ,2179
+ ,2461
+ ,2098
+ ,2621
+ ,2703
+ ,2388
+ ,3880
+ ,3310
+ ,3093
+ ,3237
+ ,3002
+ ,2670
+ ,2311
+ ,2062
+ ,2059
+ ,2465
+ ,2213
+ ,2028
+ ,2322
+ ,2825
+ ,2687
+ ,2373
+ ,2889
+ ,2708
+ ,2542
+ ,2477
+ ,2419
+ ,2977
+ ,3001
+ ,3075
+ ,2870
+ ,3756
+ ,3443
+ ,2948
+ ,3560
+ ,3257
+ ,2600
+ ,2741
+ ,2349
+ ,2783
+ ,2845
+ ,2987
+ ,2696
+ ,3874
+ ,2912
+ ,2743
+ ,3857
+ ,2660
+ ,2226
+ ,2942
+ ,2420
+ ,2516
+ ,2421
+ ,2631
+ ,2887
+ ,3328
+ ,2587
+ ,2695
+ ,3669
+ ,2773
+ ,2527
+ ,2750
+ ,2014
+ ,2763
+ ,2726
+ ,1826
+ ,2713
+ ,3040
+ ,2405
+ ,2526
+ ,2526
+ ,2529
+ ,2474
+ ,2576
+ ,2219
+ ,2900
+ ,2274
+ ,2184
+ ,2629
+ ,2739
+ ,2933
+ ,3144
+ ,3354
+ ,3357
+ ,3329)
> par10 = 'FALSE'
> par9 = '1'
> par8 = '0'
> par7 = '1'
> par6 = '2'
> par5 = '12'
> par4 = '1'
> par3 = '0'
> par2 = '0.0'
> 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 <- 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      ma1     sma1
      0.8872  0.0139  -0.5572  -0.9067
s.e.  0.2351  0.1614   0.2250   0.0739

sigma^2 estimated as 0.02227:  log likelihood = 99.93,  aic = -189.85
> (forecast <- predict(arima.out,par1))
$pred
Time Series:
Start = 241 
End = 252 
Frequency = 1 
 [1] 7.827438 7.675638 7.829063 7.802157 7.752088 7.878774 8.015195 7.865426
 [9] 7.859519 7.956593 7.792339 7.708338

$se
Time Series:
Start = 241 
End = 252 
Frequency = 1 
 [1] 0.1495198 0.1574429 0.1639790 0.1691122 0.1731815 0.1764268 0.1790265
 [8] 0.1811158 0.1827990 0.1841571 0.1852540 0.1861405

> (lb <- forecast$pred - 1.96 * forecast$se)
Time Series:
Start = 241 
End = 252 
Frequency = 1 
 [1] 7.534379 7.367050 7.507665 7.470697 7.412653 7.532977 7.664304 7.510438
 [9] 7.501233 7.595645 7.429242 7.343503
> (ub <- forecast$pred + 1.96 * forecast$se)
Time Series:
Start = 241 
End = 252 
Frequency = 1 
 [1] 8.120497 7.984226 8.150462 8.133617 8.091524 8.224570 8.366087 8.220413
 [9] 8.217805 8.317541 8.155437 8.073173
> 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] 1657.000 1418.000 1501.000 1315.000 1621.000 2308.000 3554.000 3318.000
  [9] 3252.000 2921.000 2133.000 2040.000 1858.000 1833.000 2094.000 2173.000
 [17] 2366.000 2074.000 2522.000 1822.000 1952.000 2232.000 1755.000 1791.000
 [25] 2075.000 1850.000 2137.000 2467.000 2154.000 2289.000 2628.000 2074.000
 [33] 2798.000 2194.000 2442.000 2565.000 2063.000 2070.000 2539.000 1898.000
 [41] 2139.000 2408.000 2725.000 2201.000 2311.000 2548.000 2276.000 2351.000
 [49] 2280.000 2057.000 2479.000 2379.000 2295.000 2456.000 2546.000 2844.000
 [57] 2260.000 2981.000 2678.000 3440.000 2842.000 2450.000 2669.000 2570.000
 [65] 2540.000 2318.000 2930.000 2947.000 2799.000 2695.000 2498.000 2260.000
 [73] 2160.000 2058.000 2533.000 2150.000 2172.000 2155.000 3016.000 2333.000
 [81] 2355.000 2825.000 2214.000 2360.000 2299.000 1746.000 2069.000 2267.000
 [89] 1878.000 2266.000 2282.000 2085.000 2277.000 2251.000 1828.000 1954.000
 [97] 1851.000 1570.000 1852.000 2187.000 1855.000 2218.000 2253.000 2028.000
[105] 2169.000 1997.000 2034.000 1791.000 1627.000 1631.000 2319.000 1707.000
[113] 1747.000 2397.000 2059.000 2251.000 2558.000 2406.000 2049.000 2074.000
[121] 1734.000 1983.000 2121.000 1905.000 2126.000 2363.000 2173.000 2710.000
[129] 2137.000 2742.000 2419.000 2194.000 2660.000 2189.000 2310.000 2349.000
[137] 2540.000 2434.000 2916.000 2446.000 2375.000 3032.000 2218.000 1920.000
[145] 2039.000 1889.000 2014.000 2105.000 2153.000 2309.000 2955.000 2225.000
[153] 2160.000 2386.000 1653.000 1099.000 5010.000 2672.000 2729.000 2955.000
[161] 2409.000 3086.000 3384.000 2458.000 2913.000 2448.000 2215.000 2179.000
[169] 2461.000 2098.000 2621.000 2703.000 2388.000 3880.000 3310.000 3093.000
[177] 3237.000 3002.000 2670.000 2311.000 2062.000 2059.000 2465.000 2213.000
[185] 2028.000 2322.000 2825.000 2687.000 2373.000 2889.000 2708.000 2542.000
[193] 2477.000 2419.000 2977.000 3001.000 3075.000 2870.000 3756.000 3443.000
[201] 2948.000 3560.000 3257.000 2600.000 2741.000 2349.000 2783.000 2845.000
[209] 2987.000 2696.000 3874.000 2912.000 2743.000 3857.000 2660.000 2226.000
[217] 2942.000 2420.000 2516.000 2421.000 2631.000 2887.000 3328.000 2587.000
[225] 2695.000 3669.000 2773.000 2527.000 2750.000 2014.000 2763.000 2726.000
[233] 1826.000 2713.000 3040.000 2405.000 2526.000 2526.000 2529.000 2474.000
[241] 2508.494 2155.198 2512.575 2445.872 2326.425 2640.632 3026.601 2605.619
[249] 2590.274 2854.332 2421.977 2226.838
> (perc.se <- (ub-forecast$pred)/1.96/forecast$pred)
Time Series:
Start = 241 
End = 252 
Frequency = 1 
 [1] 0.1737355 0.1844395 0.1933956 0.2005103 0.2062015 0.2107730 0.2144560
 [8] 0.2174297 0.2198341 0.2217800 0.2233554 0.2246310
> postscript(file="/var/wessaorg/rcomp/tmp/1tzgj1323202433.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/28fk11323202433.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/3tmvf1323202433.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/42giq1323202433.tab") 
> 
> try(system("convert tmp/1tzgj1323202433.ps tmp/1tzgj1323202433.png",intern=TRUE))
character(0)
> try(system("convert tmp/28fk11323202433.ps tmp/28fk11323202433.png",intern=TRUE))
character(0)
> 
> 
> proc.time()
   user  system elapsed 
  2.009   0.145   2.145