Local Principal Component Analysis (LPCA)

When to Use This Estimator

LPCA implements the local principal component analysis of Feng [LPCA]. It answers a question the rest of the factor-model family cannot: what if the outcome responds to the latent drivers nonlinearly?

Every factor estimator in mlsynth – Generalized Synthetic Control (GSYNTH), Causal Factor Model (CFM), Factor Model Approach (FMA), Matrix Completion with Nuclear Norm Minimization (MCNNM), Robust Matrix estimation with Side Information (RMSI) – writes the untreated outcome as \(y_{jt}(0) = \mathbf{f}_t^{\top}\boldsymbol{\alpha}_j + u_{jt}\). Each unit’s path is a linear combination of a few common factors, so the outcome matrix is low rank and can be recovered by principal components or by shrinking its singular values. Drop the linearity and that machinery stops working. If

\[y_{jt}(0) = \eta_t(\boldsymbol{\alpha}_j) + u_{jt},\]

with \(\boldsymbol{\alpha}_j\) low-dimensional but \(\eta_t\) an unknown, possibly nonlinear function, then the matrix is generally full rank even though the structure driving it has only a few dimensions. Global principal components and nuclear-norm completion both misread such a panel.

The way through is that a smooth surface looks linear up close. Two units with similar \(\boldsymbol{\alpha}\) have outcomes that a first-order expansion makes an approximately linear factor model, so principal components apply within a neighbourhood even where they fail globally. Since \(\boldsymbol{\alpha}_j\) is never observed, the neighbourhood has to be built by matching on the outcomes themselves.

Reach for LPCA when the panel is wide in both dimensions – many units and many periods – and the response of outcomes to the latent drivers is plausibly curved: saturating demand, threshold effects, bounded shares, anything where doubling a latent input does not double the outcome. It is a matrix-completion estimator in the same family as Matrix Completion with Nuclear Norm Minimization (MCNNM) and Synthetic Nearest Neighbors / Causal Matrix Completion (SNN), and like them it returns an imputation, not a donor-weight story.

Do not use LPCA when

• The panel is small. The neighbourhood defaults to \(K = \operatorname{round}(N^{2/3})\) and half the periods are spent on matching, so a classic thirty-donor, thirty-period study leaves nine neighbours and fifteen fitting periods. Use Vanilla Synthetic Control (VanillaSC), Forward Difference-in-Differences (FDID) or Matrix Completion with Nuclear Norm Minimization (MCNNM) there.

• A linear factor structure is credible. LPCA pays for its generality in variance; when the matrix really is low rank, Generalized Synthetic Control (GSYNTH) or Matrix Completion with Nuclear Norm Minimization (MCNNM) estimate the same object more efficiently.

• Inference is the deliverable. The paper reports no standard errors, confidence intervals or p-values – see Inference below. Use Factor Model Approach (FMA) or Counterfactual Synthetic Control with Instrumented PCA (CSC-IPCA) when a band is required.

• The post-treatment window is long. The theory assumes the number of treated periods is fixed while the panel grows; see Assumption 4.

• Interpretable weights are the result. Like other completion estimators, LPCA hands back a fitted path. Use Two-Step Synthetic Control or Synthetic Control with Multiple Outcomes (SCMO).

Notation

The outcome panel is the \(N \times T\) matrix \(\mathbf{Y} = (y_{jt})\) over units \(j \in \mathcal{N} \coloneqq \{1, \ldots, N\}\) and periods \(t \in \mathcal{T} \coloneqq \{1, \ldots, T\}\). Unit 1 is treated from \(T_0 + 1\), splitting \(\mathcal{T}\) into \(\mathcal{T}_1 \coloneqq \{t \in \mathcal{T} : t \le T_0\}\) and \(\mathcal{T}_2 \coloneqq \{t \in \mathcal{T} : t > T_0\}\), with \(T_2 \coloneqq |\mathcal{T}_2|\).

Each unit carries a latent \(\boldsymbol{\alpha}_j \in \mathbb{R}^r\) with \(r\) small, and each period a latent response function \(\eta_t\). The untreated outcome is \(y_{jt}(0) = \eta_t(\boldsymbol{\alpha}_j) + u_{jt}\) with mean-zero noise \(u_{jt}\).

LPCA needs a second split, this time of the periods, and it is not the treatment split. Write \(\mathcal{T} = \mathcal{M} \cup \mathcal{P}\) where \(\mathcal{M} \coloneqq \{1, \ldots, M\}\) is the matching block and \(\mathcal{P} \coloneqq \mathcal{T} \setminus \mathcal{M}\) the PCA block. Matching happens on \(\mathcal{M}\), fitting on \(\mathcal{P}\), and the intervention must fall inside \(\mathcal{P}\) – that is, \(M < T_0\) – so the fit has pre-treatment periods to be judged on.

The neighbourhood of unit \(j\) is \(\mathcal{K}_j \subset \mathcal{N}\) with \(|\mathcal{K}_j| = K\), containing \(j\) itself.

The estimator

Step one, matching. On the matching block, measure how far apart two units are by the pseudo-max distance of Zhang, Levina and Zhu [ZLZ2017],

\[\rho(\mathbf{y}_i, \mathbf{y}_j) \coloneqq \max_{l \neq i, j} \Bigl| \frac{1}{M} \sum_{t \in \mathcal{M}} (y_{it} - y_{jt})\, y_{lt} \Bigr| ,\]

and let \(\mathcal{K}_j\) collect the \(K\) smallest. Averaging across the other units is what does the work: the idiosyncratic noise washes out as \(M\) grows, so the distance reflects the noise-free structure \(\eta_t(\boldsymbol{\alpha})\) and units close in it are close in \(\boldsymbol{\alpha}\). The Euclidean distance does not have this property under heteroskedastic errors, which is why it is not the default.

Step two, local principal components. Stack the neighbourhood’s rows over the PCA block into \(\mathbf{Y}_{\langle j \rangle} \in \mathbb{R}^{K \times |\mathcal{P}|}\) and take its rank-\(d\) truncated singular value decomposition,

\[(\widehat{\mathbf{F}}_{\langle j \rangle}, \widehat{\boldsymbol{\Lambda}}_{\langle j \rangle}) \coloneqq \operatorname*{argmin}_{\mathbf{F}, \boldsymbol{\Lambda}} \bigl\| \mathbf{Y}_{\langle j \rangle} - \mathbf{F} \boldsymbol{\Lambda}^{\top} \bigr\|_F^2 ,\]

and read the counterfactual off unit \(j\)’s row of \(\widehat{\mathbf{F}}_{\langle j \rangle} \widehat{\boldsymbol{\Lambda}}_{\langle j \rangle}^{\top}\). The treatment effect is \(\tau_t \coloneqq y_{1t} - \widehat{y}_{1t}(0)\) for \(t \in \mathcal{T}_2\) and the ATT is \(\widehat{\tau} \coloneqq T_2^{-1} \sum_{t \in \mathcal{T}_2} \tau_t\).

Assumptions and remarks

Assumption 1 (smooth low-dimensional latent structure). The latent variables \(\boldsymbol{\alpha}_j\) are i.i.d. on a compact support and each \(\eta_t\) is smooth, with bounded derivatives up to some order (paper Assumption 2.1). Remark. Smoothness is what makes the tangent-plane approximation good; the dimension \(r\) being small is what makes neighbours exist. Neither requires the matrix to be low rank, which is the whole point.

Assumption 2 (informative matching). The chosen distance denoises the data, and distances in the noise-free structure translate into distances in \(\boldsymbol{\alpha}\) (paper Assumption 4.1). Remark. Two things can break here and they are different. If the panel is short, averaging does not kill the noise and the neighbours are wrong. If two distinct \(\boldsymbol{\alpha}\) values generate near-identical outcome paths, the observables are uninformative and no distance recovers them. The second is a completeness condition, not a sample-size problem.

Assumption 3 (independent blocks). Matching and fitting use disjoint periods, and the noise is independent (or weakly dependent) across periods (paper Section 3). Remark. This is the reason for the split, and it is not cosmetic. Choosing neighbours uses the realised noise; fitting on those same periods would correlate the estimated factors with the errors and standard PCA would no longer apply.

Assumption 4 (short treated window). The number of treated periods is fixed as the panel grows (paper Theorem 6.1). Remark. The treated cells are set to their period mean before fitting, and the theorem bounds the damage that does. The bound is only meaningful when few cells are affected: Feng’s Kansas application zeroes 16 cells of 104. A long post-period puts the treated unit’s own imputed values into the decomposition that is supposed to predict them.

Preprocessing is yours to do

The estimator takes the outcome column as given. Feng’s application works with first-differenced log GDP per capita, not levels, and that transformation is the user’s to apply – transform the column before passing it in. The counterfactual comes back on the scale of whatever outcome names.

Internally each period is centred across units before matching, and the period means are added back to the counterfactual. That round trip has to close: a counterfactual left in centred space is offset from the observed series by the average of the period means, which is a difference of the same order as the effect being measured. It is pinned by a test.

Inference

There is none. That is the paper’s position, not an omission here. Theorem 6.1 gives a uniform max-norm convergence rate for the fitted surface; Section 6.1 reports point predictions, and its figures carry no bands. LPCA therefore returns an empty inference slot and res.att_ci is None.

Four diagnostics stand in for it, all on the result. neighbourhood_size is the realised neighbourhood, which can exceed the requested \(K\) because the selection rule is a threshold and ties are all kept – discrete or binary panels tie routinely. local_rank is the number of components the singular-value ratio rule retained.

neighbour_weights gives the weight each neighbour receives. The reconstruction is literally that weighted combination of the neighbourhood’s rows: the weights are a column of the rank-\(d\) projector \(\mathbf{U}_d \mathbf{U}_d^{\top}\), so unlike a synthetic control’s they neither sum to one nor stay non-negative. They are still readable as a comparison set.

self_weight is the one to check. The treated unit belongs to its own neighbourhood, so its own row enters the decomposition that produces its counterfactual – and for the post-treatment periods that row holds the values the estimator overwrote with the period mean. This number says how much the counterfactual leans on them, and it is the practical face of Assumption 4. It lies in \([0, 1]\), and across a neighbourhood the self-weights sum to the rank, so \(d / K\) is the natural benchmark. Well above that, on a long post-period, is the configuration Theorem 6.1 stops covering.

Warning

The rank rule is inert for \(K \le 15\). It keeps components while consecutive singular values satisfy \(\sigma_i / \sigma_{i+1} < \log \log K\). Ratios of a descending spectrum are at least 1, and \(\log \log K < 1\) for \(K \le 15\) (the threshold crosses 1 at \(e^e \approx 15.15\)), so below sixteen neighbours the comparison never fires and the rank is pinned at max_components - 1 whatever the data says. Feng’s Kansas application uses \(K = 14\), so its rank of 2 is mechanical. Check res.metadata["rank_rule_active"]. Paper Remark 4.3 leaves formal rank selection to future research.

The reported window

Local PCA predicts only the periods held out from matching, so res.counterfactual and res.time_series cover \(\mathcal{P}\), not all of \(\mathcal{T}\). The first \(M\) periods are dropped from the reported series instead of being padded, and time_series.time_periods carries the labels that remain. pre_rmse is therefore computed on \(\mathcal{P} \cap \mathcal{T}_1\).

Example

The 2012 Kansas tax cuts. Feng analyses quarterly log GDP per capita growth for the 50 states, with Kansas treated from 2012Q2, so the outcome is first-differenced before it reaches the estimator.

import pandas as pd
from mlsynth import LPCA

url = ("https://raw.githubusercontent.com/jgreathouse9/mlsynth/"
       "refs/heads/main/basedata/kansas_taxcut.csv")
df = pd.read_csv(url).sort_values(["state", "year_qtr"])

# The paper works in growth rates, in percent.
df["growth"] = df.groupby("state")["lngdpcapita"].diff() * 100.0
df = df.dropna(subset=["growth"])

res = LPCA({
    "df": df, "outcome": "growth", "treat": "treated",
    "unitid": "state", "time": "year_qtr",
    "match_periods": 40,        # the paper's split
    "display_graphs": True,
}).fit()

print(f"ATT (2012Q2-2016Q1) = {res.att:+.4f} pp")
print(f"K = {res.n_neighbours}, realised = {res.neighbourhood_size}")
print(f"local rank = {res.local_rank} "
      f"(rule active: {res.metadata['rank_rule_active']})")
print(f"neighbours: {', '.join(res.neighbours)}")
print(f"self weight = {res.self_weight:.3f} "
      f"(benchmark {res.local_rank / res.n_neighbours:.3f})")

The counterfactual sits about 0.53 points above observed Kansas growth: the paper’s estimate that the tax cut cost growth.

Verification

Note

Empirical (Kansas). The ATT reproduces Feng’s Section 6.1 to four decimals (\(-0.5306\) against the reported \(-0.53\)), and the observed series falls below the LPCA path in 9 of 16 post-treatment quarters, as reported. Pinned in mlsynth/tests/test_lpca.py.

Monte Carlo. All 48 cells of the paper’s Table 1 reproduce at 500 replications, median disagreement 0.83 Monte Carlo standard errors, with local PCA beating global PCA on the two nonlinear designs and the advantage widening with the severity of the nonlinearity.

Both live in benchmarks/reference/lpca_kansas/, which also records that the synthetic-control comparison in the paper’s November 2023 version was a centring defect the author corrected in July 2024. The docs above quote the corrected version.

Core API

LPCA: Local Principal Component Analysis (Feng 2024).

Feng, Y. (2024). “Optimal Estimation of Large-Dimensional Nonlinear Factor Models.” Working paper; arXiv:2311.07243.

Every factor-model estimator in mlsynth writes the untreated outcome as \(y_{it}(0) = f_t' \alpha_i + u_{it}\) – the outcome matrix is low-rank because each unit’s path is a linear combination of a few common factors. LPCA drops the linearity. It supposes only

\[y_{it}(0) = \eta_t(\alpha_i) + u_{it},\]

where the latent variable \(\alpha_i\) is low-dimensional but the response \(\eta_t\) is an unknown, possibly nonlinear function. Such a matrix is generally full rank, so nuclear-norm completion and principal components fitted to the whole panel both misread it.

The way through is that a smooth surface looks linear up close. If two units have similar \(\alpha\), a first-order expansion makes their outcomes an approximately linear factor model, so PCA applies within a neighbourhood even though it fails globally. Since \(\alpha_i\) is unobserved, the neighbourhood is built by matching on the outcomes themselves:

1. Split the periods in two. On the first block, compute the pseudo-max distance between units and take each unit’s K nearest neighbours. The split matters: searching for neighbours uses the noise, so matching and fitting on the same periods would correlate the factors with the errors.

2. On the second block, take a truncated SVD of the neighbour submatrix and read off the treated unit’s row. That row is the estimate of \(\eta_t(\alpha_i)\) – the counterfactual.

The treated unit’s post-treatment cells are set to the period mean before any of this, and Theorem 6.1 is the statement that doing so perturbs the estimate by a vanishing amount when the post-period is short relative to the panel.

LPCA wants a wide panel: neighbours are only close in \(\alpha\) when many units are available to choose among, and the paper’s theory is asymptotic in both dimensions. It is the wrong tool for a classic thirty-donor study, where K = round(N^(2/3)) leaves a handful of neighbours and half the periods go to matching.

The paper reports no standard errors, confidence intervals or p-values, and Theorem 6.1 is a uniform max-norm rate. This estimator therefore returns a point counterfactual and no inference.

class mlsynth.estimators.lpca.LPCA(config: LPCAConfig | dict)

Bases: object

Local Principal Component Analysis estimator.

Parameters:

config (LPCAConfig or dict) – Configuration object. See mlsynth.config_models.LPCAConfig.

fit() → LPCAResults

Run local PCA and return LPCAResults.

Configuration

class mlsynth.config_models.LPCAConfig(*, df: ~pandas.DataFrame, outcome: str, treat: str, unitid: str, time: str, display_graphs: bool = True, save: bool | str = False, counterfactual_color: ~typing.List[str] = <factory>, treated_color: str = 'black', plot: ~mlsynth.config_models.PlotConfig = <factory>, n_neighbors: ~typing.Annotated[int | None, ~annotated_types.Ge(ge=2)] = None, match_periods: ~typing.Annotated[int | None, ~annotated_types.Ge(ge=1)] = None, match_fraction: ~typing.Annotated[float, ~annotated_types.Gt(gt=0.0), ~annotated_types.Lt(lt=1.0)] = 0.5, distance: ~typing.Literal['pseudo_max', 'euclidean', 'average'] = 'pseudo_max', max_components: ~typing.Annotated[int, ~annotated_types.Ge(ge=1)] = 3, demean: bool = True)

Configuration for the LPCA estimator.

Feng (2024), “Optimal Estimation of Large-Dimensional Nonlinear Factor Models” (arXiv:2311.07243). Imputes the treated unit’s counterfactual under a factor structure that may be nonlinear: units are matched to their K nearest neighbours on one block of periods, and a truncated SVD of the neighbour submatrix on the complementary block supplies the fit. Inherits the standard df / outcome / treat / unitid / time interface.

The estimator takes the outcome as given. Feng’s Kansas application works with first-differenced log GDP per capita; transform the outcome before passing it in, so the returned counterfactual is always on the scale of the column named by outcome.

Parameters:

• n_neighbors (int, optional) – Neighbourhood size K, counting the unit itself. Defaults to round(N ** (2/3)), the paper’s rule. Ties in the distance are all kept, so the realised neighbourhood can exceed this; the result reports both.

• match_periods (int, optional) – Number of leading periods reserved for nearest-neighbour matching. Defaults to round(T * match_fraction). Feng’s Kansas application uses 40 of 104 quarters. Must leave at least one pre-treatment period in the remaining block.

• match_fraction (float) – Fraction of periods used for matching when match_periods is None. Default 0.5, matching the paper’s simulations.

• distance ({“pseudo_max”, “euclidean”, “average”}) – Distance used for matching. The default separates units by their noise-free factor structure under heteroskedastic errors, which the Euclidean distance does not.

• max_components (int) – Maximum number of local components extracted (the paper’s nlam, default 3). The realised rank is chosen by a singular-value ratio rule and reported on the result.

• demean (bool) – Centre each period across units before matching, adding the period means back to the counterfactual. Default True.

demean: bool

distance: Literal['pseudo_max', 'euclidean', 'average']

match_fraction: float

match_periods: int | None

max_components: int

model_config: ClassVar[ConfigDict] = {'arbitrary_types_allowed': True, 'extra': 'forbid'}

Configuration for the model, should be a dictionary conforming to [ConfigDict][pydantic.config.ConfigDict].

n_neighbors: int | None

Result Containers

LPCA.fit() returns a LPCAResults: the ATT, the counterfactual path over the PCA block, the matched neighbours and their realised count, the neighbour_weights and the treated unit’s self_weight, the local_rank and the local singular_values, the period_means removed and restored, and the standardized sub-models.

class mlsynth.utils.lpca_helpers.structures.LPCAResults(*, effects: ~mlsynth.config_models.EffectsResults | None = None, fit_diagnostics: ~mlsynth.config_models.FitDiagnosticsResults | None = None, time_series: ~mlsynth.config_models.TimeSeriesResults | None = None, weights: ~mlsynth.config_models.WeightsResults | None = None, inference: ~mlsynth.config_models.InferenceResults | None = None, method_details: ~mlsynth.config_models.MethodDetailsResults | None = None, sub_method_results: ~typing.Dict[str, ~typing.Any] | None = None, additional_outputs: ~typing.Dict[str, ~typing.Any] | None = None, raw_results: ~typing.Dict[str, ~typing.Any] | None = None, execution_summary: ~typing.Dict[str, ~typing.Any] | None = None, plot_config: ~mlsynth.config_models.PlotConfig | None = None, inputs: ~mlsynth.utils.lpca_helpers.structures.LPCAInputs, counterfactual_matrix: ~numpy.ndarray, neighbours: ~typing.List[~typing.Any], neighbours_by_unit: ~typing.Dict[~typing.Any, ~typing.List[~typing.Any]], neighbour_weights: ~typing.Dict[~typing.Any, float], self_weight: float, n_neighbours: int, neighbourhood_size: int, local_rank: int, singular_values: ~numpy.ndarray, period_means: ~numpy.ndarray | None = None, metadata: ~typing.Dict[str, ~typing.Any] = <factory>)

Bases: BaseEstimatorResults

Top-level container returned by mlsynth.LPCA.fit().

An EffectResult: alongside the LPCA-specific fields below it exposes the standardized sub-models (effects, time_series, weights, inference, fit_diagnostics, method_details) and the flat accessors att / counterfactual / gap / att_ci / pre_rmse.

The reported series cover the PCA block only. Local PCA predicts the periods held out from matching, so the first match_periods periods have no counterfactual and are excluded from time_series instead of padded; time_series.time_periods carries the labels that remain.

Parameters:

• inputs (LPCAInputs)

• counterfactual_matrix (np.ndarray) – Reconstructed untreated outcome for the treated units on the PCA block, shape (n_treated, T - match_periods), on the outcome’s scale.

• neighbours (list) – Names of the treated unit’s matched neighbours, itself included. With several treated units these describe the first; the rest are in neighbours_by_unit.

• neighbours_by_unit (dict) – {treated unit name: list of neighbour names}.

• neighbour_weights (dict) – {neighbour name: weight} for the treated unit. The reconstruction is literally this weighted combination of the neighbourhood’s rows. Being a projection and not a simplex fit, the weights neither sum to one nor stay non-negative.

• self_weight (float) – The treated unit’s own weight in its reconstruction. The treated row belongs to its own neighbourhood, so this measures how much the counterfactual leans on the post-treatment cells the estimator set to the period mean – the perturbation Theorem 6.1 bounds. A large value on a long post-period is the warning sign.

• n_neighbours (int) – Requested neighbourhood size K.

• neighbourhood_size (int) – Realised size for the treated unit. Exceeds n_neighbours when the distance ties, which discrete panels do routinely.

• local_rank (int) – Components kept by the singular-value ratio rule. Pinned at max_components - 1 whenever n_neighbours <= 15, where the rule cannot fire – see select_rank().

• singular_values (np.ndarray) – Leading singular values of the treated unit’s neighbourhood.

• period_means (np.ndarray, optional) – Per-period means removed before matching and added back to the counterfactual; None when demean is False.

• metadata (dict)

counterfactual_matrix: np.ndarray

inputs: LPCAInputs

local_rank: int

metadata: Dict[str, Any]

model_config: ClassVar[ConfigDict] = {'arbitrary_types_allowed': True, 'extra': 'forbid', 'frozen': True, 'json_encoders': {<class 'numpy.ndarray'>: <function BaseEstimatorResults.Config.<lambda>>}}

Configuration for the model, should be a dictionary conforming to [ConfigDict][pydantic.config.ConfigDict].

n_neighbours: int

neighbour_weights: Dict[Any, float]

neighbourhood_size: int

neighbours: List[Any]

neighbours_by_unit: Dict[Any, List[Any]]

period_means: np.ndarray | None

self_weight: float

singular_values: np.ndarray

Helper Modules

The two numerical steps: the matching distances, the neighbourhood rule, the rank rule, and the local reconstruction.

The two numerical steps of local principal component analysis.

Feng (2024), Optimal Estimation of Large-Dimensional Nonlinear Factor Models, Section 3, Algorithm 1: match each unit to its K nearest neighbours on one block of periods, then take a truncated SVD of the neighbour submatrix on the complementary block. Ported from the author’s knn.index / findknn / lpca / sel.r.

The matrix convention here is units by periods, the transpose of the paper’s X.

mlsynth.utils.lpca_helpers.core.local_pca_row(block: ndarray, neighbours: ndarray, unit: int, n_neighbours: int, n_components: int = 3) → tuple[ndarray, int]

Local PCA reconstruction of one unit’s row on the PCA block.

Applies a truncated SVD to the neighbour submatrix and returns the row of \(\widehat F_{\langle i \rangle} \widehat\Lambda_{\langle i \rangle}^\top\) belonging to unit – the estimate of \(\eta_l(\alpha_i)\) in paper equation (3.3).

Only one row of the reconstruction is ever read and the right singular vectors cancel from it, so the left factor of the K x K Gram matrix suffices: the row is (U[pos, :r] U[:, :r]') B. That agrees with the K x p SVD to about 5e-14 and is several times faster.

Parameters:

• block (np.ndarray) – PCA block, shape (n_units, n_pca_periods).

• neighbours (np.ndarray) – Boolean neighbour indicator for unit (a column of neighbour_index()).

• unit (int) – Row index of the unit being reconstructed.

• n_neighbours (int) – Requested K, passed to select_rank() for its log log K.

• n_components (int) – Maximum number of components extracted (the paper’s nlam).

Returns:

• tuple – The reconstructed row, the rank actually used, and the weight each neighbour receives.

• The weights are the unit-th column of the rank-r projector

• \(U_r U_r^\top\), so the reconstruction is literally a weighted

• combination of the neighbourhood’s rows. They neither sum to one nor stay

• non-negative – this is a projection, not a simplex fit. The unit’s own

• entry is the interesting one (the treated row is part of its own)

• neighbourhood, so that weight measures how much the reconstruction leans on

• the cells the estimator zeroed, which is the perturbation Theorem 6.1

• bounds.

mlsynth.utils.lpca_helpers.core.neighbour_index(distance: ndarray, n_neighbours: int) → ndarray

Boolean neighbourhood indicator, one column per unit.

The reference rule is a threshold, tmp.d <= nth(tmp.d, K), so exact ties are all kept and the realised neighbourhood can exceed K. Binary or otherwise discrete panels tie routinely, which is why the estimator reports the realised size next to the requested one.

Parameters:

• distance (np.ndarray) – Square distance matrix from pseudo_max_distance().

• n_neighbours (int) – Requested neighbourhood size K, counting the unit itself.

Returns:

np.ndarray – Boolean (n_units, n_units); column i flags unit i’s neighbours.

mlsynth.utils.lpca_helpers.core.pseudo_max_distance(A: ndarray, metric: str = 'pseudo_max') → ndarray

Pairwise distances between units on the matching block.

Three of the paper’s Section 3 choices are available. The default is the pseudo-max distance of Zhang, Levina & Zhu (2017),

\[\rho(x_i, x_j) = \max_{l \neq i, j} \bigl| \tfrac{1}{p} (x_i - x_j)^\top x_l \bigr|,\]

which averages over features and so separates units by their noise-free factor structure even when the errors are heteroskedastic; the Euclidean distance does not. "euclidean" is \(p^{-1/2} \lVert x_i - x_j \rVert\) and "average" is \(p^{-1} | \mathbf{1}^\top (x_i - x_j) |\).

Parameters:

• A (np.ndarray) – Matching block, shape (n_units, n_match_periods).

• metric ({“pseudo_max”, “euclidean”, “average”})

Returns:

np.ndarray – Symmetric (n_units, n_units) distances with a zero diagonal, so each unit is its own nearest neighbour – paper equation (3.1) includes i in its own neighbourhood.

Raises:

MlsynthConfigError – If metric is not one of the three supported names.

mlsynth.utils.lpca_helpers.core.select_rank(singular_values: ndarray, n_neighbours: int) → int

Number of local components to keep (R: sel.r).

Consecutive singular values whose ratio falls below log log K count as indistinguishable from noise, so the rank is the position of the first such pair, floored at one; if every ratio clears the threshold the rank is len(singular_values) - 1.

The threshold has a consequence the paper does not draw out. Ratios of a descending spectrum are at least 1, and log log K < 1 for K <= 15 (it crosses 1 at e^e = 15.15), so below sixteen neighbours the comparison never fires and the rank is pinned at len(singular_values) - 1 however the data behaves. Feng’s Kansas application sets K = 14. Paper Remark 4.3 leaves formal rank selection to future work, so the estimator surfaces the realised rank as a diagnostic.

Data preparation – the DataFrame touchpoint: pivots to the outcome matrix, masks the treated cells, and enforces the block split.

Panel ingestion for the LPCA estimator.

mlsynth.utils.lpca_helpers.setup.prepare_lpca_inputs(df: DataFrame, outcome: str, treat: str, unitid: str, time: str, *, match_periods: int | None = None, match_fraction: float = 0.5) → LPCAInputs

Pivot a long panel into LPCAInputs.

Ingestion goes through mlsynth.utils.datautils.dataprep(). The treated units’ post-treatment cells become the missing entries local PCA imputes.

Raises:

MlsynthDataError – If the panel is unbalanced, has fewer than three units, is treated in the first period, or is split so that the PCA block holds no pre-treatment period.

mlsynth.utils.lpca_helpers.setup.resolve_match_periods(n_periods: int, match_periods: int | None, match_fraction: float) → int

Number of leading periods reserved for matching.

mlsynth.utils.lpca_helpers.setup.resolve_neighbours(n_units: int, n_neighbors: int | None) → int

Neighbourhood size K, defaulting to the paper’s round(N^(2/3)).

Run loop: centring, matching, the local fits, and the re-centred counterfactual.

The LPCA causal pipeline: preprocess, match, fit locally, re-centre.

mlsynth.utils.lpca_helpers.pipeline.run_lpca(inputs: LPCAInputs, *, n_neighbours: int, max_components: int = 3, distance: str = 'pseudo_max', demean: bool = True) → LPCAResults

Run local PCA and assemble LPCAResults.

Parameters:

• inputs (LPCAInputs)

• n_neighbours (int) – Requested neighbourhood size K, resolved by the caller.

• max_components (int) – Maximum local components extracted; the realised rank comes from the singular-value ratio rule.

• distance (str) – Matching distance, one of METRICS.

• demean (bool) – Centre each period across units before matching. The period means are added back to the counterfactual, so the returned path is always on the scale of the outcome column.

Raises:

MlsynthDataError – If the requested neighbourhood exceeds the number of units, or exceeds what the local SVD can use.