polars_online.corr

Correlation matrices, read and repaired (docs/ENHANCEMENTS.md E62).

polars_online.gram’s complement. Where that module solves and diagnoses a design matrix, this one takes a correlation matrix and does the arithmetic that comes after a pass: repair it to the nearest one that is a correlation matrix, shrink it towards a structured target, summarise it as an equicorrelation or a block matrix, place its eigenvalues against the Marchenko-Pastur edges, score a forecast of it, undo the Epps attenuation, and put a standard error on a single correlation.

Every function is a pure function of arrays. Where the input is a matrix it may equally be a mapping from gram() or a row from closed_groups(): matrix() reads all three.

The arithmetic is the papers’, named in each docstring, and every function is held against a longhand check in tests/test_corr.py – Higham’s own published examples for nearest(), Ledoit and Wolf’s formulae for shrink(), the closed forms for the rest.

Requires numpy, which is an optional extra of this package (pip install polars-online[numpy]) – not a dependency, as it is not one of polars’ either. Nothing here needs scipy or scikit-learn.

polars_online.corr.absorption(r: Any, k: int) → float

Kritzman, Li, Page and Rigobon’s absorption ratio: the share of total variance the top k eigenvectors explain, sum_{i<=k} lam_i / sum lam_i. High means the market is moving as one thing.

polars_online.corr.block_means(r: Any, labels: Sequence[Any]) → tuple[Any, Any]

The mean correlation within and between labelled blocks.

B[a, b] is the mean of R[i, j] over i in block a and j in block b, excluding the diagonal (a variable’s correlation with itself is 1 and says nothing about the block). Returns (B, counts) with the number of pairs behind each mean, in the order the labels first appear.

polars_online.corr.epps_invert(gram_or_row: Any, *, L: int) → Any

The correlation at scale L rows, from lagged co-moments at scale 1.

Toth and Kertesz’s equation 12: a correlation computed over fine intervals is attenuated because the two series do not move at the same instants, and the attenuation is undone by summing the lagged cross-covariances over the coarser interval. The weights are triangular, on the numerator and on both denominators:

rho_L[a, b] = sum_x (L - |x|) C_x[a, b]
              / sqrt(sum_x (L - |x|) C_x[a, a] * sum_x (L - |x|) C_x[b, b])

over x = -(L-1) .. L-1 with C_{-x}[a, b] = C_x[b, a] – so the numerator is L C_0[a, b] + sum_{l=1}^{L-1} (L - l) (C_l[a, b] + C_l[b, a]) and each auto term is L C_0[a, a] + 2 sum_l (L - l) C_l[a, a].

The input is a mapping from gram() or a closed row with lags and lag_comoments (ENHANCEMENTS E56), and must carry every lag 1 .. L-1: build it with ew_cov(lags=list(range(1, L))). A missing lag is an error naming it.

L = 1 is the plain correlation, which is the identity this reduces to.

polars_online.corr.equicorr(r: Any) → float

The mean off-diagonal correlation: the one number a deco tracks.

polars_online.corr.equicorr_loglik(row: Any, rho: float) → float

The Gaussian log-density of a standardised row under an equicorrelation matrix at rho, in closed form:

det R  = (1 - rho)**(n-1) * (1 + (n-1) rho)
r'R^-1 r = (S2 - rho S1**2 / (1 + (n-1) rho)) / (1 - rho)

which is deco’s loglik, offline. nan outside (-1/(n-1), 1), where R is not a correlation matrix.

polars_online.corr.equicorr_row(row: Any) → float

Engle and Kelly’s Lemma 2.3 on one standardised row:

u = (S1**2 - S2) / ((n - 1) * S2),  S1 = sum(r), S2 = sum(r*r)

the same closed form deco computes per row (ENHANCEMENTS E55), offline.

polars_online.corr.fisher_se(n: float, rho: float | None = None, phi_a: float | None = None, phi_b: float | None = None) → float

The standard error of a correlation estimated from n observations.

1 / sqrt(n - 3) is the standard error of Fisher’s z. With rho it is the delta-method error of the correlation itself, (1 - rho**2) / sqrt(n - 3).

With both phi it is inflated by sqrt((1 + phi_a phi_b) / (1 - phi_a phi_b)) for a pair of AR(1) series. That comes from Bartlett’s formula Var(r) ~ (1/n) sum_k rho_a(k) rho_b(k) and the geometric series 1 + 2 sum_{k>=1} (phi_a phi_b)**k. Both of its assumptions matter: it holds under a zero true cross-correlation and linear dependence, and it is not valid under ARCH-type innovations, where the variance of a sample correlation depends on the fourth moments and this understates it.

One phi without the other is an error: the inflation is a property of the pair.

polars_online.corr.from_blocks(b: Any, labels: Sequence[Any]) → Any

The block-equicorrelation matrix B describes: entry (i, j) is B[block(i), block(j)] off the diagonal, and 1 on it. The inverse of block_means() up to the within-block averaging, which is the test.

polars_online.corr.from_spectral(vals: Any, vecs: Any, *, unit_diag: bool = True) → Any

V' diag(vals) V, completed to a unit diagonal when unit_diag.

The remainder 1 - diag(V' L V) is non-negative, because the dropped components are PSD, so adding it to the diagonal leaves a correlation matrix rather than something that only looks like one.

polars_online.corr.from_z(z: Any) → Any

tanh(z), the inverse of to_z() away from the clip.

polars_online.corr.loss(fcst: Any, real: Any, kind: str = 'qlike', *, mu: Any = None, n: int | None = None) → float

How wrong a forecast correlation matrix was.

"qlike" is the Gaussian quasi-likelihood loss, shifted by the forecast-free constant so that it is zero at ``fcst == real`` and positive elsewhere:

tr(F^-1 R) - log det(F^-1 R) - k

"z_mse" is sum_{i<j} (z_ij(F) - z_ij(R))**2 * (n - 3), the squared Fisher-z error in units of its own sampling standard error; n is required.

"minvar" is Engle and Colacito’s minimum-variance loss w'Rw with w = F^-1 mu / (mu' F^-1 mu), the realised variance of the portfolio the forecast would have held; mu defaults to ones. It is minimised over F at F = R, which is what makes it a proper scoring rule for a covariance forecast.

polars_online.corr.matrix(obj: Any) → Any

The k x k correlation matrix of whatever this is.

An array (returned as a float array), a mapping from gram() or polars_online.gram.from_row() (its comoments scaled), or a one-row frame from closed_groups() (read through polars_online.gram.from_row() first).

polars_online.corr.mp_density(lam: Any, n: int, m: int, sigma2: float = 1.0) → Any

The Marchenko-Pastur density (Q / 2 pi sigma2) sqrt((lam_+ - lam) (lam - lam_-)) / lam, zero outside the edges – the curve to draw a spectrum against.

polars_online.corr.mp_edge(n: int, m: int, sigma2: float = 1.0) → tuple[float, float]

The Marchenko-Pastur edges for n observations of m series:

lam_pm = sigma2 * (1 +- 1/sqrt(Q))**2,    Q = n / m

(Laloux, Cizeau, Bouchaud and Potters 1999). Eigenvalues inside them are what pure noise produces; only those above lam_+ carry information. Q = 1 gives (0, 4 sigma2).

polars_online.corr.nearest(a: Any, w: Any = None, *, tol: float = 1e-08, max_iter: int = 100) → tuple[Any, float, int]

The nearest correlation matrix to a, by Higham (2002).

Alternating projections with Dykstra’s correction – his Algorithm 3.3, written out:

dS_0 = 0,  Y_0 = A
R_k  = Y_{k-1} - dS_{k-1}          # Dykstra's correction
X_k  = P_S(R_k)                    # onto the PSD cone
dS_k = X_k - R_k
Y_k  = P_U(X_k)                    # onto the unit diagonal

with P_U setting the diagonal to 1 and, for a diagonal weight w, P_S(A) = W^-1/2 (W^1/2 A W^1/2)_+ W^-1/2 where (.)_+ clips the eigenvalues at zero. It stops on his test (4.1): the largest of the relative infinity-norm changes in X, in Y, and between them, below tol.

Returns (X, dist, iters) with dist the weighted Frobenius distance ||W^1/2 (A - X) W^1/2||_F. Convergence is linear (a factor of about 3 an iteration on his own examples), so max_iter is a safety net rather than a knob; reaching it returns the last iterate.

The result has an exactly unit diagonal and is PSD to tol: the iteration converges to the boundary of the cone, so the smallest eigenvalue can be a small negative number of that order. Tighten tol if that matters; clipping it here would break the diagonal again.

a must be square and finite; it is symmetrised on the way in, since an “almost-correlation” matrix from two different estimates of the same pair is the case this exists for. w is the diagonal of W as a vector, all positive. max_iter must be at least 1: returning the input unprojected, with a distance of 0, would say it was already a correlation matrix.

polars_online.corr.shift(ar_fast: Any, ar_slow: Any, *, scale: float | None = None) → Any

The standardised absorption shift (fast - slow) / scale, elementwise over two aligned series of absorption ratios; scale defaults to the standard deviation of ar_slow over the sample. The two windows are the caller’s.

polars_online.corr.shrink(r: Any, target: str = 'constant', alpha: float | None = None, x: Any = None) → tuple[Any, float]

(1 - a) R + a F: Ledoit and Wolf (2004) shrinkage towards a structured target.

target="constant" (the default) is their constant-correlation target: f_ii = s_ii and f_ij = rbar * sqrt(s_ii s_jj) with rbar the mean off-diagonal correlation, which on a correlation matrix is the equicorrelation matrix at rbar. target="identity" is the identity.

alpha fixes the intensity. Left out, it is their optimal delta = max(0, min(kappa / T, 1)) with kappa = (pi - rho) / gamma, which needs the rows:

pi_ij    = (1/T) sum_t ((y_it - ybar_i)(y_jt - ybar_j) - s_ij)**2
theta_ii,ij = (1/T) sum_t ((y_it - ybar_i)**2 - s_ii)
                      * ((y_it - ybar_i)(y_jt - ybar_j) - s_ij)
rho      = sum_i pi_ii + sum_{i!=j} (rbar/2)
             * (sqrt(s_jj/s_ii) theta_ii,ij + sqrt(s_ii/s_jj) theta_jj,ij)
gamma    = sum_ij (f_ij - s_ij)**2

pi and theta are fourth-moment sums, so they cannot be recovered from R and T: pass x as the T x k sample the matrix came from (the standardised rows), or pass alpha yourself. One of the two is required.

``x`` has to be the sample ``r`` is of. s_ij above is the plain covariance of x, and gamma compares the target with it, so an r that is some other estimate of the same pairs – a decayed one out of a bank, say – makes kappa a comparison of two different matrices. The intensity is still in [0, 1] and the result still a shrunk matrix, but it is not the optimal intensity for the r handed in (docs/REVIEW-E54-E64.md K2). Pass alpha yourself when the two cannot be the same sample. x needs at least two rows.

Returns (shrunk, alpha). The result is positive definite whenever the target is and alpha > 0, which is the point: a sample correlation matrix from fewer rows than columns is singular, and this is the cheapest honest repair. nearest() is the other one, and they answer different questions – shrinkage trades bias for variance, Higham’s projection changes the matrix as little as possible.

polars_online.corr.signal_share(z_blocks: Any, n_eff_blocks: Any) → Any

How much of the between-block movement in a correlation is not sampling noise.

Per pair over B blocks, with z_b the Fisher-z of that block’s correlation and n_b its effective sample size:

clip(1 - mean_b(1 / (n_b - 3)) / var_b(z_b), 0, 1)

1 / (n - 3) is the sampling variance of z, so the ratio is the share of the observed variance the floor explains, and one minus it is what is left. 0 means the correlation moved no more than noise would. A 1-D input (one pair) returns a scalar.

polars_online.corr.spectral(r: Any, k: int) → tuple[Any, Any]

The top k eigenpairs, descending: (values, vectors) with vectors row-major k x n, each signed so its largest-magnitude entry is positive – the rule ew_cov’s pca uses on a first refresh, so the two agree up to that convention.

polars_online.corr.to_z(rho: Any) → Any

Fisher’s z = atanh(rho), elementwise, with |rho| clipped at Z_CLIP so a degenerate ±1 is finite (z ~ 7.25).