x <- c(1.3866 ,1.3582 ,1.3332 ,1.3595 ,1.3617 ,1.3684 ,1.3394 ,1.3262 ,1.3173 ,1.3085 ,1.327 ,1.3182 ,1.293 ,1.291 ,1.2984 ,1.2795 ,1.299 ,1.3174 ,1.326 ,1.3111 ,1.2816 ,1.276 ,1.2849 ,1.2818 ,1.2829 ,1.2796 ,1.3008 ,1.2967 ,1.2938 ,1.2833 ,1.2823 ,1.2765 ,1.2634 ,1.2596 ,1.2705 ,1.2591 ,1.2798 ,1.2763 ,1.2795 ,1.2782 ,1.2644 ,1.2596 ,1.2615 ,1.2555 ,1.2555 ,1.2658 ,1.2565 ,1.2783 ,1.2786 ,1.2782 ,1.2905 ,1.3042 ,1.2942 ,1.313 ,1.3671 ,1.3549 ,1.3558 ,1.3507 ,1.3494 ,1.3607 ,1.3295 ,1.3193 ,1.3308 ,1.3246 ,1.3392 ,1.3425 ,1.3496 ,1.3255 ,1.3231 ,1.3273 ,1.3276 ,1.3173 ,1.3196 ,1.3058 ,1.2966 ,1.2932 ,1.2947 ,1.305 ,1.3232 ,1.3125 ,1.2992 ,1.3266 ,1.3275 ,1.3223 ,1.3403 ,1.3322 ,1.3363 ,1.3425 ,1.3574 ,1.3683 ,1.3623 ,1.3563 ,1.3518 ,1.3494 ,1.3612 ,1.369 ,1.3771 ,1.3972 ,1.401 ,1.3908 ,1.3901 ,1.3856 ,1.4098 ,1.422 ,1.4238 ,1.4207 ,1.4095 ,1.4177 ,1.3866 ,1.3959 ,1.4102 ,1.3969 ,1.4004 ,1.385 ,1.389 ,1.384 ,1.392 ,1.3932 ,1.3858 ,1.3978 ,1.4029 ,1.394 ,1.4096 ,1.4058 ,1.4134 ,1.4096 ,1.4049 ,1.4009 ,1.3897 ,1.4019 ,1.3901 ,1.399 ,1.3901 ,1.3975 ,1.3991 ,1.4089 ,1.413 ,1.409 ,1.4217 ,1.4223 ,1.4191 ,1.4229 ,1.4227 ,1.4269 ,1.4229 ,1.4104 ,1.4053 ,1.4138 ,1.4303 ,1.4384 ,1.441 ,1.437 ,1.4357 ,1.4202 ,1.4166 ,1.417 ,1.4293 ,1.4294 ,1.4072 ,1.4101 ,1.4112 ,1.4243 ,1.433 ,1.4323 ,1.4324 ,1.427 ,1.4268 ,1.4364 ,1.4272 ,1.4314 ,1.422 ,1.4335 ,1.4262 ,1.433 ,1.4473 ,1.4522 ,1.4545 ,1.4594 ,1.4561 ,1.4611 ,1.4671 ,1.4712 ,1.4705 ,1.4658 ,1.478 ,1.4783 ,1.4768 ,1.467 ,1.465 ,1.4549 ,1.4643 ,1.4539 ,1.4537 ,1.4616 ,1.4722 ,1.4694 ,1.4763 ,1.475 ,1.4765 ,1.4864 ,1.4881 ,1.4864 ,1.4869 ,1.4918 ,1.4971 ,1.4921 ,1.5 ,1.502 ,1.5019 ,1.4874 ,1.4785 ,1.4788 ,1.48 ,1.4772 ,1.4658 ,1.4761 ,1.4867 ,1.4862 ,1.4984 ,1.4966 ,1.5037 ,1.4922 ,1.4868 ,1.4965 ,1.4875 ,1.4957 ,1.4863 ,1.4815 ,1.4968 ,1.4969 ,1.5083 ,1.5071 ,1.4918 ,1.5023 ,1.5074 ,1.509 ,1.512 ,1.5068 ,1.4787 ,1.4774 ,1.4768 ,1.473 ,1.4757 ,1.4647 ,1.4541 ,1.456 ,1.4343 ,1.4337 ,1.4368 ,1.4279 ,1.4276 ,1.4398 ,1.4405 ,1.4433 ,1.4338 ,1.4406 ,1.4389 ,1.4442 ,1.435 ,1.4304 ,1.4273 ,1.4528 ,1.4481 ,1.4563 ,1.4486 ,1.4374 ,1.4369 ,1.4279 ,1.4132 ,1.4064 ,1.4135 ,1.4151 ,1.4085 ,1.4072 ,1.3999 ,1.3966 ,1.3913 ,1.3937 ,1.3984 ,1.3847 ,1.3691 ,1.3675 ,1.376 ,1.374 ,1.3718 ,1.3572 ,1.3607 ,1.3649 ,1.3726 ,1.3567 ,1.3519 ,1.3626 ,1.3577 ,1.3547 ,1.3489 ,1.357 ,1.3525 ,1.3548 ,1.3641 ,1.3668 ,1.3582 ,1.3662 ,1.3557 ,1.361 ,1.3657 ,1.3765 ,1.3705 ,1.3723 ,1.3756 ,1.366 ,1.3548 ,1.3471 ,1.3519 ,1.3338 ,1.3356 ,1.3353 ,1.3471 ,1.3482 ,1.3479 ,1.3468 ,1.3396 ,1.334 ,1.3296 ,1.3384 ,1.3585 ,1.3583 ,1.3615 ,1.3544 ,1.3535 ,1.3432 ,1.3486 ,1.3373 ,1.3339 ,1.3311 ,1.3321 ,1.329 ,1.3245 ,1.3256 ,1.3315 ,1.3238 ,1.3089 ,1.2924 ,1.2727 ,1.2746 ,1.2969 ,1.2698 ,1.2686 ,1.2587 ,1.2492 ,1.2349 ,1.2428 ,1.227 ,1.2334 ,1.2497 ,1.236 ,1.2223 ,1.2309 ,1.2255 ,1.2384 ,1.2307 ,1.2155 ,1.2218 ,1.2268 ,1.206 ,1.1959 ,1.1942) par10 = 'FALSE' par9 = '0' par8 = '1' par7 = '0' par6 = '0' par5 = '1' par4 = '0' par3 = '1' par2 = '1.9' par1 = '15' par1 <- 30 #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 <- 5 #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/freestat/rcomp/tmp/13x9f1292887059.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/www/html/freestat/rcomp/tmp/2szsj1292887060.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/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
(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/3gi7d1292887060.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/4kinj1292887060.tab") try(system("convert tmp/13x9f1292887059.ps tmp/13x9f1292887059.png",intern=TRUE)) try(system("convert tmp/2szsj1292887060.ps tmp/2szsj1292887060.png",intern=TRUE))