Skip to content

mfe.crosssection

mfe.crosssection

mfe.crosssection — Cross-sectional econometrics.

ols OLS with White heteroskedastic SEs olsnw OLS with Newey-West HAC SEs fama_macbeth Two-pass FM regression with Shanken correction rolling_betas Rolling time-series betas for FM pass 1 pca Principal component analysis with factor interpretation

FMResult dataclass 

FMResult(lambda_mean: FloatArray, lambda_std: FloatArray, t_stats: FloatArray, p_values: FloatArray, t_stats_shanken: FloatArray, p_values_shanken: FloatArray, lambda_series: FloatArray, r_squared_mean: float, n_periods: int, n_assets: int, factor_names: list[str])

Fama-MacBeth estimation result.

PCAResult dataclass 

PCAResult(eigenvalues: FloatArray, eigenvectors: FloatArray, factors: FloatArray, loadings: FloatArray, explained_variance: FloatArray, cumulative_variance: FloatArray, n_components: int, n_obs: int, n_vars: int, mean: FloatArray)

Principal component analysis result.

reconstruct 

reconstruct(k_c: int | None = None) -> FloatArray

Reconstruct data from the first k_c PCA components.

Parameters:

Name Type Description Default
k_c number of components to use; if None uses all n_components
None

Returns:

Type Description
(T, K) reconstructed data matrix (in original scale, mean added back)
Source code in src/mfe/crosssection/pca.py 
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
def reconstruct(self, k_c: int | None = None) -> FloatArray:
    """
    Reconstruct data from the first k_c PCA components.

    Parameters
    ----------
    k_c : number of components to use; if None uses all n_components

    Returns
    -------
    (T, K) reconstructed data matrix (in original scale, mean added back)
    """
    k_c = k_c or self.n_components
    k_c = min(k_c, self.n_components)
    # factors = Xc @ eigvecs[:, :k_c]  (projection)
    # reconstruction = factors @ eigvecs[:, :k_c].T + mean
    # eigvecs[:, :k_c] from PCAResult.eigenvectors
    evecs = self.eigenvectors[:, :k_c]
    return self.factors[:, :k_c] @ evecs.T + self.mean[None, :]

olsnw 

olsnw(y: FloatArray, X: FloatArray, include_const: bool = True, nw_lags: int | None = None) -> OLSResult

OLS regression with Newey-West HAC standard errors.

Parameters:

Name Type Description Default
y FloatArray
required
X FloatArray
required
include_const prepend a constant (default True)
True
nw_lags int | None 
        Set to 0 for White-only (no serial correlation correction)
None

Returns:

Type Description
OLSResult with .vcv_robust = Newey-West VCV, .std_errors = NW standard errors.
Source code in src/mfe/crosssection/ols.py 
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
def olsnw(
    y: FloatArray,
    X: FloatArray,
    include_const: bool = True,
    nw_lags: int | None = None,
) -> OLSResult:
    """
    OLS regression with Newey-West HAC standard errors.

    Parameters
    ----------
    y             : (T,) dependent variable
    X             : (T, K) regressors — do NOT include a constant column
    include_const : prepend a constant (default True)
    nw_lags       : Newey-West bandwidth; if None uses floor(T^{1/3})
                    Set to 0 for White-only (no serial correlation correction)

    Returns
    -------
    OLSResult with .vcv_robust = Newey-West VCV, .std_errors = NW standard errors.
    """
    y = np.asarray(y, dtype=np.float64).ravel()
    X = np.asarray(X, dtype=np.float64)
    if X.ndim == 1:
        X = X[:, None]

    if include_const:
        X = np.column_stack([np.ones(len(y)), X])

    n, k = X.shape
    XTX_inv = np.linalg.inv(X.T @ X)
    beta = XTX_inv @ (X.T @ y)
    resid = y - X @ beta

    if nw_lags is None:
        nw_lags = int(np.floor(n ** (1 / 3)))

    scores = X * resid[:, None]   # (n, k)
    B_nw = newey_west(scores, bandwidth=nw_lags)
    vcv_nw = XTX_inv @ B_nw @ XTX_inv / n

    return _compute_ols_result(y, X, vcv_nw, include_const)

fama_macbeth 

fama_macbeth(returns: FloatArray, betas: FloatArray, include_intercept: bool = True, nw_lags: int = 0, shanken_correction: bool = True) -> FMResult

Two-pass Fama-MacBeth regression.

Parameters:

Name Type Description Default
returns FloatArray
required
betas FloatArray
required
include_intercept bool
True
nw_lags int
0
shanken_correction bool
True

Returns:

Type Description
FMResult
Source code in src/mfe/crosssection/fm.py 
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
def fama_macbeth(
    returns: FloatArray,             # (T, N) — T periods, N assets
    betas: FloatArray,               # (N, K) — pre-estimated factor loadings
    include_intercept: bool = True,
    nw_lags: int = 0,
    shanken_correction: bool = True,
) -> FMResult:
    """
    Two-pass Fama-MacBeth regression.

    Parameters
    ----------
    returns          : (T, N) panel of asset returns
    betas            : (N, K) pre-estimated betas (from pass 1 or exogenous)
    include_intercept: add a constant to the cross-sectional regression
    nw_lags          : Newey-West lags for FM standard errors; 0 = no HAC
    shanken_correction: apply Shanken (1992) EIV correction to t-stats

    Returns
    -------
    FMResult
    """
    R = np.asarray(returns, dtype=np.float64)   # (T, N)
    B = np.asarray(betas, dtype=np.float64)     # (N, K)
    T, N = R.shape
    _, K = B.shape

    if include_intercept:
        X = np.column_stack([np.ones(N), B])  # (N, K+1)
        k_total = K + 1
    else:
        X = B
        k_total = K

    XTX_inv = np.linalg.pinv(X.T @ X)  # (k_total, k_total)

    # Pass 2: cross-sectional OLS at each t
    lambda_t = np.empty((T, k_total), dtype=np.float64)
    r2_t = np.empty(T, dtype=np.float64)

    for t in range(T):
        r_t = R[t]  # (N,)
        lam = XTX_inv @ (X.T @ r_t)
        lambda_t[t] = lam
        fitted = X @ lam
        ss_res = np.sum((r_t - fitted) ** 2)
        ss_tot = np.sum((r_t - np.mean(r_t)) ** 2)
        r2_t[t] = 1 - ss_res / ss_tot if ss_tot > 0 else 0.0

    # FM standard errors
    lam_mean = np.mean(lambda_t, axis=0)
    lam_std = np.std(lambda_t, axis=0, ddof=1)

    if nw_lags > 0:
        # Newey-West on the lambda_t series
        centered = lambda_t - lam_mean[None, :]
        B_nw = newey_west(centered, bandwidth=nw_lags)
        fm_var = np.diag(B_nw) / T
    else:
        fm_var = lam_std ** 2 / T

    fm_se = np.sqrt(np.maximum(fm_var, 0.0))
    t_stats = np.where(fm_se > 0, lam_mean / fm_se, np.nan)
    p_values = 2 * (1 - stats.t.cdf(np.abs(t_stats), df=T - 1))

    # Shanken (1992) correction
    # c = 1 + lambda_f' Sigma_f^{-1} lambda_f  (where lambda_f are the premia on factors)
    if shanken_correction and K > 0:
        lam_f = lam_mean[1:] if include_intercept else lam_mean  # factor premia only
        try:
            Sigma_f = np.cov(lambda_t[:, 1:].T) if include_intercept else np.cov(lambda_t.T)
            c = 1.0 + float(lam_f @ np.linalg.solve(Sigma_f, lam_f))
        except np.linalg.LinAlgError:
            c = 1.0
        shanken_se = fm_se * np.sqrt(c)
        t_stats_s = np.where(shanken_se > 0, lam_mean / shanken_se, np.nan)
        p_values_s = 2 * (1 - stats.t.cdf(np.abs(t_stats_s), df=T - 1))
    else:
        t_stats_s = t_stats.copy()
        p_values_s = p_values.copy()

    factor_names = (["alpha"] if include_intercept else []) + [f"factor_{k+1}" for k in range(K)]

    return FMResult(
        lambda_mean=lam_mean,
        lambda_std=lam_std,
        t_stats=t_stats,
        p_values=p_values,
        t_stats_shanken=t_stats_s,
        p_values_shanken=p_values_s,
        lambda_series=lambda_t,
        r_squared_mean=float(np.mean(r2_t)),
        n_periods=T,
        n_assets=N,
        factor_names=factor_names,
    )

rolling_betas 

rolling_betas(returns: FloatArray, factors: FloatArray, window: int = 60) -> FloatArray

Rolling time-series OLS betas: for each asset n, regress returns on factors using a trailing window.

Returns (N, K) array of end-of-sample betas. Useful for the first pass of Fama-MacBeth.

Source code in src/mfe/crosssection/fm.py 
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
def rolling_betas(
    returns: FloatArray,    # (T, N)
    factors: FloatArray,    # (T, K)
    window: int = 60,
) -> FloatArray:
    """
    Rolling time-series OLS betas: for each asset n, regress returns on factors
    using a trailing window.

    Returns (N, K) array of end-of-sample betas.
    Useful for the first pass of Fama-MacBeth.
    """
    R = np.asarray(returns, dtype=np.float64)
    F = np.asarray(factors, dtype=np.float64)
    T, N = R.shape
    _, K = F.shape

    # Use full sample if T < window
    w = min(window, T)
    F_w = F[-w:]
    X = np.column_stack([np.ones(w), F_w])  # (w, K+1)

    betas = np.empty((N, K), dtype=np.float64)
    XTX_inv = np.linalg.pinv(X.T @ X)

    for n in range(N):
        r_w = R[-w:, n]
        coef = XTX_inv @ (X.T @ r_w)
        betas[n] = coef[1:]  # drop intercept

    return betas

fm

Fama-MacBeth two-pass cross-sectional regression.

Fama, E.F. & MacBeth, J.D. (1973): "Risk, Return, and Equilibrium: Empirical Tests", Journal of Political Economy.

Two passes: Pass 1: For each time period t, regress cross-sectional returns on factor loadings (betas) to get factor risk premia lambda_t. Pass 2: Average lambda_t across time and compute t-statistics with Shanken (1992) correction for errors-in-variables.

Also implements rolling-window beta estimation (first step of pass 1).

FMResult dataclass 

FMResult(lambda_mean: FloatArray, lambda_std: FloatArray, t_stats: FloatArray, p_values: FloatArray, t_stats_shanken: FloatArray, p_values_shanken: FloatArray, lambda_series: FloatArray, r_squared_mean: float, n_periods: int, n_assets: int, factor_names: list[str])

Fama-MacBeth estimation result.

fama_macbeth 

fama_macbeth(returns: FloatArray, betas: FloatArray, include_intercept: bool = True, nw_lags: int = 0, shanken_correction: bool = True) -> FMResult

Two-pass Fama-MacBeth regression.

Parameters:

Name Type Description Default
returns FloatArray
required
betas FloatArray
required
include_intercept bool
True
nw_lags int
0
shanken_correction bool
True

Returns:

Type Description
FMResult
Source code in src/mfe/crosssection/fm.py 
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
def fama_macbeth(
    returns: FloatArray,             # (T, N) — T periods, N assets
    betas: FloatArray,               # (N, K) — pre-estimated factor loadings
    include_intercept: bool = True,
    nw_lags: int = 0,
    shanken_correction: bool = True,
) -> FMResult:
    """
    Two-pass Fama-MacBeth regression.

    Parameters
    ----------
    returns          : (T, N) panel of asset returns
    betas            : (N, K) pre-estimated betas (from pass 1 or exogenous)
    include_intercept: add a constant to the cross-sectional regression
    nw_lags          : Newey-West lags for FM standard errors; 0 = no HAC
    shanken_correction: apply Shanken (1992) EIV correction to t-stats

    Returns
    -------
    FMResult
    """
    R = np.asarray(returns, dtype=np.float64)   # (T, N)
    B = np.asarray(betas, dtype=np.float64)     # (N, K)
    T, N = R.shape
    _, K = B.shape

    if include_intercept:
        X = np.column_stack([np.ones(N), B])  # (N, K+1)
        k_total = K + 1
    else:
        X = B
        k_total = K

    XTX_inv = np.linalg.pinv(X.T @ X)  # (k_total, k_total)

    # Pass 2: cross-sectional OLS at each t
    lambda_t = np.empty((T, k_total), dtype=np.float64)
    r2_t = np.empty(T, dtype=np.float64)

    for t in range(T):
        r_t = R[t]  # (N,)
        lam = XTX_inv @ (X.T @ r_t)
        lambda_t[t] = lam
        fitted = X @ lam
        ss_res = np.sum((r_t - fitted) ** 2)
        ss_tot = np.sum((r_t - np.mean(r_t)) ** 2)
        r2_t[t] = 1 - ss_res / ss_tot if ss_tot > 0 else 0.0

    # FM standard errors
    lam_mean = np.mean(lambda_t, axis=0)
    lam_std = np.std(lambda_t, axis=0, ddof=1)

    if nw_lags > 0:
        # Newey-West on the lambda_t series
        centered = lambda_t - lam_mean[None, :]
        B_nw = newey_west(centered, bandwidth=nw_lags)
        fm_var = np.diag(B_nw) / T
    else:
        fm_var = lam_std ** 2 / T

    fm_se = np.sqrt(np.maximum(fm_var, 0.0))
    t_stats = np.where(fm_se > 0, lam_mean / fm_se, np.nan)
    p_values = 2 * (1 - stats.t.cdf(np.abs(t_stats), df=T - 1))

    # Shanken (1992) correction
    # c = 1 + lambda_f' Sigma_f^{-1} lambda_f  (where lambda_f are the premia on factors)
    if shanken_correction and K > 0:
        lam_f = lam_mean[1:] if include_intercept else lam_mean  # factor premia only
        try:
            Sigma_f = np.cov(lambda_t[:, 1:].T) if include_intercept else np.cov(lambda_t.T)
            c = 1.0 + float(lam_f @ np.linalg.solve(Sigma_f, lam_f))
        except np.linalg.LinAlgError:
            c = 1.0
        shanken_se = fm_se * np.sqrt(c)
        t_stats_s = np.where(shanken_se > 0, lam_mean / shanken_se, np.nan)
        p_values_s = 2 * (1 - stats.t.cdf(np.abs(t_stats_s), df=T - 1))
    else:
        t_stats_s = t_stats.copy()
        p_values_s = p_values.copy()

    factor_names = (["alpha"] if include_intercept else []) + [f"factor_{k+1}" for k in range(K)]

    return FMResult(
        lambda_mean=lam_mean,
        lambda_std=lam_std,
        t_stats=t_stats,
        p_values=p_values,
        t_stats_shanken=t_stats_s,
        p_values_shanken=p_values_s,
        lambda_series=lambda_t,
        r_squared_mean=float(np.mean(r2_t)),
        n_periods=T,
        n_assets=N,
        factor_names=factor_names,
    )

rolling_betas 

rolling_betas(returns: FloatArray, factors: FloatArray, window: int = 60) -> FloatArray

Rolling time-series OLS betas: for each asset n, regress returns on factors using a trailing window.

Returns (N, K) array of end-of-sample betas. Useful for the first pass of Fama-MacBeth.

Source code in src/mfe/crosssection/fm.py 
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
def rolling_betas(
    returns: FloatArray,    # (T, N)
    factors: FloatArray,    # (T, K)
    window: int = 60,
) -> FloatArray:
    """
    Rolling time-series OLS betas: for each asset n, regress returns on factors
    using a trailing window.

    Returns (N, K) array of end-of-sample betas.
    Useful for the first pass of Fama-MacBeth.
    """
    R = np.asarray(returns, dtype=np.float64)
    F = np.asarray(factors, dtype=np.float64)
    T, N = R.shape
    _, K = F.shape

    # Use full sample if T < window
    w = min(window, T)
    F_w = F[-w:]
    X = np.column_stack([np.ones(w), F_w])  # (w, K+1)

    betas = np.empty((N, K), dtype=np.float64)
    XTX_inv = np.linalg.pinv(X.T @ X)

    for n in range(N):
        r_w = R[-w:, n]
        coef = XTX_inv @ (X.T @ r_w)
        betas[n] = coef[1:]  # drop intercept

    return betas

ols

OLS and OLSNW regression — MFE toolbox ols.m / olsnw.m equivalents.

Design rationale: statsmodels OLS exists but requires a DataFrame/array and returns a result object with a non-trivial API. These functions are thin, fast wrappers that match the MFE toolbox calling convention exactly and fit naturally into pipelines that already use mfe.utils.vcv.

ols(Y, X) — OLS with White heteroskedastic SEs olsnw(Y, X) — OLS with Newey-West HAC SEs

ols 

ols(y: FloatArray, X: FloatArray, include_const: bool = True) -> OLSResult

OLS regression with White heteroskedasticity-robust standard errors.

Parameters:

Name Type Description Default
y FloatArray
required
X FloatArray
required
include_const prepend a column of ones (default True)
True

Returns:

Type Description
OLSResult

.params — coefficient vector [const?, beta_1..beta_K] .t_stats — White-robust t-statistics .std_errors — White-robust standard errors .vcv — classical (homoskedastic) VCV .vcv_robust — White heteroskedasticity-robust VCV
Source code in src/mfe/crosssection/ols.py 
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
def ols(
    y: FloatArray,
    X: FloatArray,
    include_const: bool = True,
) -> OLSResult:
    """
    OLS regression with White heteroskedasticity-robust standard errors.

    Parameters
    ----------
    y             : (T,) dependent variable
    X             : (T, K) regressors — do NOT include a constant column
    include_const : prepend a column of ones (default True)

    Returns
    -------
    OLSResult
        .params       — coefficient vector [const?, beta_1..beta_K]
        .t_stats      — White-robust t-statistics
        .std_errors   — White-robust standard errors
        .vcv          — classical (homoskedastic) VCV
        .vcv_robust   — White heteroskedasticity-robust VCV
    """
    y = np.asarray(y, dtype=np.float64).ravel()
    X = np.asarray(X, dtype=np.float64)
    if X.ndim == 1:
        X = X[:, None]

    if include_const:
        X = np.column_stack([np.ones(len(y)), X])

    n, k = X.shape
    XTX_inv = np.linalg.inv(X.T @ X)
    beta = XTX_inv @ (X.T @ y)
    resid = y - X @ beta

    # White VCV: (X'X)^{-1} * (sum e_t^2 x_t x_t') * (X'X)^{-1}
    S = (X * resid[:, None]).T @ (X * resid[:, None]) / n
    vcv_white = XTX_inv @ S @ XTX_inv / n

    return _compute_ols_result(y, X, vcv_white, include_const)

olsnw 

olsnw(y: FloatArray, X: FloatArray, include_const: bool = True, nw_lags: int | None = None) -> OLSResult

OLS regression with Newey-West HAC standard errors.

Parameters:

Name Type Description Default
y FloatArray
required
X FloatArray
required
include_const prepend a constant (default True)
True
nw_lags int | None 
        Set to 0 for White-only (no serial correlation correction)
None

Returns:

Type Description
OLSResult with .vcv_robust = Newey-West VCV, .std_errors = NW standard errors.
Source code in src/mfe/crosssection/ols.py 
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
def olsnw(
    y: FloatArray,
    X: FloatArray,
    include_const: bool = True,
    nw_lags: int | None = None,
) -> OLSResult:
    """
    OLS regression with Newey-West HAC standard errors.

    Parameters
    ----------
    y             : (T,) dependent variable
    X             : (T, K) regressors — do NOT include a constant column
    include_const : prepend a constant (default True)
    nw_lags       : Newey-West bandwidth; if None uses floor(T^{1/3})
                    Set to 0 for White-only (no serial correlation correction)

    Returns
    -------
    OLSResult with .vcv_robust = Newey-West VCV, .std_errors = NW standard errors.
    """
    y = np.asarray(y, dtype=np.float64).ravel()
    X = np.asarray(X, dtype=np.float64)
    if X.ndim == 1:
        X = X[:, None]

    if include_const:
        X = np.column_stack([np.ones(len(y)), X])

    n, k = X.shape
    XTX_inv = np.linalg.inv(X.T @ X)
    beta = XTX_inv @ (X.T @ y)
    resid = y - X @ beta

    if nw_lags is None:
        nw_lags = int(np.floor(n ** (1 / 3)))

    scores = X * resid[:, None]   # (n, k)
    B_nw = newey_west(scores, bandwidth=nw_lags)
    vcv_nw = XTX_inv @ B_nw @ XTX_inv / n

    return _compute_ols_result(y, X, vcv_nw, include_const)

pca

Principal Component Analysis for financial returns.

Standalone PCA module with the interface matching the MFE MATLAB pca.m: - eigendecomposition of the sample covariance - proportion of variance explained per component - factor scores (principal components) - loadings matrix - reconstruction of the original data from K_c components

This is a clean public API wrapping the internals already used in mfe.multivariate.gogarch. The MFE MATLAB pca.m is documented but rarely exposed directly — we make it first-class here.

Distinct from sklearn.decomposition.PCA in that: - we expose the covariance structure (not correlation by default) - we match financial conventions: K_c components from covariance, not correlation, and we return eigenvalues in variance units (not explained variance ratio) alongside the standard stats

PCAResult dataclass 

PCAResult(eigenvalues: FloatArray, eigenvectors: FloatArray, factors: FloatArray, loadings: FloatArray, explained_variance: FloatArray, cumulative_variance: FloatArray, n_components: int, n_obs: int, n_vars: int, mean: FloatArray)

Principal component analysis result.

reconstruct 
reconstruct(k_c: int | None = None) -> FloatArray

Reconstruct data from the first k_c PCA components.

Parameters:

Name Type Description Default
k_c number of components to use; if None uses all n_components
None

Returns:

Type Description
(T, K) reconstructed data matrix (in original scale, mean added back)
Source code in src/mfe/crosssection/pca.py 
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
def reconstruct(self, k_c: int | None = None) -> FloatArray:
    """
    Reconstruct data from the first k_c PCA components.

    Parameters
    ----------
    k_c : number of components to use; if None uses all n_components

    Returns
    -------
    (T, K) reconstructed data matrix (in original scale, mean added back)
    """
    k_c = k_c or self.n_components
    k_c = min(k_c, self.n_components)
    # factors = Xc @ eigvecs[:, :k_c]  (projection)
    # reconstruction = factors @ eigvecs[:, :k_c].T + mean
    # eigvecs[:, :k_c] from PCAResult.eigenvectors
    evecs = self.eigenvectors[:, :k_c]
    return self.factors[:, :k_c] @ evecs.T + self.mean[None, :]

pca 

pca(data: FloatArray, n_components: int | None = None, demean: bool = True, standardize: bool = False) -> PCAResult

Principal component analysis of a (T, K) data matrix.

Parameters:

Name Type Description Default
data FloatArray
required
n_components number of components to retain; if None, keeps all K
None
demean bool
True
standardize bool 
       if False (default), operates on covariance matrix
False

Returns:

Type Description
PCAResult

.eigenvalues — (K,) eigenvalues of sample covariance/correlation matrix .eigenvectors — (K, K) eigenvector matrix (columns sorted descending) .factors — (T, K_c) factor scores = demeaned_data @ eigenvectors[:, :K_c] .loadings — (K, K_c) factor loadings scaled by sqrt(eigenvalue) .explained_variance — proportion of total variance per component
Source code in src/mfe/crosssection/pca.py 
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
def pca(
    data: FloatArray,
    n_components: int | None = None,
    demean: bool = True,
    standardize: bool = False,
) -> PCAResult:
    """
    Principal component analysis of a (T, K) data matrix.

    Parameters
    ----------
    data         : (T, K) matrix (e.g. asset returns)
    n_components : number of components to retain; if None, keeps all K
    demean       : subtract column means before decomposition (default True)
    standardize  : divide by column std after demeaning (correlation PCA);
                   if False (default), operates on covariance matrix

    Returns
    -------
    PCAResult
        .eigenvalues  — (K,) eigenvalues of sample covariance/correlation matrix
        .eigenvectors — (K, K) eigenvector matrix (columns sorted descending)
        .factors      — (T, K_c) factor scores = demeaned_data @ eigenvectors[:, :K_c]
        .loadings     — (K, K_c) factor loadings scaled by sqrt(eigenvalue)
        .explained_variance — proportion of total variance per component
    """
    X = np.asarray(data, dtype=np.float64)
    T, K = X.shape
    K_c = K if n_components is None else min(n_components, K)

    # Demean
    mu = X.mean(axis=0) if demean else np.zeros(K, dtype=np.float64)
    Xc = X - mu[None, :]

    # Optionally standardize
    if standardize:
        std = Xc.std(axis=0, ddof=1)
        std = np.where(std > 0, std, 1.0)
        Xc = Xc / std[None, :]
    else:
        std = np.ones(K, dtype=np.float64)

    # Sample covariance (no Bessel correction to match MATLAB mfe convention)
    S = Xc.T @ Xc / T

    # Eigendecomposition (eigh is stable for symmetric matrices)
    eigvals, eigvecs = np.linalg.eigh(S)
    # Sort descending
    idx = np.argsort(eigvals)[::-1]
    eigvals = eigvals[idx]
    eigvecs = eigvecs[:, idx]

    # Clip near-zero eigenvalues
    eigvals_safe = np.maximum(eigvals, 0.0)

    # Factor scores: projection of centered data onto eigenvectors
    factors = Xc @ eigvecs[:, :K_c]  # (T, K_c) — unit-variance if we scale

    # Loadings: eigenvectors scaled by sqrt(eigenvalue) → covariance structure
    loadings = eigvecs[:, :K_c] * np.sqrt(eigvals_safe[:K_c])[None, :]  # (K, K_c)

    # Explained variance
    total_var = float(np.sum(eigvals_safe))
    exp_var = eigvals_safe[:K_c] / max(total_var, 1e-30)
    cum_var = np.cumsum(exp_var)

    return PCAResult(
        eigenvalues=eigvals,
        eigenvectors=eigvecs,
        factors=factors,
        loadings=loadings,
        explained_variance=exp_var,
        cumulative_variance=cum_var,
        n_components=K_c,
        n_obs=T,
        n_vars=K,
        mean=mu,
    )