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(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
+ ,632.6
+ ,643.8
+ ,593.1
+ ,579.7
+ ,546
+ ,562.9
+ ,572.5)
> par10 = 'FALSE'
> par9 = '1'
> par8 = '0'
> par7 = '1'
> par6 = '2'
> par5 = '12'
> par4 = '1'
> par3 = '1'
> par2 = '0.5'
> 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     ar2      ma1     sma1
      0.4963  0.1698  -0.3959  -0.7181
s.e.  0.1716  0.0684   0.1711   0.0417

sigma^2 estimated as 0.2869:  log likelihood = -280.17,  aic = 570.35
> (forecast <- predict(arima.out,par1))
$pred
Time Series:
Start = 361 
End = 372 
Frequency = 1 
 [1] 26.10305 26.08018 25.24227 23.87661 23.04841 25.46023 24.82645 24.17909
 [9] 23.82369 23.35966 23.74394 23.47967

$se
Time Series:
Start = 361 
End = 372 
Frequency = 1 
 [1] 0.5356629 0.7964708 1.0650526 1.3169505 1.5556715 1.7805271 1.9923956
 [8] 2.1922062 2.3810184 2.5598723 2.7297500 2.8915500

> (lb <- forecast$pred - 1.96 * forecast$se)
Time Series:
Start = 361 
End = 372 
Frequency = 1 
 [1] 25.05315 24.51910 23.15477 21.29538 19.99929 21.97039 20.92136 19.88236
 [9] 19.15690 18.34231 18.39363 17.81223
> (ub <- forecast$pred + 1.96 * forecast$se)
Time Series:
Start = 361 
End = 372 
Frequency = 1 
 [1] 27.15295 27.64127 27.32977 26.45783 26.09753 28.95006 28.73155 28.47581
 [9] 28.49049 28.37701 29.09425 29.14711
> 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] 235.1000 280.7000 264.6000 240.7000 201.4000 240.8000 241.1000 223.8000
  [9] 206.1000 174.7000 203.3000 220.5000 299.5000 347.4000 338.3000 327.7000
 [17] 351.6000 396.6000 438.8000 395.6000 363.5000 378.8000 357.0000 369.0000
 [25] 464.8000 479.1000 431.3000 366.5000 326.3000 355.1000 331.6000 261.3000
 [33] 249.0000 205.5000 235.6000 240.9000 264.9000 253.8000 232.3000 193.8000
 [41] 177.0000 213.2000 207.2000 180.6000 188.6000 175.4000 199.0000 179.6000
 [49] 225.8000 234.0000 200.2000 183.6000 178.2000 203.2000 208.5000 191.8000
 [57] 172.8000 148.0000 159.4000 154.5000 213.2000 196.4000 182.8000 176.4000
 [65] 153.6000 173.2000 171.0000 151.2000 161.9000 157.2000 201.7000 236.4000
 [73] 356.1000 398.3000 403.7000 384.6000 365.8000 368.1000 367.9000 347.0000
 [81] 343.3000 292.9000 311.5000 300.9000 366.9000 356.9000 329.7000 316.2000
 [89] 269.0000 289.3000 266.2000 253.6000 233.8000 228.4000 253.6000 260.1000
 [97] 306.6000 309.2000 309.5000 271.0000 279.9000 317.9000 298.4000 246.7000
[105] 227.3000 209.1000 259.9000 266.0000 320.6000 308.5000 282.2000 262.7000
[113] 263.5000 313.1000 284.3000 252.6000 250.3000 246.5000 312.7000 333.2000
[121] 446.4000 511.6000 515.5000 506.4000 483.2000 522.3000 509.8000 460.7000
[129] 405.8000 375.0000 378.5000 406.8000 467.8000 469.8000 429.8000 355.8000
[137] 332.7000 378.0000 360.5000 334.7000 319.5000 323.1000 363.6000 352.1000
[145] 411.9000 388.6000 416.4000 360.7000 338.0000 417.2000 388.4000 371.1000
[153] 331.5000 353.7000 396.7000 447.0000 533.5000 565.4000 542.3000 488.7000
[161] 467.1000 531.3000 496.1000 444.0000 403.4000 386.3000 394.1000 404.1000
[169] 462.1000 448.1000 432.3000 386.3000 395.2000 421.9000 382.9000 384.2000
[177] 345.5000 323.4000 372.6000 376.0000 462.7000 487.0000 444.2000 399.3000
[185] 394.9000 455.4000 414.0000 375.5000 347.0000 339.4000 385.8000 378.8000
[193] 451.8000 446.1000 422.5000 383.1000 352.8000 445.3000 367.5000 355.1000
[201] 326.2000 319.8000 331.8000 340.9000 394.1000 417.2000 369.9000 349.2000
[209] 321.4000 405.7000 342.9000 316.5000 284.2000 270.9000 288.8000 278.8000
[217] 324.4000 310.9000 299.0000 273.0000 279.3000 359.2000 305.0000 282.1000
[225] 250.3000 246.5000 257.9000 266.5000 315.9000 318.4000 295.4000 266.4000
[233] 245.8000 362.8000 324.9000 294.2000 289.5000 295.2000 290.3000 272.0000
[241] 307.4000 328.7000 292.9000 249.1000 230.4000 361.5000 321.7000 277.2000
[249] 260.7000 251.0000 257.6000 241.8000 287.5000 292.3000 274.7000 254.2000
[257] 230.0000 339.0000 318.2000 287.0000 295.8000 284.0000 271.0000 262.7000
[265] 340.6000 379.4000 373.3000 355.2000 338.4000 466.9000 451.0000 422.0000
[273] 429.2000 425.9000 460.7000 463.6000 541.4000 544.2000 517.5000 469.4000
[281] 439.4000 549.0000 533.0000 506.1000 484.0000 457.0000 481.5000 469.5000
[289] 544.7000 541.2000 521.5000 469.7000 434.4000 542.6000 517.3000 485.7000
[297] 465.8000 447.0000 426.6000 411.6000 467.5000 484.5000 451.2000 417.4000
[305] 379.9000 484.7000 455.0000 420.8000 416.5000 376.3000 405.6000 405.8000
[313] 500.8000 514.0000 475.5000 430.1000 414.4000 538.0000 526.0000 488.5000
[321] 520.2000 504.4000 568.5000 610.6000 818.0000 830.9000 835.9000 782.0000
[329] 762.3000 856.9000 820.9000 769.6000 752.2000 724.4000 723.1000 719.5000
[337] 817.4000 803.3000 752.5000 689.0000 630.4000 765.5000 757.7000 732.2000
[345] 702.6000 683.3000 709.5000 702.2000 784.8000 810.9000 755.6000 656.8000
[353] 615.1000 745.3000 694.1000 675.7000 643.7000 622.1000 634.6000 588.0000
[361] 681.3690 680.1760 637.1722 570.0923 531.2292 648.2231 616.3527 584.6282
[369] 567.5684 545.6735 563.7748 551.2949
> (perc.se <- (ub-forecast$pred)/1.96/forecast$pred)
Time Series:
Start = 361 
End = 372 
Frequency = 1 
 [1] 0.04186756 0.06290661 0.08787576 0.11627583 0.14392079 0.14945315
 [7] 0.17312931 0.19744239 0.21946436 0.24270785 0.25583809 0.27602822
> postscript(file="/var/www/html/freestat/rcomp/tmp/1sk5p1291063597.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/2g32i1291063597.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/3n40u1291063597.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/495y01291063597.tab") 
> 
> try(system("convert tmp/1sk5p1291063597.ps tmp/1sk5p1291063597.png",intern=TRUE))
character(0)
> try(system("convert tmp/2g32i1291063597.ps tmp/2g32i1291063597.png",intern=TRUE))
character(0)
> 
> 
> proc.time()
   user  system elapsed 
  2.708   0.471   5.447