Replication: Weiß and Testik (2023)

This page reproduces Table 1 and Table 2 of Weiß and Testik (2023) with StatsOrdinalPatterns.jl. Both tables are Monte Carlo results, so replicating them checks the ordinal-pattern coding, the EWMA recursion, the alarm rules and the root finder that calibrates the charts, all at once.

Everything below is computed from scratch; the published numbers are hard-coded only so that they can be printed next to the replicated ones.

Libraries

using StatsOrdinalPatterns, Random, Statistics, Printf, DataFrames

Statistics of the paper and their types in the package

The paper monitors six statistics of the EWMA-smoothed pattern distribution \(\hat{\mathbf{p}}_t\) (Equation (5) of Weiß and Testik (2023)): the three omnibus measures of Equation (3) and the three directed measures of Bandt (2019) given in Equation (4). Each of them is selected in the package through the chart_choice keyword:

Paper Statistic chart_choice Alarm rule
Eq. (3) \(\widehat{H}^{(d)}\), permutation entropy Shannon(base=exp(1)) \(\widehat{H} < cl\)
Eq. (3) \(\widehat{H}_{\text{ex}}^{(d)}\), permutation extropy ShannonExtropy(base=exp(1)) \(\widehat{H}_{\text{ex}} < cl\)
Eq. (3) \(\widehat{\Delta}^{(d)}\), distance to white noise DistanceToWhiteNoise() \(\widehat{\Delta} > cl\)
Eq. (4) \(\hat{\beta}^{(d)}\), up-down balance UpDownBalance() \(\lvert\hat{\beta}\rvert > cl\)
Eq. (4) \(\hat{\tau}^{(d)}\), persistence Persistence() \(\lvert\hat{\tau}\rvert > cl\)
Eq. (4) \(\hat{\delta}^{(d)}\), up-down scaling UpDownScaling() \(\lvert\hat{\delta}\rvert > cl\)

Setup

The in-control process is i.i.d. standard normal. Since all six statistics depend on the data only through ranks, the resulting designs apply to any continuously distributed i.i.d. process — the distribution-free property noted below Table 1 of the paper.

ic = ContinuousDGPIC(Normal(0, 1))     # in-control process
L0 = 370.0                             # target in-control ARL
m, d = 3, 1                            # pattern length and delay
lams = (0.25, 0.10, 0.05)              # smoothing parameters of Table 1
reps = 10_000                          # replications per ARL estimate

charts = [
    ("H", Shannon(base=exp(1))),
    ("Hex", ShannonExtropy(base=exp(1))),
    ("Δ", DistanceToWhiteNoise()),
    ("β", UpDownBalance()),
    ("τ", Persistence()),
    ("δ", UpDownScaling()),
]
Note

With reps = 10_000 the whole page runs in well under a minute on a multicore machine, and the standard error of an ARL estimate is roughly 1 % of the ARL. That is enough to confirm every entry of both tables, but the ARLs move by about a percent from render to render. Raise reps for tighter agreement.

Table 1: in-control designs

cl_op searches for the critical limit at which the in-control ARL equals L0: it first brackets the root with a coarse Monte Carlo estimate and then refines the bracket with the ITP algorithm. Only a rough starting value is needed — the ones below are deliberately round numbers that differ from the published limits.

# (λ, chart) => (starting value, step size of the bracket search)
start = Dict(
    (0.25, "H") => (1.05, 0.01), (0.25, "Hex") => (0.70, 0.01), (0.25, "Δ") => (0.30, 0.01),
    (0.25, "β") => (0.60, 0.01), (0.25, "τ") => (0.40, 0.01), (0.25, "δ") => (0.70, 0.02),
    (0.10, "H") => (1.50, 0.01), (0.10, "Hex") => (0.85, 0.005), (0.10, "Δ") => (0.10, 0.005),
    (0.10, "β") => (0.35, 0.005), (0.10, "τ") => (0.25, 0.005), (0.10, "δ") => (0.45, 0.01),
    (0.05, "H") => (1.70, 0.02), (0.05, "Hex") => (0.89, 0.005), (0.05, "Δ") => (0.045, 0.002),
    (0.05, "β") => (0.22, 0.005), (0.05, "τ") => (0.16, 0.005), (0.05, "δ") => (0.30, 0.01),
)

cls = Dict(
    (lam, name) => cl_op(
        ic, lam, L0, start[(lam, name)][1];
        chart_choice=chart, m=m, d=d,
        reps_bracket=reps, reps_final=reps,
        bracket_step=start[(lam, name)][2], seed=42
    )
    for lam in lams for (name, chart) in charts
)

The limits are then verified by an independent run of arl_op_ic.

arls = Dict(
    (lam, name) => arl_op_ic(ic, lam, cls[(lam, name)], reps; chart_choice=chart, m=m, d=d)
    for lam in lams for (name, chart) in charts
)

Replicated table

The layout follows Table 1: for every smoothing parameter, the critical limit of each chart and the in-control ARL it achieves.

3×13 DataFrame
Row λ l_H ARL_H l_Hex ARL_Hex l_Δ ARL_Δ l_β ARL_β l_τ ARL_τ l_δ ARL_δ
Float64 Float64 Float64 Float64 Float64 Float64 Float64 Float64 Float64 Float64 Float64 Float64 Float64
1 0.25 1.0128 376.5 0.66182 372.7 0.33442 372.9 0.64435 377.0 0.42556 370.7 0.76511 362.5
2 0.1 1.4605 368.3 0.84063 367.8 0.11138 364.3 0.36399 371.1 0.2532 374.5 0.48763 376.2
3 0.05 1.6358 365.4 0.88015 365.0 0.051294 372.4 0.23319 373.8 0.16772 363.3 0.32447 362.5

Comparison with the published limits

Published limits and comparison table 
published_cl = Dict(
    (0.25, "H") => 1.014, (0.25, "Hex") => 0.6621, (0.25, "Δ") => 0.3338,
    (0.25, "β") => 0.6437, (0.25, "τ") => 0.4253, (0.25, "δ") => 0.7656,
    (0.10, "H") => 1.4601, (0.10, "Hex") => 0.8405, (0.10, "Δ") => 0.1115,
    (0.10, "β") => 0.3638, (0.10, "τ") => 0.2529, (0.10, "δ") => 0.4876,
    (0.05, "H") => 1.6356, (0.05, "Hex") => 0.88017, (0.05, "Δ") => 0.05125,
    (0.05, "β") => 0.2330, (0.05, "τ") => 0.16775, (0.05, "δ") => 0.3246,
)

cmp1 = DataFrame(
    λ=[lam for lam in lams for _ in charts],
    chart=[name for _ in lams for (name, _) in charts],
    cl_paper=[published_cl[(lam, name)] for lam in lams for (name, _) in charts],
    cl_replicated=[round(cls[(lam, name)], sigdigits=5) for lam in lams for (name, _) in charts],
    ARL=[round(arls[(lam, name)][1], digits=1) for lam in lams for (name, _) in charts],
    se=[round(arls[(lam, name)][2], digits=1) for lam in lams for (name, _) in charts],
)
cmp1.dev_pct = round.(100 .* (cmp1.cl_replicated .- cmp1.cl_paper) ./ cmp1.cl_paper, digits=2)
cmp1
18×7 DataFrame
Row λ chart cl_paper cl_replicated ARL se dev_pct
Float64 String Float64 Float64 Float64 Float64 Float64
1 0.25 H 1.014 1.0128 376.5 3.7 -0.12
2 0.25 Hex 0.6621 0.66182 372.7 3.6 -0.04
3 0.25 Δ 0.3338 0.33442 372.9 3.7 0.19
4 0.25 β 0.6437 0.64435 377.0 3.7 0.1
5 0.25 τ 0.4253 0.42556 370.7 3.7 0.06
6 0.25 δ 0.7656 0.76511 362.5 3.5 -0.06
7 0.1 H 1.4601 1.4605 368.3 3.5 0.03
8 0.1 Hex 0.8405 0.84063 367.8 3.6 0.02
9 0.1 Δ 0.1115 0.11138 364.3 3.5 -0.11
10 0.1 β 0.3638 0.36399 371.1 3.6 0.05
11 0.1 τ 0.2529 0.2532 374.5 3.7 0.12
12 0.1 δ 0.4876 0.48763 376.2 3.7 0.01
13 0.05 H 1.6356 1.6358 365.4 3.5 0.01
14 0.05 Hex 0.88017 0.88015 365.0 3.5 -0.0
15 0.05 Δ 0.05125 0.051294 372.4 3.6 0.09
16 0.05 β 0.233 0.23319 373.8 3.6 0.08
17 0.05 τ 0.16775 0.16772 363.3 3.5 -0.02
18 0.05 δ 0.3246 0.32447 362.5 3.5 -0.04 
largest relative deviation of a critical limit: 0.19 %
in-control ARLs: min 362.5, max 377.0 (target 370)

All eighteen limits agree with the published ones to within a few tenths of a percent, and every chart holds the target ARL of 370 within Monte Carlo error.

Table 2: out-of-control performance under DGP 1

DGP 1 of the paper is the Gaussian AR(1) process \(X_t = \alpha X_{t-1} + \epsilon_t\) with \(\epsilon_t \sim \text{i.i.d. } N(0,1)\) and \(\alpha \in \{\pm 0.2, \pm 0.4, \pm 0.6, \pm 0.8\}\). It is constructed by AR1(α, Normal(0, 1)) and started from its stationary distribution.

Table 2 reports the four charts that react to this alternative — \(H\), \(H_{\text{ex}}\), \(\Delta\) and \(\tau\). The paper omits the \(\beta\)- and \(\delta\)-charts here because both measure asymmetries that a Gaussian AR(1) process does not create.

The charts keep the designs computed above, so this section also propagates any error in the critical limits.

charts2 = [("H", Shannon(base=exp(1))), ("Hex", ShannonExtropy(base=exp(1))),
    ("Δ", DistanceToWhiteNoise()), ("τ", Persistence())]
alphas = (0.2, 0.4, 0.6, 0.8, -0.2, -0.4, -0.6, -0.8)

oc = Dict(
    (lam, α, name) => arl_op_oc(
        AR1(α, Normal(0, 1)), lam, cls[(lam, name)], reps;
        chart_choice=chart, m=m, d=d
    )[1]
    for lam in lams for α in alphas for (name, chart) in charts2
)
tab2 (generic function with 1 method)

\(\lambda = 0.25\)

8×5 DataFrame
Row α H Hex Δ τ
Float64 Float64 Float64 Float64 Float64
1 0.2 237.1 217.2 220.6 189.1
2 0.4 144.9 127.0 129.4 103.6
3 0.6 87.4 76.2 76.0 60.2
4 0.8 50.5 44.0 44.9 39.6
5 -0.2 595.4 649.9 648.1 800.7
6 -0.4 942.4 1177.7 1163.5 1835.0
7 -0.6 1289.3 2338.2 2375.7 4656.2
8 -0.8 512.9 5629.8 5265.0 13596.5

\(\lambda = 0.10\)

8×5 DataFrame
Row α H Hex Δ τ
Float64 Float64 Float64 Float64 Float64
1 0.2 300.5 231.6 240.0 192.9
2 0.4 212.4 144.1 153.7 100.2
3 0.6 134.7 88.4 94.9 59.0
4 0.8 81.3 54.7 57.4 40.2
5 -0.2 410.1 553.0 517.9 412.5
6 -0.4 388.4 766.8 669.9 202.4
7 -0.6 240.5 675.8 563.3 83.6
8 -0.8 84.7 229.0 189.2 37.5

\(\lambda = 0.05\)

8×5 DataFrame
Row α H Hex Δ τ
Float64 Float64 Float64 Float64 Float64
1 0.2 298.4 260.0 270.7 202.8
2 0.4 206.0 165.8 172.8 105.3
3 0.6 133.1 103.0 107.9 63.6
4 0.8 84.2 67.8 70.1 43.5
5 -0.2 348.1 420.7 407.6 273.1
6 -0.4 236.1 344.0 327.3 118.6
7 -0.6 118.9 183.8 166.9 57.1
8 -0.8 57.0 75.5 71.1 30.8

Comparison with the published ARLs

Published ARLs and comparison table 
# published values, ordered (H, Hex, Δ, τ)
published_oc = Dict(
    (0.25, 0.2) => (231.9, 217.8, 217.6, 185.4), (0.25, 0.4) => (143.7, 128.0, 128.5, 101.8),
    (0.25, 0.6) => (85.9, 75.2, 75.3, 60.4), (0.25, 0.8) => (50.6, 44.7, 44.8, 39.9),
    (0.25, -0.2) => (588.0, 647.8, 640.2, 798.2), (0.25, -0.4) => (942.4, 1174.4, 1162.0, 1840.3),
    (0.25, -0.6) => (1255.2, 2324.6, 2266.1, 4634.4), (0.25, -0.8) => (491.8, 5508.4, 5242.5, 13649.7),
    (0.10, 0.2) => (299.8, 233.7, 242.2, 191.3), (0.10, 0.4) => (212.4, 145.1, 152.9, 99.9),
    (0.10, 0.6) => (135.3, 89.3, 94.2, 59.7), (0.10, 0.8) => (81.7, 54.7, 57.6, 40.5),
    (0.10, -0.2) => (411.1, 570.5, 531.3, 411.2), (0.10, -0.4) => (392.8, 786.8, 681.8, 203.7),
    (0.10, -0.6) => (242.3, 692.3, 565.1, 83.8), (0.10, -0.8) => (85.7, 231.6, 190.8, 37.4),
    (0.05, 0.2) => (303.3, 261.2, 269.5, 206.2), (0.05, 0.4) => (208.3, 166.2, 172.9, 105.7),
    (0.05, 0.6) => (132.4, 104.0, 108.7, 63.8), (0.05, 0.8) => (84.5, 67.2, 69.4, 43.9),
    (0.05, -0.2) => (349.0, 424.2, 411.2, 271.4), (0.05, -0.4) => (237.8, 347.7, 326.2, 118.0),
    (0.05, -0.6) => (120.6, 180.6, 168.2, 56.2), (0.05, -0.8) => (57.0, 75.1, 71.3, 30.7),
)

dev = DataFrame(λ=Float64[], α=Float64[], chart=String[],
    ARL_paper=Float64[], ARL_replicated=Float64[], dev_pct=Float64[])
for lam in lams, α in alphas
    p = published_oc[(lam, α)]
    for (i, (name, _)) in enumerate(charts2)
        r = oc[(lam, α, name)]
        push!(dev, (lam, α, name, p[i], round(r, digits=1), round(100 * (r - p[i]) / p[i], digits=2)))
    end
end

DataFrame(
    λ=collect(lams),
    max_abs_dev_pct=[round(maximum(abs, dev[dev.λ.==lam, :dev_pct]), digits=2) for lam in lams],
    mean_abs_dev_pct=[round(mean(abs, dev[dev.λ.==lam, :dev_pct]), digits=2) for lam in lams],
)
3×3 DataFrame
Row λ max_abs_dev_pct mean_abs_dev_pct
Float64 Float64 Float64
1 0.25 4.84 1.15
2 0.1 3.06 0.9
3 0.05 1.79 0.74

The five largest deviations across all 96 entries:

5×6 DataFrame
Row λ α chart ARL_paper ARL_replicated dev_pct
Float64 Float64 String Float64 Float64 Float64
1 0.25 -0.6 Δ 2266.1 2375.7 4.84
2 0.25 -0.8 H 491.8 512.9 4.3
3 0.1 -0.2 Hex 570.5 553.0 -3.06
4 0.25 -0.6 H 1255.2 1289.3 2.72
5 0.1 -0.4 Hex 786.8 766.8 -2.54

Every cell is reproduced to within a few percent, which is the resolution that 10 000 replications provide. The largest deviations occur exactly where they should: in the cells with ARLs of several thousand, where the run-length distribution is most dispersed.

References

Bandt, Christoph. 2019. “Small Order Patterns in Big Time Series: A Practical Guide.” Entropy 21 (6): 613. https://doi.org/10.3390/e21060613.
Weiß, Christian H, and Murat Caner Testik. 2023. “Nonparametric Control Charts for Monitoring Serial Dependence Based on Ordinal Patterns.” Technometrics 65 (3): 340–50.