Replication: Weiß and Kim (2024)

This page reproduces Table B.8 of Weiß and Kim (2024). The table reports simulated rejection rates of the four asymptotic spatial ordinal pattern (SOP) tests at level \(\alpha = 0.05\), together with the spatial autocorrelation at lag \(\mathbf{1}\) as a benchmark, for DGP 8 and growing field size.

Libraries

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

Statistics of the paper and their types in the package

Each column of Table B.8 corresponds to one chart_choice of test_sop, which returns the statistic, the asymptotic critical value, the \(p\) value and the reject decision. The benchmark column is test_sacf.

Table B.8 column	Package call	Rejects \(H_0\) when	Label below
\(\hat{\tau}\)	test_sop(data, 1, 1; chart_choice=TauHat())	\(\lvert\hat{\tau}\rvert > c\)	tau_hat
\(\hat{\kappa}\)	test_sop(data, 1, 1; chart_choice=KappaHat())	\(\lvert\hat{\kappa}\rvert > c\)	kappa_hat
\(\tilde{\tau}\)	test_sop(data, 1, 1; chart_choice=TauTilde())	\(\lvert\tilde{\tau}\rvert > c\)	tau_tilde
\(\tilde{\kappa}\)	test_sop(data, 1, 1; chart_choice=KappaTilde())	\(\lvert\tilde{\kappa}\rvert > c\)	kappa_tilde
\(\hat{\rho}(\mathbf{1})\)	test_sacf(data, 1, 1)	\(\lvert\hat{\rho}\rvert > c\)	rho_1

All four SOP statistics are linear in the same vector of type frequencies \(\hat{\mathbf{p}} = (\hat p_1, \hat p_2, \hat p_3)\) of the classical SOP classification, which stat_sop returns as its second output.

The header of the table asks for critical values computed with the exact variance. crit_val_sop uses exactly that: its variances carry the finite sample factor \(1 - 1/(2m) - 1/(2n)\), so no asymptotic simplification enters the critical values below.

Setup

DGP 8 is the spatial quadratic moving average process SQMA(1,1) with \((\beta_1, \beta_2, \beta_3) = (0.8, 0.8, 0.8)\) and i.i.d. standard normal errors. The SQMA11 type of the package holds it: 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 simulated field.

m and n in the table are the dimensions of the SOP matrix, so the field is \((m+1) \times (n+1)\).

β = (0.8, 0.8, 0.8)                                # dependence parameters of Table B.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 = [
    ("DGP8-1", "1² 2² 3²", (2, 2, 2)),
    ("DGP8-2", "1² 2¹ 3²", (2, 1, 2)),
    ("DGP8-3", "1¹ 2¹ 3²", (1, 1, 2)),
    ("DGP8-4", "1² 2¹ 3¹", (2, 1, 1)),
]

charts = [
    ("tau_hat", TauHat()),
    ("kappa_hat", KappaHat()),
    ("tau_tilde", TauTilde()),
    ("kappa_tilde", KappaTilde()),
]

# 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 8 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

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=TauTilde(), alpha=level)

SOPTestResult
  Chart:            TauTilde()
  Statistic:        -0.0633
  ─────────────────────────────
  Asymptotic test
    Critical value: 0.1004
    p-value:        0.2162
    Reject H₀:      false

test_sacf(Y, 1, 1; alpha=level)

SACFTestResult
  Statistic:        -0.1425
  ─────────────────────────────
  Asymptotic test
    Critical value: 0.1782
    p-value:        0.1171
    Reject H₀:      false

Running the power study

crit_val_sop and crit_val_sacf return the critical values alone. They depend on the field size and the level, not on the data, so they are computed once per row of the table. test_sop additionally computes the \(p\) value, which a power study does not use. Each replication then costs one pass over the field: stat_sop counts the SOP frequencies once, and the four statistics are read off the resulting \(\hat{\mathbf{p}}\) with chart_stat_sop.

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, alpha=level) for (_, cc) in charts]
    crit_rho = crit_val_sacf(M, N; alpha=level)

    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(charts) + 1)
            for _ in idx
                Y = draw_field!(dgp, bufs...)
                p̂ = stat_sop(Y, 1, 1; chart_choice=Shannon(base=exp(1)))[2]
                for (j, (_, cc)) in enumerate(charts)
                    hits[j] += abs(chart_stat_sop(p̂, cc)) > crits[j]
                end
                hits[end] += abs(stat_sacf(Y, 1, 1)) > crit_rho
            end
            hits
        end
    end
    return sum(fetch.(tasks)) ./ reps
end

The shortcut has to give the same statistic and the same decision as test_sop and test_sacf, not merely a similar one:

let dgp = sqma11((1, 1, 2), 20, 20), M = 21, N = 21
    Y = collect(draw_field!(dgp, buffers(dgp)...))   # copy out of the buffer
    p̂ = stat_sop(Y, 1, 1; chart_choice=Shannon(base=exp(1)))[2]
    ok = all(
        begin
            r = test_sop(Y, 1, 1; chart_choice=cc, alpha=level)
            s = chart_stat_sop(p̂, cc)
            s ≈ r.stat && crit_val_sop(M, N, 1, 1; chart_choice=cc, alpha=level) ≈ r.asymp_crit &&
                (abs(s) > r.asymp_crit) == r.asymp_reject
        end for (_, cc) in charts
    )
    ok &= (abs(stat_sacf(Y, 1, 1)) > crit_val_sacf(M, N; alpha=level)) ==
          test_sacf(Y, 1, 1; alpha=level).asymp_reject
    (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 B.8: DGP 8, \((\beta_1, \beta_2, \beta_3) = (0.8, 0.8, 0.8)\)

res = Dict(
    (name, m, n) => power(m, n, eps, reps; level=level)
    for (name, _, eps) in models for (m, n) in sizes
)

16×9 DataFrame

Row	model	label	m	n	tau_hat	kappa_hat	tau_tilde	kappa_tilde	rho_1
	String	String	Int64	Int64	Float64	Float64	Float64	Float64	Float64
1	DGP8-1	1² 2² 3²	10	10	0.142	0.165	0.25	0.059	0.068
2	DGP8-1	1² 2² 3²	15	15	0.223	0.362	0.469	0.071	0.091
3	DGP8-1	1² 2² 3²	20	20	0.352	0.574	0.766	0.064	0.094
4	DGP8-1	1² 2² 3²	40	25	0.716	0.93	0.991	0.079	0.107
5	DGP8-2	1² 2¹ 3²	10	10	0.249	0.412	0.591	0.076	0.083
6	DGP8-2	1² 2¹ 3²	15	15	0.412	0.806	0.926	0.075	0.103
7	DGP8-2	1² 2¹ 3²	20	20	0.65	0.967	0.997	0.073	0.115
8	DGP8-2	1² 2¹ 3²	40	25	0.958	1.0	1.0	0.08	0.126
9	DGP8-3	1¹ 2¹ 3²	10	10	0.269	0.196	0.382	0.088	0.077
10	DGP8-3	1¹ 2¹ 3²	15	15	0.447	0.446	0.702	0.12	0.096
11	DGP8-3	1¹ 2¹ 3²	20	20	0.696	0.711	0.943	0.145	0.097
12	DGP8-3	1¹ 2¹ 3²	40	25	0.975	0.984	1.0	0.295	0.113
13	DGP8-4	1² 2¹ 3¹	10	10	0.161	0.126	0.215	0.064	0.468
14	DGP8-4	1² 2¹ 3¹	15	15	0.238	0.289	0.419	0.075	0.883
15	DGP8-4	1² 2¹ 3¹	20	20	0.399	0.453	0.694	0.09	0.991
16	DGP8-4	1² 2¹ 3¹	40	25	0.753	0.848	0.974	0.13	1.0

The paper prints the maximal rejection rate of each row in bold:

16×5 DataFrame

Row	model	m	n	max_stat	rate
	String	Int64	Int64	String	Float64
1	DGP8-1	10	10	tau_tilde	0.25
2	DGP8-1	15	15	tau_tilde	0.469
3	DGP8-1	20	20	tau_tilde	0.766
4	DGP8-1	40	25	tau_tilde	0.991
5	DGP8-2	10	10	tau_tilde	0.591
6	DGP8-2	15	15	tau_tilde	0.926
7	DGP8-2	20	20	tau_tilde	0.997
8	DGP8-2	40	25	kappa_hat, tau_tilde	1.0
9	DGP8-3	10	10	tau_tilde	0.382
10	DGP8-3	15	15	tau_tilde	0.702
11	DGP8-3	20	20	tau_tilde	0.943
12	DGP8-3	40	25	tau_tilde	1.0
13	DGP8-4	10	10	rho_1	0.468
14	DGP8-4	15	15	rho_1	0.883
15	DGP8-4	20	20	rho_1	0.991
16	DGP8-4	40	25	rho_1	1.0

Comparison with the published rates

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

The five largest deviations across all 80 entries:

5×7 DataFrame

Row	model	m	n	statistic	paper	replicated	dev
	String	Int64	Int64	String	Float64	Float64	Float64
1	DGP8-3	20	20	tau_hat	0.708	0.696	-0.012
2	DGP8-3	15	15	kappa_hat	0.457	0.446	-0.011
3	DGP8-2	15	15	tau_hat	0.422	0.412	-0.01
4	DGP8-2	10	10	kappa_hat	0.403	0.412	0.009
5	DGP8-3	20	20	kappa_tilde	0.154	0.145	-0.009

References

Weiß, Christian H, and Hee-Young Kim. 2024. “Using Spatial Ordinal Patterns for Non-Parametric Testing of Spatial Dependence.” Spatial Statistics 59: 100800.