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(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 = '1'
> 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.2057
s.e.  0.0518

sigma^2 estimated as 0.04782:  log likelihood = 36.02,  aic = -68.04
> (forecast <- predict(arima.out,par1))
$pred
Time Series:
Start = 358 
End = 369 
Frequency = 1 
 [1] 7.329161 7.328140 7.327930 7.327887 7.327878 7.327877 7.327876 7.327876
 [9] 7.327876 7.327876 7.327876 7.327876

$se
Time Series:
Start = 358 
End = 369 
Frequency = 1 
 [1] 0.2186722 0.3425418 0.4379712 0.5170549 0.5857373 0.6472056 0.7033285
 [8] 0.7552939 0.8039074 0.8497443 0.8932322 0.9346989

> (lb <- forecast$pred - 1.96 * forecast$se)
Time Series:
Start = 358 
End = 369 
Frequency = 1 
 [1] 6.900563 6.656759 6.469507 6.314460 6.179833 6.059354 5.949352 5.847500
 [9] 5.752218 5.662377 5.577141 5.495866
> (ub <- forecast$pred + 1.96 * forecast$se)
Time Series:
Start = 358 
End = 369 
Frequency = 1 
 [1] 7.757758 7.999522 8.186354 8.341315 8.475924 8.596399 8.706400 8.808252
 [9] 8.903535 8.993375 9.078611 9.159886
> 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.4238 145.3732 145.3628
[361] 145.3607 145.3602 145.3601 145.3601 145.3601 145.3601 145.3601 145.3601
[369] 145.3601
> (perc.se <- (ub-forecast$pred)/1.96/forecast$pred)
Time Series:
Start = 358 
End = 369 
Frequency = 1 
 [1] 0.07789285 0.12500925 0.16279873 0.19511114 0.22391391 0.25027984
 [7] 0.27484204 0.29800360 0.32003803 0.34114052 0.36145603 0.38109570
> postscript(file="/var/www/html/freestat/rcomp/tmp/1ycuo1292267308.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/2c4sf1292267308.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/3bxam1292267309.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/4wgqa1292267309.tab") 
> 
> try(system("convert tmp/1ycuo1292267308.ps tmp/1ycuo1292267308.png",intern=TRUE))
character(0)
> try(system("convert tmp/2c4sf1292267308.ps tmp/2c4sf1292267308.png",intern=TRUE))
character(0)
> 
> 
> proc.time()
   user  system elapsed 
  0.929   0.450   1.043