x <- c(97.4,97,105.4,102.7,98.1,104.5,87.4,89.9,109.8,111.7,98.6,96.9,95.1,97,112.7,102.9,97.4,111.4,87.4,96.8,114.1,110.3,103.9,101.6,94.6,95.9,104.7,102.8,98.1,113.9,80.9,95.7,113.2,105.9,108.8,102.3,99,100.7,115.5,100.7,109.9,114.6,85.4,100.5,114.8,116.5,112.9,102,106,105.3,118.8,106.1,109.3,117.2,92.5,104.2,112.5,122.4,113.3) par20 = '' par19 = '' par18 = '' par17 = '' par16 = '' par15 = '' par14 = '' par13 = '' par12 = '' par11 = '' par10 = 'FALSE' par9 = '0' par8 = '0' par7 = '0' par6 = '3' par5 = '12' par4 = '1' par3 = '0' par2 = '1' par1 = '12' ylab = '' xlab = '' main = '' #'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/www/html/rcomp/tmp/19otu1260531435.ps",horizontal=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/www/html/rcomp/tmp/2oek51260531436.ps",horizontal=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/www/html/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab load(file="/var/www/html/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