R version 2.9.0 (2009-04-17) Copyright (C) 2009 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. 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(87.28 + ,87.28 + ,87.09 + ,86.92 + ,87.59 + ,90.72 + ,90.69 + ,90.3 + ,89.55 + ,88.94 + ,88.41 + ,87.82 + ,87.07 + ,86.82 + ,86.4 + ,86.02 + ,85.66 + ,85.32 + ,85 + ,84.67 + ,83.94 + ,82.83 + ,81.95 + ,81.19 + ,80.48 + ,78.86 + ,69.47 + ,68.77 + ,70.06 + ,73.95 + ,75.8 + ,77.79 + ,81.57 + ,83.07 + ,84.34 + ,85.1 + ,85.25 + ,84.26 + ,83.63 + ,86.44 + ,85.3 + ,84.1 + ,83.36 + ,82.48 + ,81.58 + ,80.47 + ,79.34 + ,82.13 + ,81.69 + ,80.7 + ,79.88 + ,79.16 + ,78.38 + ,77.42 + ,76.47 + ,75.46 + ,74.48 + ,78.27 + ,80.7 + ,79.91 + ,78.75 + ,77.78 + ,81.14 + ,81.08 + ,80.03 + ,78.91 + ,78.01 + ,76.9 + ,75.97 + ,81.93 + ,80.27 + ,78.67 + ,77.42 + ,76.16 + ,74.7 + ,76.39 + ,76.04 + ,74.65 + ,73.29 + ,71.79 + ,74.39 + ,74.91 + ,74.54 + ,73.08 + ,72.75 + ,71.32 + ,70.38 + ,70.35 + ,70.01 + ,69.36 + ,67.77 + ,69.26 + ,69.8 + ,68.38 + ,67.62 + ,68.39 + ,66.95 + ,65.21 + ,66.64 + ,63.45 + ,60.66 + ,62.34 + ,60.32 + ,58.64 + ,60.46 + ,58.59 + ,61.87 + ,61.85 + ,67.44 + ,77.06 + ,91.74 + ,93.15 + ,94.15 + ,93.11 + ,91.51 + ,89.96 + ,88.16 + ,86.98 + ,88.03 + ,86.24 + ,84.65 + ,83.23 + ,81.7 + ,80.25 + ,78.8 + ,77.51 + ,76.2 + ,75.04 + ,74 + ,75.49 + ,77.14 + ,76.15 + ,76.27 + ,78.19 + ,76.49 + ,77.31 + ,76.65 + ,74.99 + ,73.51 + ,72.07 + ,70.59 + ,71.96 + ,76.29 + ,74.86 + ,74.93 + ,71.9 + ,71.01 + ,77.47 + ,75.78 + ,76.6 + ,76.07 + ,74.57 + ,73.02 + ,72.65 + ,73.16 + ,71.53 + ,69.78 + ,67.98 + ,69.96 + ,72.16 + ,70.47 + ,68.86 + ,67.37 + ,65.87 + ,72.16 + ,71.34 + ,69.93 + ,68.44 + ,67.16 + ,66.01 + ,67.25 + ,70.91 + ,69.75 + ,68.59 + ,67.48 + ,66.31 + ,64.81 + ,66.58 + ,65.97 + ,64.7 + ,64.7 + ,60.94 + ,59.08 + ,58.42 + ,57.77 + ,57.11 + ,53.31 + ,49.96 + ,49.4 + ,48.84 + ,48.3 + ,47.74 + ,47.24 + ,46.76 + ,46.29 + ,48.9 + ,49.23 + ,48.53 + ,48.03 + ,54.34 + ,53.79 + ,53.24 + ,52.96 + ,52.17 + ,51.7 + ,58.55 + ,78.2 + ,77.03 + ,76.19 + ,77.15 + ,75.87 + ,95.47 + ,109.67 + ,112.28 + ,112.01 + ,107.93 + ,105.96 + ,105.06 + ,102.98 + ,102.2 + ,105.23 + ,101.85 + ,99.89 + ,96.23 + ,94.76 + ,91.51 + ,91.63 + ,91.54 + ,85.23 + ,87.83 + ,87.38 + ,84.44 + ,85.19 + ,84.03 + ,86.73 + ,102.52 + ,104.45 + ,106.98 + ,107.02 + ,99.26 + ,94.45 + ,113.44 + ,157.33 + ,147.38 + ,171.89 + ,171.95 + ,132.71 + ,126.02 + ,121.18 + ,115.45 + ,110.48 + ,117.85 + ,117.63 + ,124.65 + ,109.59 + ,111.27 + ,99.78 + ,98.21 + ,99.2 + ,97.97 + ,89.55 + ,87.91 + ,93.34 + ,94.42 + ,93.2 + ,90.29 + ,91.46 + ,89.98 + ,88.35 + ,88.41 + ,82.44 + ,79.89 + ,75.69 + ,75.66 + ,84.5 + ,96.73 + ,87.48 + ,82.39 + ,83.48 + ,79.31 + ,78.16 + ,72.77 + ,72.45 + ,68.46 + ,67.62 + ,68.76 + ,70.07 + ,68.55 + ,65.3 + ,58.96 + ,59.17 + ,62.37 + ,66.28 + ,55.62 + ,55.23 + ,55.85 + ,56.75 + ,50.89 + ,53.88 + ,52.95 + ,55.08 + ,53.61 + ,58.78 + ,61.85 + ,55.91 + ,53.32 + ,46.41 + ,44.57 + ,50 + ,50 + ,53.36 + ,46.23 + ,50.45 + ,49.07 + ,45.85 + ,48.45 + ,49.96 + ,46.53 + ,50.51 + ,47.58 + ,48.05 + ,46.84 + ,47.67 + ,49.16 + ,55.54 + ,55.82 + ,58.22 + ,56.19 + ,57.77 + ,63.19 + ,54.76 + ,55.74 + ,62.54 + ,61.39 + ,69.6 + ,79.23 + ,80 + ,93.68 + ,107.63 + ,100.18 + ,97.3 + ,90.45 + ,80.64 + ,80.58 + ,75.82 + ,85.59 + ,89.35 + ,89.42 + ,104.73 + ,95.32 + ,89.27 + ,90.44 + ,86.97 + ,79.98 + ,81.22 + ,87.35 + ,83.64 + ,82.22) > par10 = 'FALSE' > par9 = '0' > par8 = '0' > par7 = '0' > par6 = '0' > par5 = '12' > par4 = '0' > par3 = '1' > par2 = '1' > par1 = '24' > par1 <- as.numeric(par1) #cut off periods > par1 <- 28 > 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 > par6 <- 3 > par7 <- as.numeric(par7) #q > par7 <- 3 > 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 ar2 ar3 ma1 ma2 ma3 -0.1056 -0.4009 -0.6392 0.2425 0.4914 0.8027 s.e. 0.1484 0.1044 0.1388 0.1254 0.0768 0.1219 sigma^2 estimated as 24.51: log likelihood = -993.95, aic = 2001.9 > (forecast <- predict(arima.out,par1)) $pred Time Series: Start = 331 End = 358 Frequency = 1 [1] 63.82758 65.27299 65.19499 64.21627 63.42702 63.95259 64.83905 65.03921 [9] 64.32678 63.75517 63.97320 64.63469 64.84279 64.41627 63.95508 64.04176 [17] 64.49011 64.70279 64.44519 64.10056 64.10429 64.40670 64.59355 64.45020 [25] 64.19714 64.16190 64.35870 64.51379 $se Time Series: Start = 331 End = 358 Frequency = 1 [1] 4.950619 7.495610 9.604162 11.598525 12.997599 14.110672 15.133533 [8] 16.289984 17.436100 18.435791 19.261599 20.057364 20.907781 21.786601 [15] 22.589505 23.297264 23.971990 24.678114 25.404590 26.094582 26.723492 [22] 27.325378 27.941311 28.572160 29.183545 29.755683 30.305310 30.859958 > (lb <- forecast$pred - 1.96 * forecast$se) Time Series: Start = 331 End = 358 Frequency = 1 [1] 54.124363 50.581592 46.370837 41.483163 37.951729 36.295669 35.177326 [8] 33.110843 30.152021 27.621019 26.220466 25.322260 23.863537 21.714532 [15] 19.679650 18.379126 17.505011 16.333688 14.652195 12.955180 11.726240 [22] 10.848959 9.828580 8.448769 6.997392 5.840764 4.960288 4.028273 > (ub <- forecast$pred + 1.96 * forecast$se) Time Series: Start = 331 End = 358 Frequency = 1 [1] 73.53079 79.96439 84.01915 86.94938 88.90232 91.60950 94.50077 [8] 96.96758 98.50153 99.88932 101.72593 103.94713 105.82204 107.11801 [15] 108.23051 109.70440 111.47521 113.07190 114.23819 115.24594 116.48233 [22] 117.96444 119.35852 120.45164 121.39689 122.48304 123.75710 124.99931 > 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] 87.28000 87.28000 87.09000 86.92000 87.59000 90.72000 90.69000 [8] 90.30000 89.55000 88.94000 88.41000 87.82000 87.07000 86.82000 [15] 86.40000 86.02000 85.66000 85.32000 85.00000 84.67000 83.94000 [22] 82.83000 81.95000 81.19000 80.48000 78.86000 69.47000 68.77000 [29] 70.06000 73.95000 75.80000 77.79000 81.57000 83.07000 84.34000 [36] 85.10000 85.25000 84.26000 83.63000 86.44000 85.30000 84.10000 [43] 83.36000 82.48000 81.58000 80.47000 79.34000 82.13000 81.69000 [50] 80.70000 79.88000 79.16000 78.38000 77.42000 76.47000 75.46000 [57] 74.48000 78.27000 80.70000 79.91000 78.75000 77.78000 81.14000 [64] 81.08000 80.03000 78.91000 78.01000 76.90000 75.97000 81.93000 [71] 80.27000 78.67000 77.42000 76.16000 74.70000 76.39000 76.04000 [78] 74.65000 73.29000 71.79000 74.39000 74.91000 74.54000 73.08000 [85] 72.75000 71.32000 70.38000 70.35000 70.01000 69.36000 67.77000 [92] 69.26000 69.80000 68.38000 67.62000 68.39000 66.95000 65.21000 [99] 66.64000 63.45000 60.66000 62.34000 60.32000 58.64000 60.46000 [106] 58.59000 61.87000 61.85000 67.44000 77.06000 91.74000 93.15000 [113] 94.15000 93.11000 91.51000 89.96000 88.16000 86.98000 88.03000 [120] 86.24000 84.65000 83.23000 81.70000 80.25000 78.80000 77.51000 [127] 76.20000 75.04000 74.00000 75.49000 77.14000 76.15000 76.27000 [134] 78.19000 76.49000 77.31000 76.65000 74.99000 73.51000 72.07000 [141] 70.59000 71.96000 76.29000 74.86000 74.93000 71.90000 71.01000 [148] 77.47000 75.78000 76.60000 76.07000 74.57000 73.02000 72.65000 [155] 73.16000 71.53000 69.78000 67.98000 69.96000 72.16000 70.47000 [162] 68.86000 67.37000 65.87000 72.16000 71.34000 69.93000 68.44000 [169] 67.16000 66.01000 67.25000 70.91000 69.75000 68.59000 67.48000 [176] 66.31000 64.81000 66.58000 65.97000 64.70000 64.70000 60.94000 [183] 59.08000 58.42000 57.77000 57.11000 53.31000 49.96000 49.40000 [190] 48.84000 48.30000 47.74000 47.24000 46.76000 46.29000 48.90000 [197] 49.23000 48.53000 48.03000 54.34000 53.79000 53.24000 52.96000 [204] 52.17000 51.70000 58.55000 78.20000 77.03000 76.19000 77.15000 [211] 75.87000 95.47000 109.67000 112.28000 112.01000 107.93000 105.96000 [218] 105.06000 102.98000 102.20000 105.23000 101.85000 99.89000 96.23000 [225] 94.76000 91.51000 91.63000 91.54000 85.23000 87.83000 87.38000 [232] 84.44000 85.19000 84.03000 86.73000 102.52000 104.45000 106.98000 [239] 107.02000 99.26000 94.45000 113.44000 157.33000 147.38000 171.89000 [246] 171.95000 132.71000 126.02000 121.18000 115.45000 110.48000 117.85000 [253] 117.63000 124.65000 109.59000 111.27000 99.78000 98.21000 99.20000 [260] 97.97000 89.55000 87.91000 93.34000 94.42000 93.20000 90.29000 [267] 91.46000 89.98000 88.35000 88.41000 82.44000 79.89000 75.69000 [274] 75.66000 84.50000 96.73000 87.48000 82.39000 83.48000 79.31000 [281] 78.16000 72.77000 72.45000 68.46000 67.62000 68.76000 70.07000 [288] 68.55000 65.30000 58.96000 59.17000 62.37000 66.28000 55.62000 [295] 55.23000 55.85000 56.75000 50.89000 53.88000 52.95000 55.08000 [302] 53.61000 58.78000 61.85000 55.91000 53.32000 46.41000 44.57000 [309] 50.00000 50.00000 53.36000 46.23000 50.45000 49.07000 45.85000 [316] 48.45000 49.96000 46.53000 50.51000 47.58000 48.05000 46.84000 [323] 47.67000 49.16000 55.54000 55.82000 58.22000 56.19000 57.77000 [330] 63.19000 63.82758 65.27299 65.19499 64.21627 63.42702 63.95259 [337] 64.83905 65.03921 64.32678 63.75517 63.97320 64.63469 64.84279 [344] 64.41627 63.95508 64.04176 64.49011 64.70279 64.44519 64.10056 [351] 64.10429 64.40670 64.59355 64.45020 64.19714 64.16190 64.35870 [358] 64.51379 > (perc.se <- (ub-forecast$pred)/1.96/forecast$pred) Time Series: Start = 331 End = 358 Frequency = 1 [1] 0.07756239 0.11483480 0.14731440 0.18061660 0.20492211 0.22064271 [7] 0.23340152 0.25046405 0.27105509 0.28916543 0.30108856 0.31031885 [13] 0.32243803 0.33821581 0.35320893 0.36378236 0.37171574 0.38140726 [19] 0.39420459 0.40708820 0.41687529 0.42426297 0.43257123 0.44332149 [25] 0.45459260 0.46375936 0.47088136 0.47834669 > postscript(file="/var/www/html/rcomp/tmp/1lmfk1260395143.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() 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/rcomp/tmp/29zel1260395143.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() null device 1 > > #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
(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/rcomp/tmp/3jfqu1260395143.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/rcomp/tmp/4fuhb1260395143.tab") > > system("convert tmp/1lmfk1260395143.ps tmp/1lmfk1260395143.png") > system("convert tmp/29zel1260395143.ps tmp/29zel1260395143.png") > > > proc.time() user system elapsed 1.048 0.357 1.178