R version 2.12.0 (2010-10-15)
Copyright (C) 2010 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(168.67
+ ,164.83
+ ,184.38
+ ,180.81
+ ,190.54
+ ,181.41
+ ,155.67
+ ,135.99
+ ,125.88
+ ,126.09
+ ,114.86
+ ,127.98
+ ,127.98
+ ,125.11
+ ,125.93
+ ,128.2
+ ,125.93
+ ,111.94
+ ,120.01
+ ,124.09
+ ,126.02
+ ,136.41
+ ,143.79
+ ,141.67
+ ,143.9
+ ,155
+ ,144.83
+ ,141.4
+ ,137
+ ,141.02
+ ,131.11
+ ,132.83
+ ,136.73
+ ,141.18
+ ,137.86
+ ,133.79
+ ,128.53
+ ,125.87
+ ,124.27
+ ,123.96
+ ,128.15
+ ,126.4
+ ,127.86
+ ,129.31
+ ,132.56
+ ,141.28
+ ,145.55
+ ,146.54
+ ,143.14
+ ,145.72
+ ,148.21
+ ,150.4
+ ,149.94
+ ,146.66
+ ,143.37
+ ,145.29
+ ,140.24
+ ,136.12
+ ,140.25
+ ,140.64
+ ,145.58
+ ,143.73
+ ,141.27
+ ,140.66
+ ,141.94
+ ,141.16
+ ,134.31
+ ,132.93
+ ,133.07
+ ,140.48
+ ,154.85
+ ,196.77
+ ,235.3
+ ,226.52
+ ,237.62
+ ,224.07
+ ,208.74
+ ,174.54
+ ,170.63
+ ,172.23
+ ,198.36
+ ,175.91
+ ,154.63
+ ,134.31
+ ,121.75
+ ,119.6
+ ,102.04
+ ,106.3
+ ,116.38
+ ,103.72
+ ,98.56
+ ,100.9
+ ,110
+ ,118.26
+ ,124.77
+ ,125.22
+ ,126.38
+ ,137.14
+ ,134.74
+ ,134.3
+ ,136.39
+ ,141.83
+ ,139.24
+ ,128.89
+ ,134.83
+ ,130.43
+ ,132.09
+ ,144.95
+ ,149.5
+ ,137.57
+ ,139.38
+ ,143.06
+ ,138.65
+ ,123.21
+ ,85.91
+ ,77.4
+ ,77.84
+ ,67.76
+ ,70.72
+ ,72.55
+ ,75.83
+ ,84.01
+ ,93.96
+ ,93.73
+ ,92.02
+ ,88.26
+ ,86.48
+ ,94.42
+ ,94.92
+ ,91.41
+ ,84.84
+ ,89.89
+ ,86.32
+ ,89.57
+ ,93.72
+ ,92.27
+ ,87.59
+ ,85.5
+ ,82.81
+ ,81.62
+ ,87.45
+ ,79.86
+ ,78.52
+ ,75.1
+ ,72.99
+ ,67.88
+ ,70.14
+ ,65.43
+ ,60.26
+ ,58.38
+ ,57.68
+ ,52.42
+ ,52.73
+ ,61.4
+ ,67.13
+ ,77.46
+ ,68.66
+ ,67.46
+ ,62.77
+ ,56.88
+ ,61.48
+ ,61.99
+ ,71.56
+ ,76.56
+ ,79.82
+ ,75.05
+ ,77.07
+ ,80
+ ,77.21
+ ,82.16
+ ,85.57
+ ,89.23
+ ,121.98
+ ,142.56
+ ,217.67
+ ,198.07
+ ,220.1
+ ,198.68
+ ,181.64
+ ,167.47
+ ,172.33
+ ,168.71
+ ,178.22
+ ,172.81
+ ,168.83
+ ,152.25
+ ,143.83
+ ,151.41
+ ,131.87
+ ,125.38
+ ,123.23
+ ,103.99
+ ,109.38
+ ,123.79
+ ,119.05
+ ,122.01
+ ,128.56
+ ,127.91
+ ,120.47
+ ,122.49
+ ,114.05
+ ,120.62
+ ,119.61
+ ,115.01
+ ,131.83
+ ,167.2
+ ,193.82
+ ,204.43
+ ,264.5
+ ,212.55
+ ,186.52
+ ,185.17
+ ,184.38
+ ,161.45
+ ,154.15
+ ,174.25
+ ,175.04
+ ,175.87
+ ,154.82
+ ,147.08
+ ,134.35
+ ,121.56
+ ,113.86
+ ,119.89
+ ,108.07
+ ,107.07
+ ,115.14
+ ,116.03
+ ,111.48
+ ,103.24
+ ,103.23
+ ,99.69
+ ,108.91
+ ,104.21
+ ,90.85
+ ,87.64
+ ,81.06
+ ,92.2
+ ,114.02
+ ,123.56
+ ,109.17
+ ,101.65
+ ,97.95
+ ,92.56
+ ,91.76
+ ,84.1
+ ,84.67
+ ,74.52
+ ,73.83
+ ,75.37
+ ,70.47
+ ,64.5
+ ,64.98
+ ,66.94
+ ,65.93
+ ,65.51
+ ,68.94
+ ,63.67
+ ,58.47
+ ,59.68
+ ,57.71
+ ,56.53
+ ,58.96
+ ,55.6
+ ,57.34
+ ,60.51
+ ,66.38
+ ,65.78
+ ,58.43
+ ,55.16
+ ,53.09
+ ,52.02
+ ,57.58
+ ,64.05
+ ,70.18
+ ,63.86
+ ,65.22
+ ,67.6
+ ,61.66
+ ,65.32
+ ,66.18
+ ,61.34
+ ,62.29
+ ,63.6
+ ,65.51
+ ,62.58
+ ,62.36
+ ,64.88
+ ,73.73
+ ,77.51
+ ,77.47
+ ,74.34
+ ,75.81
+ ,82.16
+ ,73.96
+ ,73.17
+ ,80.99
+ ,79.81
+ ,89.51
+ ,102.57
+ ,107.11
+ ,122.23
+ ,134.69
+ ,128.79
+ ,126.16
+ ,119.98
+ ,108.45
+ ,108.43
+ ,98.17
+ ,106.09
+ ,108.81
+ ,103.03
+ ,124.36
+ ,118.52
+ ,112.2
+ ,114.71
+ ,107.96
+ ,101.21
+ ,102.77
+ ,112.13
+ ,109.36
+ ,110.91
+ ,123.57
+ ,129.95
+ ,124.46
+ ,122.34
+ ,116.61
+ ,114.59
+ ,112.52
+ ,118.67
+ ,116.8
+ ,123.63
+ ,128.04
+ ,134.57
+ ,130.33
+ ,136.47
+ ,139.05
+ ,158.21
+ ,148.07
+ ,137.74
+ ,139.74
+ ,144.08
+ ,145.35
+ ,145.77
+ ,140.56
+ ,121.41
+ ,120.44
+ ,116.97
+ ,128.03
+ ,128.51
+ ,127.76
+ ,134.58
+ ,147.64
+ ,144.46
+ ,137.6
+ ,146.87
+ ,145.67
+ ,151.95
+ ,150.23
+ ,155.86
+ ,154.4
+ ,156.36
+ ,162.13
+ ,171.06
+ ,174.01
+ ,193.52
+ ,205.26
+ ,212.8
+ ,222.1)
> par10 = 'FALSE'
> par9 = '0'
> par8 = '0'
> par7 = '0'
> par6 = '1'
> par5 = '12'
> par4 = '0'
> par3 = '1'
> par2 = '-0.4'
> 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
      0.2225
s.e.  0.0516

sigma^2 estimated as 2.193e-05:  log likelihood = 1404.32,  aic = -2804.64
> (forecast <- predict(arima.out,par1))
$pred
Time Series:
Start = 358 
End = 369 
Frequency = 1 
 [1] 0.1364484 0.1364705 0.1364754 0.1364765 0.1364768 0.1364768 0.1364768
 [8] 0.1364768 0.1364768 0.1364768 0.1364768 0.1364768

$se
Time Series:
Start = 358 
End = 369 
Frequency = 1 
 [1] 0.004683437 0.007396986 0.009497617 0.011238830 0.012749742 0.014100755
 [7] 0.015333412 0.016474139 0.017540847 0.018546306 0.019499991 0.020409160

> (lb <- forecast$pred - 1.96 * forecast$se)
Time Series:
Start = 358 
End = 369 
Frequency = 1 
 [1] 0.12726885 0.12197241 0.11786009 0.11444841 0.11148726 0.10883933
 [7] 0.10642334 0.10418752 0.10209677 0.10012607 0.09825685 0.09647487
> (ub <- forecast$pred + 1.96 * forecast$se)
Time Series:
Start = 358 
End = 369 
Frequency = 1 
 [1] 0.1456279 0.1509686 0.1550908 0.1585046 0.1614663 0.1641143 0.1665303
 [8] 0.1687661 0.1708569 0.1728276 0.1746968 0.1764788
> 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] 168.6700 164.8300 184.3800 180.8100 190.5400 181.4100 155.6700 135.9900
  [9] 125.8800 126.0900 114.8600 127.9800 127.9800 125.1100 125.9300 128.2000
 [17] 125.9300 111.9400 120.0100 124.0900 126.0200 136.4100 143.7900 141.6700
 [25] 143.9000 155.0000 144.8300 141.4000 137.0000 141.0200 131.1100 132.8300
 [33] 136.7300 141.1800 137.8600 133.7900 128.5300 125.8700 124.2700 123.9600
 [41] 128.1500 126.4000 127.8600 129.3100 132.5600 141.2800 145.5500 146.5400
 [49] 143.1400 145.7200 148.2100 150.4000 149.9400 146.6600 143.3700 145.2900
 [57] 140.2400 136.1200 140.2500 140.6400 145.5800 143.7300 141.2700 140.6600
 [65] 141.9400 141.1600 134.3100 132.9300 133.0700 140.4800 154.8500 196.7700
 [73] 235.3000 226.5200 237.6200 224.0700 208.7400 174.5400 170.6300 172.2300
 [81] 198.3600 175.9100 154.6300 134.3100 121.7500 119.6000 102.0400 106.3000
 [89] 116.3800 103.7200  98.5600 100.9000 110.0000 118.2600 124.7700 125.2200
 [97] 126.3800 137.1400 134.7400 134.3000 136.3900 141.8300 139.2400 128.8900
[105] 134.8300 130.4300 132.0900 144.9500 149.5000 137.5700 139.3800 143.0600
[113] 138.6500 123.2100  85.9100  77.4000  77.8400  67.7600  70.7200  72.5500
[121]  75.8300  84.0100  93.9600  93.7300  92.0200  88.2600  86.4800  94.4200
[129]  94.9200  91.4100  84.8400  89.8900  86.3200  89.5700  93.7200  92.2700
[137]  87.5900  85.5000  82.8100  81.6200  87.4500  79.8600  78.5200  75.1000
[145]  72.9900  67.8800  70.1400  65.4300  60.2600  58.3800  57.6800  52.4200
[153]  52.7300  61.4000  67.1300  77.4600  68.6600  67.4600  62.7700  56.8800
[161]  61.4800  61.9900  71.5600  76.5600  79.8200  75.0500  77.0700  80.0000
[169]  77.2100  82.1600  85.5700  89.2300 121.9800 142.5600 217.6700 198.0700
[177] 220.1000 198.6800 181.6400 167.4700 172.3300 168.7100 178.2200 172.8100
[185] 168.8300 152.2500 143.8300 151.4100 131.8700 125.3800 123.2300 103.9900
[193] 109.3800 123.7900 119.0500 122.0100 128.5600 127.9100 120.4700 122.4900
[201] 114.0500 120.6200 119.6100 115.0100 131.8300 167.2000 193.8200 204.4300
[209] 264.5000 212.5500 186.5200 185.1700 184.3800 161.4500 154.1500 174.2500
[217] 175.0400 175.8700 154.8200 147.0800 134.3500 121.5600 113.8600 119.8900
[225] 108.0700 107.0700 115.1400 116.0300 111.4800 103.2400 103.2300  99.6900
[233] 108.9100 104.2100  90.8500  87.6400  81.0600  92.2000 114.0200 123.5600
[241] 109.1700 101.6500  97.9500  92.5600  91.7600  84.1000  84.6700  74.5200
[249]  73.8300  75.3700  70.4700  64.5000  64.9800  66.9400  65.9300  65.5100
[257]  68.9400  63.6700  58.4700  59.6800  57.7100  56.5300  58.9600  55.6000
[265]  57.3400  60.5100  66.3800  65.7800  58.4300  55.1600  53.0900  52.0200
[273]  57.5800  64.0500  70.1800  63.8600  65.2200  67.6000  61.6600  65.3200
[281]  66.1800  61.3400  62.2900  63.6000  65.5100  62.5800  62.3600  64.8800
[289]  73.7300  77.5100  77.4700  74.3400  75.8100  82.1600  73.9600  73.1700
[297]  80.9900  79.8100  89.5100 102.5700 107.1100 122.2300 134.6900 128.7900
[305] 126.1600 119.9800 108.4500 108.4300  98.1700 106.0900 108.8100 103.0300
[313] 124.3600 118.5200 112.2000 114.7100 107.9600 101.2100 102.7700 112.1300
[321] 109.3600 110.9100 123.5700 129.9500 124.4600 122.3400 116.6100 114.5900
[329] 112.5200 118.6700 116.8000 123.6300 128.0400 134.5700 130.3300 136.4700
[337] 139.0500 158.2100 148.0700 137.7400 139.7400 144.0800 145.3500 145.7700
[345] 140.5600 121.4100 120.4400 116.9700 128.0300 128.5100 127.7600 134.5800
[353] 147.6400 144.4600 137.6000 146.8700 145.6700 145.4049 145.3460 145.3329
[361] 145.3300 145.3293 145.3292 145.3291 145.3291 145.3291 145.3291 145.3291
[369] 145.3291
> (perc.se <- (ub-forecast$pred)/1.96/forecast$pred)
Time Series:
Start = 358 
End = 369 
Frequency = 1 
 [1] 0.09703489 0.16539165 0.22594048 0.28204975 0.33571249 0.38810633
 [7] 0.43996088 0.49175969 0.54384601 0.59647914 0.64986588 0.70417891
> postscript(file="/var/www/rcomp/tmp/1i4bz1291722826.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/rcomp/tmp/2wwr71291722826.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/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/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/rcomp/tmp/3v6nm1291722826.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/rcomp/tmp/4z63a1291722826.tab") 
> try(system("convert tmp/1i4bz1291722826.ps tmp/1i4bz1291722826.png",intern=TRUE))
character(0)
> try(system("convert tmp/2wwr71291722826.ps tmp/2wwr71291722826.png",intern=TRUE))
character(0)
> 
> 
> proc.time()
   user  system elapsed 
  0.660   0.500   1.139