x <- c(235.1
,280.7
,264.6
,240.7
,201.4
,240.8
,241.1
,223.8
,206.1
,174.7
,203.3
,220.5
,299.5
,347.4
,338.3
,327.7
,351.6
,396.6
,438.8
,395.6
,363.5
,378.8
,357
,369
,464.8
,479.1
,431.3
,366.5
,326.3
,355.1
,331.6
,261.3
,249
,205.5
,235.6
,240.9
,264.9
,253.8
,232.3
,193.8
,177
,213.2
,207.2
,180.6
,188.6
,175.4
,199
,179.6
,225.8
,234
,200.2
,183.6
,178.2
,203.2
,208.5
,191.8
,172.8
,148
,159.4
,154.5
,213.2
,196.4
,182.8
,176.4
,153.6
,173.2
,171
,151.2
,161.9
,157.2
,201.7
,236.4
,356.1
,398.3
,403.7
,384.6
,365.8
,368.1
,367.9
,347
,343.3
,292.9
,311.5
,300.9
,366.9
,356.9
,329.7
,316.2
,269
,289.3
,266.2
,253.6
,233.8
,228.4
,253.6
,260.1
,306.6
,309.2
,309.5
,271
,279.9
,317.9
,298.4
,246.7
,227.3
,209.1
,259.9
,266
,320.6
,308.5
,282.2
,262.7
,263.5
,313.1
,284.3
,252.6
,250.3
,246.5
,312.7
,333.2
,446.4
,511.6
,515.5
,506.4
,483.2
,522.3
,509.8
,460.7
,405.8
,375
,378.5
,406.8
,467.8
,469.8
,429.8
,355.8
,332.7
,378
,360.5
,334.7
,319.5
,323.1
,363.6
,352.1
,411.9
,388.6
,416.4
,360.7
,338
,417.2
,388.4
,371.1
,331.5
,353.7
,396.7
,447
,533.5
,565.4
,542.3
,488.7
,467.1
,531.3
,496.1
,444
,403.4
,386.3
,394.1
,404.1
,462.1
,448.1
,432.3
,386.3
,395.2
,421.9
,382.9
,384.2
,345.5
,323.4
,372.6
,376
,462.7
,487
,444.2
,399.3
,394.9
,455.4
,414
,375.5
,347
,339.4
,385.8
,378.8
,451.8
,446.1
,422.5
,383.1
,352.8
,445.3
,367.5
,355.1
,326.2
,319.8
,331.8
,340.9
,394.1
,417.2
,369.9
,349.2
,321.4
,405.7
,342.9
,316.5
,284.2
,270.9
,288.8
,278.8
,324.4
,310.9
,299
,273
,279.3
,359.2
,305
,282.1
,250.3
,246.5
,257.9
,266.5
,315.9
,318.4
,295.4
,266.4
,245.8
,362.8
,324.9
,294.2
,289.5
,295.2
,290.3
,272
,307.4
,328.7
,292.9
,249.1
,230.4
,361.5
,321.7
,277.2
,260.7
,251
,257.6
,241.8
,287.5
,292.3
,274.7
,254.2
,230
,339
,318.2
,287
,295.8
,284
,271
,262.7
,340.6
,379.4
,373.3
,355.2
,338.4
,466.9
,451
,422
,429.2
,425.9
,460.7
,463.6
,541.4
,544.2
,517.5
,469.4
,439.4
,549
,533
,506.1
,484
,457
,481.5
,469.5
,544.7
,541.2
,521.5
,469.7
,434.4
,542.6
,517.3
,485.7
,465.8
,447
,426.6
,411.6
,467.5
,484.5
,451.2
,417.4
,379.9
,484.7
,455
,420.8
,416.5
,376.3
,405.6
,405.8
,500.8
,514
,475.5
,430.1
,414.4
,538
,526
,488.5
,520.2
,504.4
,568.5
,610.6
,818
,830.9
,835.9
,782
,762.3
,856.9
,820.9
,769.6
,752.2
,724.4
,723.1
,719.5
,817.4
,803.3
,752.5
,689
,630.4
,765.5
,757.7
,732.2
,702.6
,683.3
,709.5
,702.2
,784.8
,810.9
,755.6
,656.8
,615.1
,745.3
,694.1
,675.7
,643.7
,622.1
,634.6
,588
,689.7
,673.9
,647.9
,568.8
,545.7)
par10 = 'FALSE'
par9 = '0'
par8 = '0'
par7 = '0'
par6 = '0'
par5 = '12'
par4 = '0'
par3 = '0'
par2 = '1'
par1 = '12'
par10 <- 'FALSE'
par9 <- '0'
par8 <- '0'
par7 <- '0'
par6 <- '0'
par5 <- '12'
par4 <- '0'
par3 <- '0'
par2 <- '1'
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'))
(forecast <- predict(arima.out,par1))
(lb <- forecast$pred - 1.96 * forecast$se)
(ub <- forecast$pred + 1.96 * forecast$se)
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))
(perc.se <- (ub-forecast$pred)/1.96/forecast$pred)
postscript(file="/var/fisher/rcomp/tmp/19mx71352745233.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()
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/fisher/rcomp/tmp/2jeal1352745233.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()

#Note: the /var/fisher/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
load(file="/var/fisher/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/fisher/rcomp/tmp/3qizt1352745233.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/fisher/rcomp/tmp/44dlx1352745233.tab") 

try(system("convert tmp/19mx71352745233.ps tmp/19mx71352745233.png",intern=TRUE))
try(system("convert tmp/2jeal1352745233.ps tmp/2jeal1352745233.png",intern=TRUE))