Replication: Weiß and Kim (2025)

This page reproduces Table 8 of Weiß and Kim (2025). The table reports simulated rejection rates at level \(\alpha = 0.05\) for DGP 4, an SQMA(1,1) model with \((\beta_1, \beta_2, \beta_3) = (0.8, 0.8, 0.8)\), for the classical spatial ordinal pattern (SOP) classification and for the three refined ones.

Libraries

using StatsOrdinalPatterns, Random, Statistics, Printf, DataFrames, Distributions

Statistics of the paper and their types in the package

The four column groups of Table 8 are the refinement argument of test_sop:

Table 8 group refinement
Ordinary types OrdinaryType()
Rotation types RotationType()
Direction types DirectionType()
Diagonal types DiagonalType()

The columns within a group are its chart_choice:

Table 8 column chart_choice Rejects \(H_0\) when Label below
\(\tilde{\tau}\) TauTilde() \(\lvert\tilde{\tau}\rvert > c\) tau_tilde
\(\widehat{H}\) Shannon(base=exp(1)) rescaled \(\widehat{H} > c\) H
\(\widehat{H}_{ex}\) ShannonExtropy(base=exp(1)) rescaled \(\widehat{H}_{ex} > c\) H_ex
\(\widehat{\Delta}\) DistanceToWhiteNoise() \(\widehat{\Delta} > c\) Delta

\(\tilde{\tau}\) appears in the ordinary group only; the refined classifications have asymptotic theory for the entropy type charts. The rescaling of the entropy statistics and the generalized \(\chi^2\) critical value are what test_sop applies internally, and crit_val_sop returns the critical value for a given scheme, level and field size.

Setup

DGP 4 is the SQMA11 type of the package: dgp_params are the \(\beta_j\), eps_params are the exponents \((a, b, c)\) of the three error terms, written “\(1^a\,2^b\,3^c\)” in the paper, and M_rows, N_cols are the dimensions of the field. m and n in the table are the dimensions of the SOP matrix, so the field is \((m+1) \times (n+1)\).

Parameters of Table 8 and one field of DGP 4 
β = (0.8, 0.8, 0.8)                                # dependence parameters of Table 8
level = 0.05                                       # nominal significance level
sizes = ((10, 10), (15, 15), (20, 20), (40, 25))   # SOP matrix dimensions (m, n)
reps = 10_000                                      # Monte Carlo replications

models = [
    ("DGP4-1", "1² 2² 3²", (2, 2, 2)),
    ("DGP4-2", "1² 2¹ 3²", (2, 1, 2)),
    ("DGP4-3", "1¹ 2¹ 3²", (1, 1, 2)),
    ("DGP4-4", "1² 2¹ 3¹", (2, 1, 1)),
]

entropies = [("H", Shannon(base=exp(1))), ("H_ex", ShannonExtropy(base=exp(1))), ("Delta", DistanceToWhiteNoise())]

schemes = [
    ("ord", OrdinaryType(), vcat([("tau_tilde", TauTilde())], entropies)),
    ("rot", RotationType(), entropies),
    ("dir", DirectionType(), entropies),
    ("diag", DiagonalType(), entropies),
]

stats = ["$(sn)_$(cn)" for (sn, _, cs) in schemes for (cn, _) in cs]   # 13 columns

# One field of the DGP. `fill_mat_dgp_sop!` is the routine the ARL functions of the
# package use internally; it writes into the caller's buffers and returns a view of the
# field. `mat` receives the field, `mat_ma` the errors, `mat_ao` the additive outliers,
# which DGP 4 does not have.
sqma11(eps, m, n) = SQMA11(β, eps, m + 1, n + 1, Normal(0, 1), nothing)

buffers(dgp) = ntuple(_ -> zeros(dgp.M_rows + 1, dgp.N_cols + 1), 3)

draw_field!(dgp, mat, mat_ao, mat_ma) =
    StatsOrdinalPatterns.fill_mat_dgp_sop!(dgp, dgp.dist, dgp.dist_ao, mat, mat_ao, mat_ma)

A single test

One field of DGP 4 and one SOP test 
Random.seed!(1)
dgp = sqma11((2, 2, 2), 10, 10)
Y = collect(draw_field!(dgp, buffers(dgp)...))   # copy out of the buffer
test_sop(Y, 1, 1; chart_choice=ShannonExtropy(base=exp(1)), refinement=DiagonalType(), alpha=level)
SOPTestResult
  Chart:            ShannonExtropy(base = 2.718281828459045)
  Statistic:        0.0169
  ─────────────────────────────
  Asymptotic test
    Critical value: 0.015
    p-value:        0.0309
    Reject H₀:      true

Running the power study

crit_val_sop returns the critical value alone, which depends on the field size, the level and the classification scheme, not on the data, so it is computed once per row of the table. test_sop additionally computes the \(p\) value, which for the entropy charts means evaluating a generalized \(\chi^2\) distribution numerically and which a power study does not use. Each replication then costs one stat_sop pass per classification scheme, and the charts of that scheme are read off the resulting \(\hat{\mathbf{p}}\) with chart_stat_sop.

Rejection rates of the 13 tests for one field size 
# two sided for τ̃, upper tail on the rescaled statistic for the entropy charts
rejects(::TauTilde, p̂, crit) = abs(chart_stat_sop(p̂, TauTilde())) > crit
rejects(cc, p̂, crit) =
    StatsOrdinalPatterns.rescale_sop(chart_stat_sop(p̂, cc), length(p̂), cc) > crit

function power(m, n, eps, reps; level=0.05)
    dgp = sqma11(eps, m, n)
    M, N = dgp.M_rows, dgp.N_cols
    crits = [crit_val_sop(M, N, 1, 1; chart_choice=cc, refinement=ref, alpha=level)
             for (_, ref, cs) in schemes for (_, cc) in cs]

    tasks = map(Iterators.partition(1:reps, cld(reps, Threads.nthreads()))) do idx
        Threads.@spawn begin
            bufs = buffers(dgp)                     # one set of buffers per task
            hits = zeros(Int, length(crits))
            for _ in idx
                Y = draw_field!(dgp, bufs...)
                k = 0
                for (_, ref, cs) in schemes
                    p̂ = stat_sop(Y, 1, 1; chart_choice=Shannon(base=exp(1)), refinement=ref)[2]
                    for (_, cc) in cs
                        k += 1
                        hits[k] += rejects(cc, p̂, crits[k])
                    end
                end
            end
            hits
        end
    end
    return sum(fetch.(tasks)) ./ reps
end

The shortcut has to give the same decision as test_sop in all 13 columns, not merely a similar one:

Check against test_sop 
let dgp = sqma11((1, 1, 2), 20, 20), M = 21, N = 21
    Y = collect(draw_field!(dgp, buffers(dgp)...))
    ok = true
    for (_, ref, cs) in schemes
        p̂ = stat_sop(Y, 1, 1; chart_choice=Shannon(base=exp(1)), refinement=ref)[2]
        for (_, cc) in cs
            r = test_sop(Y, 1, 1; chart_choice=cc, refinement=ref, alpha=level)
            c = crit_val_sop(M, N, 1, 1; chart_choice=cc, refinement=ref, alpha=level)
            ok &= c ≈ r.asymp_crit && rejects(cc, p̂, c) == r.asymp_reject
        end
    end
    (shortcut_matches_test_sop=ok,)
end
(shortcut_matches_test_sop = true,)
Note

With reps = 10_000 the page runs in a few seconds and the standard error of a rejection rate is at most \(\sqrt{0.25 / 10^4} = 0.005\).

Table 8: DGP 4, \((\beta_1, \beta_2, \beta_3) = (0.8, 0.8, 0.8)\)

Run the study for every model and field size 
res = Dict(
    (name, m, n) => power(m, n, eps, reps; level=level)
    for (name, _, eps) in models for (m, n) in sizes
)

Ordinary types

16×8 DataFrame
Row model label m n tau_tilde H H_ex Delta
String String Int64 Int64 Float64 Float64 Float64 Float64
1 DGP4-1 1² 2² 3² 10 10 0.25 0.196 0.185 0.194
2 DGP4-1 1² 2² 3² 15 15 0.469 0.424 0.41 0.409
3 DGP4-1 1² 2² 3² 20 20 0.766 0.687 0.672 0.675
4 DGP4-1 1² 2² 3² 40 25 0.991 0.98 0.978 0.978
5 DGP4-2 1² 2¹ 3² 10 10 0.591 0.491 0.469 0.485
6 DGP4-2 1² 2¹ 3² 15 15 0.926 0.889 0.876 0.876
7 DGP4-2 1² 2¹ 3² 20 20 0.997 0.994 0.993 0.993
8 DGP4-2 1² 2¹ 3² 40 25 1.0 1.0 1.0 1.0
9 DGP4-3 1¹ 2¹ 3² 10 10 0.382 0.312 0.301 0.31
10 DGP4-3 1¹ 2¹ 3² 15 15 0.702 0.666 0.651 0.65
11 DGP4-3 1¹ 2¹ 3² 20 20 0.943 0.913 0.907 0.908
12 DGP4-3 1¹ 2¹ 3² 40 25 1.0 1.0 1.0 1.0
13 DGP4-4 1² 2¹ 3¹ 10 10 0.215 0.17 0.162 0.167
14 DGP4-4 1² 2¹ 3¹ 15 15 0.419 0.382 0.37 0.368
15 DGP4-4 1² 2¹ 3¹ 20 20 0.694 0.62 0.607 0.609
16 DGP4-4 1² 2¹ 3¹ 40 25 0.974 0.957 0.955 0.955

Rotation types

16×7 DataFrame
Row model label m n H H_ex Delta
String String Int64 Int64 Float64 Float64 Float64
1 DGP4-1 1² 2² 3² 10 10 0.192 0.173 0.172
2 DGP4-1 1² 2² 3² 15 15 0.379 0.356 0.358
3 DGP4-1 1² 2² 3² 20 20 0.633 0.614 0.616
4 DGP4-1 1² 2² 3² 40 25 0.967 0.965 0.965
5 DGP4-2 1² 2¹ 3² 10 10 0.453 0.414 0.418
6 DGP4-2 1² 2¹ 3² 15 15 0.845 0.823 0.826
7 DGP4-2 1² 2¹ 3² 20 20 0.987 0.984 0.985
8 DGP4-2 1² 2¹ 3² 40 25 1.0 1.0 1.0
9 DGP4-3 1¹ 2¹ 3² 10 10 0.293 0.275 0.276
10 DGP4-3 1¹ 2¹ 3² 15 15 0.608 0.587 0.589
11 DGP4-3 1¹ 2¹ 3² 20 20 0.877 0.868 0.869
12 DGP4-3 1¹ 2¹ 3² 40 25 1.0 1.0 1.0
13 DGP4-4 1² 2¹ 3¹ 10 10 0.174 0.156 0.155
14 DGP4-4 1² 2¹ 3¹ 15 15 0.35 0.332 0.335
15 DGP4-4 1² 2¹ 3¹ 20 20 0.567 0.552 0.554
16 DGP4-4 1² 2¹ 3¹ 40 25 0.938 0.934 0.935

Direction types

16×7 DataFrame
Row model label m n H H_ex Delta
String String Int64 Int64 Float64 Float64 Float64
1 DGP4-1 1² 2² 3² 10 10 0.124 0.129 0.124
2 DGP4-1 1² 2² 3² 15 15 0.274 0.271 0.269
3 DGP4-1 1² 2² 3² 20 20 0.504 0.5 0.498
4 DGP4-1 1² 2² 3² 40 25 0.934 0.931 0.931
5 DGP4-2 1² 2¹ 3² 10 10 0.854 0.868 0.862
6 DGP4-2 1² 2¹ 3² 15 15 0.998 0.998 0.998
7 DGP4-2 1² 2¹ 3² 20 20 1.0 1.0 1.0
8 DGP4-2 1² 2¹ 3² 40 25 1.0 1.0 1.0
9 DGP4-3 1¹ 2¹ 3² 10 10 0.204 0.208 0.201
10 DGP4-3 1¹ 2¹ 3² 15 15 0.481 0.474 0.471
11 DGP4-3 1¹ 2¹ 3² 20 20 0.796 0.788 0.788
12 DGP4-3 1¹ 2¹ 3² 40 25 0.998 0.998 0.998
13 DGP4-4 1² 2¹ 3¹ 10 10 0.124 0.13 0.124
14 DGP4-4 1² 2¹ 3¹ 15 15 0.269 0.269 0.265
15 DGP4-4 1² 2¹ 3¹ 20 20 0.481 0.478 0.477
16 DGP4-4 1² 2¹ 3¹ 40 25 0.912 0.91 0.91

Diagonal types

16×7 DataFrame
Row model label m n H H_ex Delta
String String Int64 Int64 Float64 Float64 Float64
1 DGP4-1 1² 2² 3² 10 10 0.265 0.269 0.261
2 DGP4-1 1² 2² 3² 15 15 0.543 0.551 0.548
3 DGP4-1 1² 2² 3² 20 20 0.819 0.825 0.823
4 DGP4-1 1² 2² 3² 40 25 0.998 0.998 0.998
5 DGP4-2 1² 2¹ 3² 10 10 0.462 0.444 0.438
6 DGP4-2 1² 2¹ 3² 15 15 0.853 0.843 0.845
7 DGP4-2 1² 2¹ 3² 20 20 0.989 0.988 0.988
8 DGP4-2 1² 2¹ 3² 40 25 1.0 1.0 1.0
9 DGP4-3 1¹ 2¹ 3² 10 10 0.334 0.33 0.325
10 DGP4-3 1¹ 2¹ 3² 15 15 0.676 0.674 0.675
11 DGP4-3 1¹ 2¹ 3² 20 20 0.92 0.921 0.921
12 DGP4-3 1¹ 2¹ 3² 40 25 1.0 1.0 1.0
13 DGP4-4 1² 2¹ 3¹ 10 10 0.228 0.23 0.223
14 DGP4-4 1² 2¹ 3¹ 15 15 0.468 0.476 0.474
15 DGP4-4 1² 2¹ 3¹ 20 20 0.737 0.745 0.744
16 DGP4-4 1² 2¹ 3¹ 40 25 0.99 0.991 0.99

The paper prints the largest power value of each row in bold:

16×5 DataFrame
Row model m n max_stat rate
String Int64 Int64 String Float64
1 DGP4-1 10 10 diag_H_ex 0.269
2 DGP4-1 15 15 diag_H_ex 0.551
3 DGP4-1 20 20 diag_H_ex 0.825
4 DGP4-1 40 25 diag_H, diag_H_ex, diag_Delta 0.998
5 DGP4-2 10 10 dir_H_ex 0.868
6 DGP4-2 15 15 dir_H, dir_H_ex, dir_Delta 0.998
7 DGP4-2 20 20 dir_H, dir_H_ex, dir_Delta 1.0
8 DGP4-2 40 25 13 columns tied 1.0
9 DGP4-3 10 10 ord_tau_tilde 0.382
10 DGP4-3 15 15 ord_tau_tilde 0.702
11 DGP4-3 20 20 ord_tau_tilde 0.943
12 DGP4-3 40 25 10 columns tied 1.0
13 DGP4-4 10 10 diag_H_ex 0.23
14 DGP4-4 15 15 diag_H_ex 0.476
15 DGP4-4 20 20 diag_H_ex 0.745
16 DGP4-4 40 25 diag_H_ex 0.991

Comparison with the published rates

entries compared:           208
largest absolute deviation: 0.012
mean absolute deviation:    0.003
Monte Carlo standard error: at most 0.005

The five largest deviations:

5×7 DataFrame
Row model m n statistic paper replicated dev
String Int64 Int64 String Float64 Float64 Float64
1 DGP4-4 20 20 dir_H 0.469 0.481 0.012
2 DGP4-3 10 10 ord_H_ex 0.29 0.301 0.011
3 DGP4-4 15 15 dir_H 0.258 0.269 0.011
4 DGP4-2 15 15 ord_tau_tilde 0.916 0.926 0.01
5 DGP4-3 10 10 ord_Delta 0.301 0.31 0.01

References

Weiß, Christian H, and Hee-Young Kim. 2025. “Non-Parametric Entropy Tests for Spatial Dependence.” Computational Statistics, 1–38.