FreeStatistics.org

Free Statistics

of Irreproducible Research!

Description of Statistical Computation

Author's title

Author*Unverified author*
R Software ModuleRscript (source code is shown below)
Title produced by softwareR console
Date of computationFri, 04 Jun 2010 21:32:51 +0000
Cite this page as followsStatistical Computations at FreeStatistics.org, Office for Research Development and Education, URL https://freestatistics.org/blog/index.php?v=date/2010/Jun/04/t1275687169nhrdwje2h5ruocx.htm/, Retrieved Sat, 23 May 2026 23:43:55 +0000
Statistical Computations at FreeStatistics.org, Office for Research Development and Education, URL https://freestatistics.org/blog/index.php?pk=77498, Retrieved Sat, 23 May 2026 23:43:55 +0000
QR Codes:

Original text written by user:
IsPrivate?No (this computation is public)
User-defined keywordsbullwhip effect, master thesis
Estimated Impact310

Tree of Dependent Computations

Family? (F = Feedback message, R = changed R code, M = changed R Module, P = changed Parameters, D = changed Data)
-       [R console] [Bullwhip effect] [2010-06-04 21:32:51] [d41d8cd98f00b204e9800998ecf8427e] [Current]
Feedback Forum

Post a new message

Dataset

Tables (Output of Computation)

> {
+     library(forecast)
+     iteratie = 6
+     iteratie2 = 6
+     nInst <- 100
+     nObs <- 150
+     mu <- 0
+     sig <- iteratie
+     s <- .... [TRUNCATED] 
[1] 'sigma = 6   horizon = 6   i = 3'
[1] 'mean Ft =65.1998925787576'
[1] 'mean Yt =19.435612426345'
[1] 'mean et = 45.7642801524126'
[1] '%error = 2.35466108031559'
[1] 'sigma = 6   horizon = 6   i = 29'
[1] 'mean Ft =671.881856432923'
[1] 'mean Yt =629.439363132271'
[1] 'mean et = 42.4424933006518'
[1] '%error = 0.0674290420755476'
[1] 'sigma = 6   horizon = 6   i = 46'
[1] 'mean Ft =84.2553561721589'
[1] 'mean Yt =32.3419622977136'
[1] 'mean et = 51.9133938744453'
[1] '%error = 1.60514051054086'
[1] 'sigma = 6   horizon = 6   i = 93'
[1] 'mean Ft =-216.569826031084'
[1] 'mean Yt =-262.142699387163'
[1] 'mean et = 45.5728733560792'
[1] '%error = 0.173847577913173'
[1] 'sigma = 6   horizon = 6   i = 100'
[1] 'mean Ft =1044.57157982021'
[1] 'mean Yt =1000.30660606767'
[1] 'mean et = 44.2649737525372'
[1] '%error = 0.0442514060029538'

Figures (Output of Computation)

Input Parameters & R Code

Parameters (Session):
Parameters (R input):
R code (body of R function):
{
    library(forecast)
    iteratie = 6
    iteratie2 = 6
    nInst <- 100
    nObs <- 150
    mu <- 0
    sig <- iteratie
    s <- 12
    d <- 1
    D <- 1
    horizon1 <- iteratie2
    horizon2 <- 4
    horizon3 <- 5
    x <- array(NA, dim = c(nObs, nInst))
    for (i in 1:nInst) {
        dum <- rnorm(nObs, mu, sig)
        for (j in 1:d) {
            dum = cumsum(dum)
        }
        for (j in 1:D) {
            for (h in (s + 1):nObs) {
                dum[h] = dum[h] + dum[h - s]
            }
        }
        x[, i] = dum
    }
    numseries <- length(x[1, ])
    n <- length(x[, 1])
    nmh1 <- n - horizon1
    nmh1p1 <- nmh1 + 1
    nmh2 <- n - horizon2
    nmh2p1 <- nmh2 + 1
    nmh3 <- n - horizon3
    nmh3p1 <- nmh3 + 1
    nmh12 <- n - horizon1 - horizon2
    nmh12p1 <- nmh12 + 1
    nmh123 <- n - horizon1 - horizon2 - horizon3
    nmh123p1 <- nmh123 + 1
    nmh23 <- n - horizon2 - horizon3
    nmh23p1 <- nmh23 + 1
    rarr1 <- array(NA, dim = c(horizon1, numseries))
    farr1 <- array(NA, dim = c(horizon1, numseries))
    parr1 <- array(NA, dim = c(numseries, 8))
    rarr2 <- array(NA, dim = c(horizon2, numseries))
    farr2 <- array(NA, dim = c(horizon2, numseries))
    parr2 <- array(NA, dim = c(numseries, 8))
    rarr3 <- array(NA, dim = c(horizon3, numseries))
    farr3 <- array(NA, dim = c(horizon3, numseries))
    parr3 <- array(NA, dim = c(numseries, 8))
    colnames(parr1) = list("ME", "RMSE", "MAE", "MPE", "MAPE", 
        "MASE", "ACF1", "TheilU")
    colnames(parr2) = list("ME", "RMSE", "MAE", "MPE", "MAPE", 
        "MASE", "ACF1", "TheilU")
    colnames(parr3) = list("ME", "RMSE", "MAE", "MPE", "MAPE", 
        "MASE", "ACF1", "TheilU")
    for (j in 1:numseries) {
        fit <- auto.arima(ts(x[1:nmh123, j], freq = s), d = d, 
            D = D)
        forec <- forecast(fit, horizon1)
        for (i in 1:horizon1) rarr1[i, j] <- forec$residuals[i]
        forec$mean <- forec$mean + 0.5 * sd(rarr1[, j])
        forec$resid <- forec$resid - 0.5 * sd(rarr1[, j])
        for (i in 1:horizon1) farr1[i, j] <- forec$mean[i]
        parr1[j, ] <- accuracy(forec, x[nmh123p1:nmh23, j])
    }
    x1 <- x
    for (j in 1:numseries) {
        for (i in 1:horizon1) {
            x1[nmh123 + i, j] = farr1[i, j]
        }
    }
    for (j in 1:numseries) {
        fit <- auto.arima(ts(x1[1:nmh23, j], freq = s), d = d, 
            D = D)
        forec <- forecast(fit, horizon2)
        for (i in 1:horizon2) rarr2[i, j] <- forec$residuals[i]
        forec$mean <- forec$mean + 0.5 * sd(rarr2[, j])
        forec$resid <- forec$resid - 0.5 * sd(rarr2[, j])
        for (i in 1:horizon2) farr2[i, j] <- forec$mean[i]
        parr2[j, ] <- accuracy(forec, x[nmh23p1:nmh3, j])
    }
    x2 <- x1
    for (j in 1:numseries) {
        for (i in 1:horizon2) {
            x2[nmh23 + i, j] = farr2[i, j]
        }
    }
    for (j in 1:numseries) {
        fit <- auto.arima(ts(x2[1:nmh3, j], freq = s), d = d, 
            D = D)
        forec <- forecast(fit, horizon3)
        for (i in 1:horizon3) rarr3[i, j] <- forec$residuals[i]
        forec$mean <- forec$mean + 0.5 * sd(rarr3[, j])
        forec$resid <- forec$resid - 0.5 * sd(rarr3[, j])
        for (i in 1:horizon3) farr3[i, j] <- forec$mean[i]
        parr3[j, ] <- accuracy(forec, x[nmh3p1:n, j])
    }
    x3 <- x2
    for (j in 1:numseries) {
        for (i in 1:horizon3) {
            x3[nmh3 + i, j] = farr3[i, j]
        }
    }
    mydum <- array(NA, dim = c(numseries, 3))
    myerror <- array(NA, dim = numseries)
    for (i in 1:numseries) {
        mydum[i, 1] = (mean(farr1[, i]) - mean(x[nmh123p1:nmh23, 
            i]))
        mydum[i, 2] = (mean(farr2[, i]) - mean(x[nmh23p1:nmh3, 
            i]))
        mydum[i, 3] = (mean(farr3[, i]) - mean(x[nmh3p1:n, i]))
        myerror[i] = mydum[i, 3]
    }
    for (i in 1:numseries) {
        if (myerror[i] > quantile(myerror, 0.95)) {
            print(paste("sigma = ", iteratie, "   horizon = ", 
                iteratie2, "   i = ", i, sep = ""))
            print(paste("mean Ft =", mean(farr3[, i]), sep = ""))
            print(paste("mean Yt =", mean(x[nmh3p1:n, i]), sep = ""))
            print(paste("mean et = ", mean(farr3[, i]) - mean(x[nmh3p1:n, 
                i]), sep = ""))
            print(paste("%error = ", abs((mean(farr3[, i]) - 
                mean(x[nmh3p1:n, i]))/mean(x[nmh3p1:n, i])), 
                sep = ""))
        }
    }
    i = 1
    
    plot(x[nmh123:n, i], type = "l", notch = T, ylab = "Order quantity", 
        xlab = "Period")
    lines(x1[nmh123:n, i], lty = 2)
    lines(x2[nmh123:n, i], lty = 3)
    lines(x3[nmh123:n, i], lty = 4)
    
    dum <- x[nmh123:n, ] - x3[nmh123:n, ]
    dum <- x[1:n, ] - x3[1:n, ]
    dum1 <- array(NA, dim = c(horizon1, numseries))
    dum1 <- dum[nmh123p1:nmh23, ]
    dum2 <- array(NA, dim = c(horizon2, numseries))
    dum2 <- dum[nmh23p1:nmh3, ]
    dum3 <- array(NA, dim = c(horizon3, numseries))
    dum3 <- dum[nmh3p1:n, ]
    
    boxplot(cbind(t(dum1), t(dum2), t(dum3)), notch = T, ylab = "Demand error", 
        xlab = "Period")
    abline(v = horizon1 + 0.5)
    abline(v = horizon1 + horizon2 + 0.5)
    abline(v = horizon1 + horizon2 + horizon3 + 0.5)
    
}