Given an OTS of length \(T\) with range \(\mathcal{S}=\{s_0, s_1, s_2, \ldots, s_n\}\) (\(s_0 < s_1 < s_2 < \ldots < s_n\)),
\(\overline{X}_t=\{\overline{X}_1,\ldots, \overline{X}_T\}\), the function computes the
matrix \(\widehat{\boldsymbol P}(l) = \big(\widehat{p}_{i-1j-1}(l)\big)_{1 \le i, j \le n+1}\),
with \(\widehat{p}_{ij}(l)=\frac{N_{ij}(l)}{T-l}\), where \(N_{ij}(l)\) is the number
of pairs \((\overline{X}_t, \overline{X}_{t-l})=(s_i,s_j)\) in the realization \(\overline{X}_t\).