R version 2.8.0 (2008-10-20)
Copyright (C) 2008 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.

  Natural language support but running in an English locale

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(100.21
+ ,100.36
+ ,100.62
+ ,100.78
+ ,100.93
+ ,100.70
+ ,100.00
+ ,100.20
+ ,99.68
+ ,99.56
+ ,100.06
+ ,100.50
+ ,99.30
+ ,99.37
+ ,99.20
+ ,98.11
+ ,97.60
+ ,97.76
+ ,98.06
+ ,98.25
+ ,98.50
+ ,97.39
+ ,98.09
+ ,97.78
+ ,98.12
+ ,97.50
+ ,97.30
+ ,97.64
+ ,96.88
+ ,97.40
+ ,98.27
+ ,97.94
+ ,98.61
+ ,98.72
+ ,98.62
+ ,98.56
+ ,98.06
+ ,97.40
+ ,97.76
+ ,97.05
+ ,97.85
+ ,97.40
+ ,97.27
+ ,97.93
+ ,98.60
+ ,98.70
+ ,98.88
+ ,98.27
+ ,97.85
+ ,97.70
+ ,96.97
+ ,97.72
+ ,97.66
+ ,99.00
+ ,98.86
+ ,99.56
+ ,100.19
+ ,100.37
+ ,100.01
+ ,99.68
+ ,99.78
+ ,99.36
+ ,99.21
+ ,99.26
+ ,99.26
+ ,100.43
+ ,101.50
+ ,102.27
+ ,102.69
+ ,103.47
+ ,104.02
+ ,103.55
+ ,103.77
+ ,104.19
+ ,103.64
+ ,103.68
+ ,105.39
+ ,106.61
+ ,108.12
+ ,109.22
+ ,110.17
+ ,110.31
+ ,111.06
+ ,111.14
+ ,111.39
+ ,112.51
+ ,111.28
+ ,112.22
+ ,113.19
+ ,114.32
+ ,115.34
+ ,116.61
+ ,117.83
+ ,117.70
+ ,118.51
+ ,118.82
+ ,119.49
+ ,119.57
+ ,120.00
+ ,121.96
+ ,121.45
+ ,123.41
+ ,124.44
+ ,126.25
+ ,127.41
+ ,127.63
+ ,129.19
+ ,129.82
+ ,130.45
+ ,132.02
+ ,132.72
+ ,132.96
+ ,135.06
+ ,137.04
+ ,137.83
+ ,139.17
+ ,140.35
+ ,141.01
+ ,141.89
+ ,143.28
+ ,142.90
+ ,143.37
+ ,145.03
+ ,146.05
+ ,147.39
+ ,149.58
+ ,151.02
+ ,153.57
+ ,155.60
+ ,157.18
+ ,158.77
+ ,159.95
+ ,161.34
+ ,161.95
+ ,163.36
+ ,165.00
+ ,166.65
+ ,168.65
+ ,170.29
+ ,172.70
+ ,173.79
+ ,176.45
+ ,177.58
+ ,179.19
+ ,181.01
+ ,184.08
+ ,185.63
+ ,188.51
+ ,190.18
+ ,192.19
+ ,193.47
+ ,196.73
+ ,200.39
+ ,203.24
+ ,205.53
+ ,208.21
+ ,208.88
+ ,212.85
+ ,216.41
+ ,216.23
+ ,219.27
+ ,222.02
+ ,224.89
+ ,230.37
+ ,232.29
+ ,235.53
+ ,236.92
+ ,242.37
+ ,242.75
+ ,244.19
+ ,247.94
+ ,248.80
+ ,250.18
+ ,251.55
+ ,254.40
+ ,255.72
+ ,257.69
+ ,258.37
+ ,258.22
+ ,258.59
+ ,257.45
+ ,257.45
+ ,256.73
+ ,258.82
+ ,257.99
+ ,262.85
+ ,262.58
+ ,261.55
+ ,261.25
+ ,259.78
+ ,256.26
+ ,254.29
+ ,248.50
+ ,241.88
+ ,238.53
+ ,232.24
+ ,232.46
+ ,225.79
+ ,221.63
+ ,219.62
+ ,215.94
+ ,211.81
+ ,205.57
+ ,201.25
+ ,194.70
+ ,187.94
+ ,185.61
+ ,181.15
+ ,186.50
+ ,183.21
+ ,182.61
+ ,187.09
+ ,189.10
+ ,191.25
+ ,190.74
+ ,190.79)
> par10 = 'FALSE'
> par9 = '1'
> par8 = '1'
> par7 = '1'
> par6 = '3'
> par5 = '12'
> par4 = '0'
> par3 = '2'
> par2 = '-1.0'
> par1 = '0'
> par1 <- 33
> 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      ma1    sar1     sma1
      -0.1866  0.0586  0.1212  -0.9233  0.7954  -0.5815
s.e.   0.0896  0.0891  0.0807   0.0336  0.1067   0.1522

sigma^2 estimated as 1.982e-09:  log likelihood = 1554.76,  aic = -3095.52
> (forecast <- predict(arima.out,par1))
$pred
Time Series:
Start = 184 
End = 216 
Frequency = 1 
 [1] 0.003882153 0.003869044 0.003851445 0.003829355 0.003801381 0.003781857
 [7] 0.003765020 0.003757319 0.003739014 0.003740516 0.003729020 0.003717390
[13] 0.003710283 0.003698253 0.003683864 0.003668244 0.003646994 0.003632926
[19] 0.003621136 0.003616499 0.003603513 0.003606275 0.003598691 0.003591013
[25] 0.003586928 0.003578928 0.003569054 0.003558199 0.003542867 0.003533246
[31] 0.003525439 0.003523320 0.003514560

$se
Time Series:
Start = 184 
End = 216 
Frequency = 1 
 [1] 4.452233e-05 5.960655e-05 7.564149e-05 9.282209e-05 1.082414e-04
 [6] 1.237630e-04 1.391500e-04 1.543666e-04 1.696536e-04 1.849935e-04
[11] 2.004335e-04 2.160156e-04 2.353686e-04 2.541230e-04 2.733445e-04
[16] 2.930440e-04 3.128076e-04 3.328457e-04 3.531118e-04 3.735896e-04
[21] 3.943103e-04 4.152677e-04 4.364669e-04 4.579131e-04 4.819132e-04
[26] 5.057613e-04 5.300973e-04 5.549223e-04 5.799578e-04 6.053293e-04
[31] 6.309984e-04 6.569468e-04 6.831897e-04

> (lb <- forecast$pred - 1.96 * forecast$se)
Time Series:
Start = 184 
End = 216 
Frequency = 1 
 [1] 0.003794890 0.003752215 0.003703187 0.003647424 0.003589228 0.003539281
 [7] 0.003492286 0.003454760 0.003406493 0.003377929 0.003336170 0.003293999
[13] 0.003248960 0.003200172 0.003148109 0.003093878 0.003033891 0.002980548
[19] 0.002929037 0.002884263 0.002830665 0.002792350 0.002743216 0.002693503
[25] 0.002642378 0.002587636 0.002530063 0.002470551 0.002406149 0.002346801
[31] 0.002288682 0.002235704 0.002175508
> (ub <- forecast$pred + 1.96 * forecast$se)
Time Series:
Start = 184 
End = 216 
Frequency = 1 
 [1] 0.003969417 0.003985873 0.003999702 0.004011286 0.004013534 0.004024432
 [7] 0.004037753 0.004059877 0.004071535 0.004103103 0.004121869 0.004140781
[13] 0.004171605 0.004196334 0.004219619 0.004242610 0.004260097 0.004285303
[19] 0.004313235 0.004348734 0.004376361 0.004420200 0.004454166 0.004488522
[25] 0.004531478 0.004570221 0.004608044 0.004645847 0.004679584 0.004719692
[31] 0.004762195 0.004810935 0.004853612
> 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] 100.2100 100.3600 100.6200 100.7800 100.9300 100.7000 100.0000 100.2000
  [9]  99.6800  99.5600 100.0600 100.5000  99.3000  99.3700  99.2000  98.1100
 [17]  97.6000  97.7600  98.0600  98.2500  98.5000  97.3900  98.0900  97.7800
 [25]  98.1200  97.5000  97.3000  97.6400  96.8800  97.4000  98.2700  97.9400
 [33]  98.6100  98.7200  98.6200  98.5600  98.0600  97.4000  97.7600  97.0500
 [41]  97.8500  97.4000  97.2700  97.9300  98.6000  98.7000  98.8800  98.2700
 [49]  97.8500  97.7000  96.9700  97.7200  97.6600  99.0000  98.8600  99.5600
 [57] 100.1900 100.3700 100.0100  99.6800  99.7800  99.3600  99.2100  99.2600
 [65]  99.2600 100.4300 101.5000 102.2700 102.6900 103.4700 104.0200 103.5500
 [73] 103.7700 104.1900 103.6400 103.6800 105.3900 106.6100 108.1200 109.2200
 [81] 110.1700 110.3100 111.0600 111.1400 111.3900 112.5100 111.2800 112.2200
 [89] 113.1900 114.3200 115.3400 116.6100 117.8300 117.7000 118.5100 118.8200
 [97] 119.4900 119.5700 120.0000 121.9600 121.4500 123.4100 124.4400 126.2500
[105] 127.4100 127.6300 129.1900 129.8200 130.4500 132.0200 132.7200 132.9600
[113] 135.0600 137.0400 137.8300 139.1700 140.3500 141.0100 141.8900 143.2800
[121] 142.9000 143.3700 145.0300 146.0500 147.3900 149.5800 151.0200 153.5700
[129] 155.6000 157.1800 158.7700 159.9500 161.3400 161.9500 163.3600 165.0000
[137] 166.6500 168.6500 170.2900 172.7000 173.7900 176.4500 177.5800 179.1900
[145] 181.0100 184.0800 185.6300 188.5100 190.1800 192.1900 193.4700 196.7300
[153] 200.3900 203.2400 205.5300 208.2100 208.8800 212.8500 216.4100 216.2300
[161] 219.2700 222.0200 224.8900 230.3700 232.2900 235.5300 236.9200 242.3700
[169] 242.7500 244.1900 247.9400 248.8000 250.1800 251.5500 254.4000 255.7200
[177] 257.6900 258.3700 258.2200 258.5900 257.4500 257.4500 256.7300 257.5890
[185] 258.4618 259.6428 261.1406 263.0623 264.4204 265.6029 266.1472 267.4502
[193] 267.3428 268.1670 269.0059 269.5212 270.3980 271.4541 272.6100 274.1984
[201] 275.2603 276.1564 276.5105 277.5070 277.2944 277.8788 278.4730 278.7901
[209] 279.4132 280.1863 281.0410 282.2573 283.0258 283.6526 283.8232 284.5306
> (perc.se <- (ub-forecast$pred)/1.96/forecast$pred)
Time Series:
Start = 184 
End = 216 
Frequency = 1 
 [1] 0.01173218 0.01588569 0.02042605 0.02544867 0.03015729 0.03496841
 [7] 0.03984495 0.04468228 0.04980300 0.05476537 0.06007893 0.06557852
[13] 0.07244428 0.07940916 0.08682815 0.09471735 0.10310443 0.11167266
[19] 0.12055562 0.12952689 0.13929957 0.14871620 0.15910772 0.17000649
[25] 0.18237859 0.19545301 0.20951944 0.22461476 0.24103150 0.25793807
[31] 0.27570386 0.29384341 0.31403681
> postscript(file="/var/www/html/freestat/rcomp/tmp/1kadt1292683708.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/freestat/rcomp/tmp/2g1tj1292683708.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/freestat/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/freestat/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/freestat/rcomp/tmp/38m9t1292683709.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/freestat/rcomp/tmp/4cm8h1292683709.tab") 
> 
> try(system("convert tmp/1kadt1292683708.ps tmp/1kadt1292683708.png",intern=TRUE))
character(0)
> try(system("convert tmp/2g1tj1292683708.ps tmp/2g1tj1292683708.png",intern=TRUE))
character(0)
> 
> 
> proc.time()
   user  system elapsed 
  2.205   0.444   2.303