Tutorial for sequential ordinal pattern tests

Libraries

using StatsOrdinalPatterns, Random, Statistics, DataFrames, CairoMakie

Sequential testing with ordinal patterns

A classical test judges a complete series at once. A sequential test (control chart) instead monitors an incoming stream: after each new observation the ordinal pattern of the last \(m\) values is determined and its indicator vector \(\mathbf{e}_t\) (length \(m!\)) is smoothed exponentially (EWMA),

\[ \hat{\mathbf{p}}_t = \lambda\,\mathbf{e}_t + (1 - \lambda)\,\hat{\mathbf{p}}_{t-1}, \qquad \hat{\mathbf{p}}_0 = (1/m!, \ldots, 1/m!), \]

with smoothing parameter \(\lambda \in (0, 1]\). The chart statistic is computed from \(\hat{\mathbf{p}}_t\) and the chart signals as soon as the statistic leaves the critical limit cl. The number of observations until the first signal is the run length (RL).

Designing such a chart means fixing \(\lambda\) and then choosing cl so that the in-control average run length (ARL) matches a target L0 (typically 370 or 500). Chart performance is then judged by the out-of-control ARL: the smaller, the faster dependence is detected. See Weiß and Testik (2023).

Key functions

• Critical limit
  • cl_op(op_dgp, lam, L0, cl_init; chart_choice, m=3, d=1, reps_bracket, reps_final, bracket_step, seed)
    • Returns the cl for which the in-control ARL equals L0.
    • op_dgp: in-control process, ContinuousDGPIC(dist) or DiscreteDGPIC(dist, add_noise).
    • cl_init: rough starting value; the search first brackets the root with reps_bracket replications in steps of bracket_step, then refines it with reps_final replications.
    • seed: fix an integer for reproducible results.
• Chart performance
  • arl_op_ic(op_dgp, lam, cl, reps; chart_choice, m=3, d=1) — in-control ARL, returns (ARL, standard error).
  • arl_op_oc(op_dgp, lam, cl, reps; chart_choice, m=3, d=1) — out-of-control ARL for a dependent process: AR1, MA1, MA2, TEAR1, AAR1, QAR1 for continuous data, INAR1, DAR1, BAR1, … for counts.
• Monitoring observed data
  • stat_op(data, lam; chart_choice, m=3, d=1) — returns (stats, p_rel), where stats holds the sequentially computed EWMA chart statistics of the series.
  • monitor_op(data, lam, cl; chart_choice, m=3, d=1) — applies the alarm rule of the chart to these statistics and returns a ControlChartResult with the statistics, the limit and the first alarm. With Makie loaded, plot draws it.

Alarm rules

The same seven statistics as in the classical tests are available; the alarm rule follows the direction of the statistic.

Statistic	Constructor	Under \(H_0\)	Signal when
Persistence (\(\hat{\tau}\))	Persistence()	\(\hat{\tau} = 0\)	\(\lvert\hat{\tau}\rvert > cl\)
UpDownBalance (\(\hat{\beta}\))	UpDownBalance()	\(\hat{\beta} = 0\)	\(\lvert\hat{\beta}\rvert > cl\)
RotationalAsymmetry (\(\hat{\gamma}\))	RotationalAsymmetry()	\(\hat{\gamma} = 0\)	\(\lvert\hat{\gamma}\rvert > cl\)
UpDownScaling (\(\hat{\delta}\))	UpDownScaling()	\(\hat{\delta} = 0\)	\(\lvert\hat{\delta}\rvert > cl\)
DistanceToWhiteNoise (\(\hat{\Delta}\))	DistanceToWhiteNoise()	\(\hat{\Delta} = 0\)	\(\hat{\Delta} > cl\)
Shannon entropy (\(\hat{H}\))	Shannon(base=exp(1))	\(\hat{H} = \log(m!)\)	\(\hat{H} < cl\)
Shannon extropy (\(\hat{H}_\text{ex}\))	ShannonExtropy(base=exp(1))	\(\hat{H}_\text{ex} = 5\log(6/5)\)	\(\hat{H}_\text{ex} < cl\)

cl_op handles the direction internally: for the two entropy charts the ARL decreases in cl, for all others it increases.

Setup

ic = ContinuousDGPIC(Normal(0, 1))   # in-control: i.i.d. standard normal
lam = 0.1                            # EWMA smoothing parameter
L0 = 370.0                           # target in-control ARL
m, d = 3, 1

(3, 1)

Computing the critical limit

cl = cl_op(
    ic, lam, L0, 0.20;
    chart_choice=Persistence(), m=m, d=d,
    reps_bracket=1_000, reps_final=5_000, bracket_step=0.02, seed=42,
    verbose=true
)

============================================================
 STEP 1: Bracket Search
 reps = 1000 | truncation cap = 18500 | step size = 0.02
============================================================
  cl = 0.2 | ARL = 83.176
  cl = 0.22 | ARL = 142.344
  cl = 0.24 | ARL = 245.676
  cl = 0.26 | ARL = 479.155
------------------------------------------------------------
 Bracket found: [0.2, 0.26]
------------------------------------------------------------

============================================================
 STEP 2: Root Finding via ITP
 reps = 5000 | no truncation | tolerance = 0.0001
============================================================
  cl = 0.2 | ARL = 83.753
  cl = 0.26 | ARL = 467.8048
  cl = 0.244 | ARL = 283.64
  cl = 0.251554 | ARL = 359.2096
  cl = 0.252408 | ARL = 367.785
  cl = 0.252587 | ARL = 369.6772
  cl = 0.2536 | ARL = 381.9754
  cl = 0.252614 | ARL = 369.7836
  cl = 0.252632 | ARL = 369.7836
  cl = 0.2528 | ARL = 371.0074
  cl = 0.252661 | ARL = 369.9982
  cl = 0.252662 | ARL = 369.9982
  cl = 0.2527 | ARL = 370.2616

Results of univariate zero finding:

* Converged to: 0.25266172651450286
* Algorithm: Roots.ITP{Float64, Int64}(0.2, 2, 1)
* iterations: 11
* function evaluations ≈ 13
* stopped as x_n ≈ x_{n-1} using atol=xatol, rtol=xrtol

Trace:
(a₀, b₀) = ( 0.20000000000000001, 0.26000000000000001 )
(a₁, b₁) = ( 0.24400006120007772, 0.26000000000000001 )
(a₂, b₂) = ( 0.25155407738768609, 0.26000000000000001 )
(a₃, b₃) = ( 0.25240756060663316, 0.26000000000000001 )
(a₄, b₄) = ( 0.25258722887481488, 0.26000000000000001 )
(a₅, b₅) = ( 0.25258722887481488, 0.25359999999999999 )
(a₆, b₆) = ( 0.25261401697267649, 0.25359999999999999 )
(a₇, b₇) = ( 0.2526317122440343, 0.25359999999999999 )
(a₈, b₈) = ( 0.2526317122440343, 0.25279999999999997 )
(a₉, b₉) = ( 0.25266147560584368, 0.25279999999999997 )
(a₁₀, b₁₀) = ( 0.25266172651450286, 0.25279999999999997 )
(a₁₁, b₁₁) = ( 0.25266172651450286, 0.25269999999999998 )

0.25266172651450286

verbose=true prints cl and the associated ARL at every evaluation, so both phases of the search can be followed: the bracketing steps away from cl_init until the ARL crosses L0, then the ITP algorithm refines the bracket. It is worth checking that the bracketing does not crawl — cl_init = 0.20 only has to be in the right ballpark, but it should sit below the root and bracket_step should be of the order of the remaining distance.

Increase reps_final for a more precise limit — the accuracy of cl is limited by the Monte Carlo error of the ARL estimate, not by the root finder. Pass verbose=false (the default) to get the limit alone.

Verifying the in-control ARL

arl_op_ic(ic, lam, cl, 10_000; chart_choice=Persistence(), m=m, d=d)

(365.2791, 3.552213651021184)

The returned tuple is (ARL, standard error); the ARL is within a few standard errors of L0.

Out-of-control performance

The same cl is used to monitor dependent processes. The stronger the dependence, the faster the chart signals.

oc_dgps = [
    ("AR(1), α = 0.1", AR1(0.1, Normal(0, 1))),
    ("AR(1), α = 0.25", AR1(0.25, Normal(0, 1))),
    ("AR(1), α = 0.5", AR1(0.5, Normal(0, 1))),
    ("MA(1), α = 0.5", MA1(0.5, Normal(0, 1))),
    ("TEAR(1), α = 0.5", TEAR1(0.5, Normal(0, 1))),
    ("INAR(1), α = 0.5", INAR1(0.5, Poisson(5), true)),
]
oc_res = [
    arl_op_oc(dgp, lam, cl, 10_000; chart_choice=Persistence(), m=m, d=d)
    for (_, dgp) in oc_dgps
]

DataFrame(
    Process=[n for (n, _) in oc_dgps],
    ARL=round.([r[1] for r in oc_res], digits=2),
    se=round.([r[2] for r in oc_res], digits=3),
)

6×3 DataFrame

Row	Process	ARL	se
	String	Float64	Float64
1	AR(1), α = 0.1	268.5	2.593
2	AR(1), α = 0.25	160.56	1.501
3	AR(1), α = 0.5	76.74	0.693
4	MA(1), α = 0.5	259.28	2.521
5	TEAR(1), α = 0.5	73.39	0.658
6	INAR(1), α = 0.5	77.79	0.701

The count process INAR1(α, dist, add_noise) uses add_noise=true, which breaks ties by adding uniform noise — otherwise equal values would make the ordinal pattern ambiguous.

Comparing charts

Each chart needs its own critical limit, because the statistics live on different scales. Once all charts are calibrated to the same L0, their out-of-control ARLs are directly comparable.

chart_list = [
    ("Persistence (τ)", Persistence(), 0.20, 0.02),
    ("UpDownBalance (β)", UpDownBalance(), 0.30, 0.03),
    ("RotationalAsymmetry (γ)", RotationalAsymmetry(), 0.30, 0.03),
    ("DistanceToWhiteNoise (Δ)", DistanceToWhiteNoise(), 0.08, 0.01),
    ("Shannon entropy (H)", Shannon(base=exp(1)), 1.40, 0.02),
]

oc_dgp = AR1(0.25, Normal(0, 1))
cls = [
    cl_op(ic, lam, L0, init; chart_choice=c, m=m, d=d,
        reps_bracket=1_000, reps_final=5_000, bracket_step=step, seed=42)
    for (_, c, init, step) in chart_list
]
ic_arls = [
    arl_op_ic(ic, lam, cls[i], 10_000; chart_choice=chart_list[i][2], m=m, d=d)[1]
    for i in eachindex(chart_list)
]
oc_arls = [
    arl_op_oc(oc_dgp, lam, cls[i], 10_000; chart_choice=chart_list[i][2], m=m, d=d)[1]
    for i in eachindex(chart_list)
]

DataFrame(
    Chart=[t[1] for t in chart_list],
    cl=round.(cls, digits=4),
    ARL_ic=round.(ic_arls, digits=1),
    ARL_oc=round.(oc_arls, digits=1),
)

5×4 DataFrame

Row	Chart	cl	ARL_ic	ARL_oc
	String	Float64	Float64	Float64
1	Persistence (τ)	0.2527	365.3	158.0
2	UpDownBalance (β)	0.3623	361.5	251.3
3	RotationalAsymmetry (γ)	0.4057	367.5	413.7
4	DistanceToWhiteNoise (Δ)	0.1111	365.9	210.1
5	Shannon entropy (H)	1.4614	357.2	269.1

All charts hold the in-control ARL near L0 = 370, so the ARL_oc column ranks them: the smaller the value, the faster this particular AR(1) dependence is detected. The ranking is specific to the alternative — \(\hat{\gamma}\) measures rotational asymmetry, which an AR(1) process does not induce, so its chart stays at roughly the in-control level.

Monitoring an observed series

monitor_op(data, lam, cl; ...) computes the EWMA chart statistics of an observed series (the same ones stat_op(data, lam; ...) returns), applies the alarm rule of the chosen chart and records the first signal. Here the series switches from i.i.d. to AR(1) at \(t = 150\).

Random.seed!(4)
t_change = 150
x = randn(400)
for t in (t_change+1):400
    x[t] = 0.5 * x[t-1] + randn()
end

res = monitor_op(x, lam, cl; chart_choice=Persistence(), m=m, d=d)

ControlChartResult
  Family:           ordinal patterns (OP)
  Chart:            Persistence()
  λ:                0.1
  Control limit:    0.2527
  Alarm rule:       |stat| > cl
  Statistics:       398 (observations 3 to 400)
  First alarm:      observation 178 (statistic 0.2616)

The statistic at position \(i\) of res.stats is computed from the pattern covering the observations \(i, \ldots, i + (m-1)d\), so it becomes available at the last of them. This observation index is stored in res.time, and res.alarm is the position of the first signal in res.stats:

(first_alarm=res.time[res.alarm], delay=res.time[res.alarm] - t_change,
 stat_at_alarm=round(res.stats[res.alarm], digits=4))

(first_alarm = 178, delay = 28, stat_at_alarm = 0.2616)

Here the chart passes the i.i.d. part without a false alarm and signals shortly after the change. A different series may well signal earlier: with L0 = 370, a false alarm occurs on average once every 370 in-control observations.

With CairoMakie loaded, plot shows the statistics against the control limits, shades the rejection region and marks the first alarm. See the plotting tutorial for the options.

plot(res)

References

Weiß, Christian H, and Murat Caner Testik. 2023. “Nonparametric Control Charts for Monitoring Serial Dependence Based on Ordinal Patterns.” Technometrics 65 (3): 340–50.