Линейная регрессия в факторных моделях

Привет, Хабр!

Когда мы говорим «факторная модель», многие вспоминают Python‑ноутбуки. Но если отмотать плёнку, бóльшая часть индустриальных движков для риска и ценообразования десятилетиями писалась на C++ поверх BLAS/LAPACK. Там же удобно делать устойчивые разложения: QR с переупорядочиванием столбцов, SVD, регуляризацию. Библиотеки вроде Eigen дали нормальный интерфейс к этим штукам, и регрессия перестала быть болью «Ax = b» руками. QR с перестановками колонок вообще стандарт для переобусловленных задач.

Сама идея факторной модели пришла не из тетрадки с pandas, а из арбитражной теории ценообразования Россa и последующей эмпирики Fama‑French. В терминах работы это выглядит как линейная регрессия доходностей на набор общих факторов. Дальше есть два пути проверки: тайм‑серия для бета‑нагрузок и кросс‑секция для премий за риск. Это конвейер, а не разовая регрессия.

Что именно считаем в факторной модели

Разделяем процесс на четыре шага.

Выбор факторов и тест‑активов. В примерах берут Fama‑French: MKT‑RF, SMB, HML, а ещё часто добавляют momentum из Carhart.
Тайм‑серийная регрессия для бета‑нагрузок βᵢ. Классика: для каждого актива i оцениваем
rᵢₜ − r_fₜ = αᵢ + βᵢ⊤ fₜ + εᵢₜ,
где fₜ это вектор факторных доходностей.
Кросс‑секционная регрессия премий λₜ (Fama‑MacBeth): в каждый момент времени t решаем
r̄ₜ = 1·γ₀ₜ + B·λₜ + uₜ,
где r̄ₜ это вектор доходностей активов в момент t, B матрица ранее оценённых β. Потом усредняем λₜ и считаем t‑статистики.
Построение risk‑модели: ковариация активов Σ ≈ B Σ_f B⊤ + Δ, где Σ_f ковариация факторов, Δ диагональная специфическая дисперсия.

Дальше посмотрим Python‑пайплайн, затем кусок C++ на Eigen для ядра OLS/риджа, потом PCA‑вариант статистической факторной модели и финишируем частью про диагностику, устойчивые ковариации и тесты.

Подготовка данных

Часть, от которой зависят метрики. Несколько правил гигиены:

• Синхронизация по датам. Факторы и активы должны быть выровнены по частоте и календарю, без look‑ahead.

• Доходности: решите, какие используете — простые или логарифмические, будьте последовательны.

• Отсев активов с короткой историей, но не так, чтобы словить survivorship bias.

• Контроль пропусков и знакомых праздников/пересчётов.

Создадим минимальный каркас на pandas/statsmodels. Намеренно не тянем данные из сети — кладём CSV рядом: factors.csv (MKT_RF, SMB, HML, UMD, RF) и returns.csv (типа 100–300 акций с колонками тикеров). Премии считаются на избыточных доходностях.

# pyproject.toml: statsmodels==0.14.*, pandas>=2.0, numpy>=1.26, scikit-learn>=1.4 import pandas as pd import numpy as np  factors = (     pd.read_csv("factors.csv", parse_dates=["date"])       .set_index("date")       .rename(columns=str.upper) )  # обязательные столбцы: MKT_RF, SMB, HML, RF; опционально UMD (momentum) assert {"MKT_RF", "SMB", "HML", "RF"}.issubset(factors.columns)  rets = (     pd.read_csv("returns.csv", parse_dates=["date"])       .set_index("date") )  # оставим только общие даты и отсортируем df = rets.join(factors, how="inner").sort_index() # избыточные доходности активов excess = df[rets.columns].sub(df["RF"], axis=0).dropna(how="all") F = df[["MKT_RF", "SMB", "HML"]].loc[excess.index].copy() if "UMD" in df:     F["UMD"] = df["UMD"].loc[excess.index]  # фильтр активов с достаточной историей, например >= 60 наблюдений enough_hist = excess.count() >= 60 excess = excess.loc[:, enough_hist]

Тайм-серийные β с устойчивыми стандартными ошибками

OLS как есть даст смещённые стандартные ошибки, потому что есть гетероскедастичность и автокорреляция остатков. В тайм‑сериях берут HAC оценки ковариации по Newey‑West. В statsmodels это делается одной строчкой.

import statsmodels.api as sm  fcols = F.columns X = sm.add_constant(F.values)  # колонка 1 для альфы betas = {} alphas = {} nw_tstats = {}  for ticker in excess.columns:     y = excess[ticker].values     model = sm.OLS(y, X, missing="drop")     res = model.fit()     # HAC (Newey-West) с лагом, например floor(4*(T/100)^(2/9)) или по правилу Шварца     nlags = int(np.floor(4 * (len(y)/100)**(2/9))) or 1     res_hac = res.get_robustcov_results(cov_type="HAC", maxlags=nlags)     alphas[ticker] = res_hac.params[0]     betas[ticker]  = res_hac.params[1:]     # сохраним t-статистики по альфе для диагностики     nw_tstats[ticker] = res_hac.tvalues[0]  B = pd.DataFrame(betas, index=fcols).T  # shape: [assets x factors] alpha_series = pd.Series(alphas)

Почему HAC, а не White HC? Потому что в финансовых рядах автокорреляция и условная гетероскедастичность — это норма, и Newey‑West под это настроен. HAC корректирует ковариацию оценок, а не сами коэффициенты.

Кросс-секционные премии и корректные ошибки

Теперь оценим цепочку Fama‑MacBeth: на каждом t решаем cross‑sectional OLS rₜ = γ₀ₜ + B·λₜ + uₜ, затем усредняем λₜ. Для стандартных ошибок есть два классических подхода: «наивное» ст. отклонение по времени и поправка Shanken на ошибочность»β, которые мы внесли из первой стадии.

# временной цикл кросс-секций Rb = excess.loc[:, B.index]  # совпадаем по активам Bt = B.values  # фиксируем беты по времени (упрощение, можно рефитить в окне)  lam_t = [] for t, row in Rb.iterrows():     y = row.values  # доходности активов в момент t     X = np.c_[np.ones(len(Bt)), Bt]  # константа + беты     # решаем малую регрессию устойчивым способом     coef, *_ = np.linalg.lstsq(X, y, rcond=None)     lam_t.append(coef[1:])  # без константы  lam = np.mean(np.vstack(lam_t), axis=0) lam_series = pd.Series(lam, index=fcols, name="lambda_mean")

Наивные t‑статы по Fama‑MacBeth — это среднее λ делённое на стандартное отклонение λₜ по времени с поправкой на число наблюдений. Если хотите аккуратно, используйте HAC по времени для λₜ или сделайте Shanken‑коррекцию (ошибка‑в‑переменных). В statsmodels можно обернуть результаты cov_hac и получить аккуратные p‑values.

lam_T = np.vstack(lam_t)  # [T x K] T = lam_T.shape[0] lam_mean = lam_T.mean(axis=0) lam_std  = lam_T.std(axis=0, ddof=1) t_naive  = lam_mean / (lam_std / np.sqrt(T)) fm_summary = pd.DataFrame({"lambda": lam_mean, "t_naive": t_naive}, index=fcols)

На этом этапе полезно прогнать GRS‑тест на нулевые альфы из первой стадии. Если модель хорошая, α близки к нулю статистически и экономически.

Из регрессии в risk-модель

Практическая цель в портфеле — предсказать риск. Берём ковариацию факторов Σ_f по тайм‑серии, умножаем на B, добавляем специфический шум Δ по остаткам.

# оценка ковариации факторов и специфических дисперсий F_mat = F.values Sigma_f = np.cov(F_mat, rowvar=False, ddof=1)  resid = {} for i, ticker in enumerate(B.index):     # восстановим остатки, чтобы получить специфическую дисперсию     y = excess[ticker].values     y_hat = sm.add_constant(F_mat) @ np.r_[alpha_series[ticker], B.loc[ticker].values]     resid[ticker] = y - y_hat  spec_var = pd.Series({k: np.var(v, ddof=1) for k, v in resid.items()}) Delta = np.diag(spec_var.reindex(B.index).values)  Sigma_assets = B.values @ Sigma_f @ B.values.T + Delta  # ковариация активов

Для больших размерностей голая np.cov шумная. В реальном мире применяют shrinkage к Σ_f и к Δ. Ledoit‑Wolf — разумная дефолтная усадка к диагонали.

C++-ядро: устойчивый OLS и Ridge на Eigen

Сейчас короткий фрагмент платформенного ядра. Идея простая: используем ColPivHouseholderQR для устойчивого решения, а при подозрении на мультиколлинеарность — Ridge. Это то, что легко переносится в высоконагруженный бэкенд.

// CMake: find_package(Eigen3 REQUIRED) #include <Eigen/Dense> #include <vector>  struct OLSResult {     Eigen::VectorXd coef;   // [k+1]: intercept + k betas     double sigma2;          // residual variance     bool ok; };  OLSResult ols_qr(const Eigen::MatrixXd& X_in, const Eigen::VectorXd& y) {     OLSResult out{};     if (X_in.rows() != y.size() || X_in.rows() < X_in.cols()) return out;      Eigen::MatrixXd X = X_in; // копия, если нужно безопасно мутировать     Eigen::ColPivHouseholderQR<Eigen::MatrixXd> qr(X);     if (qr.rank() < X.cols()) {         // деградация: вернёмся позже к ridge     }     out.coef = qr.solve(y);     Eigen::VectorXd r = y - X * out.coef;     out.sigma2 = (r.array().square().sum()) / std::max<int>(1, X.rows() - X.cols());     out.ok = true;     return out; }

// Ridge: (X^T X + alpha I)^{-1} X^T y Eigen::VectorXd ridge(const Eigen::MatrixXd& X, const Eigen::VectorXd& y, double alpha) {     const int p = X.cols();     Eigen::MatrixXd XtX = X.transpose() * X;     XtX += alpha * Eigen::MatrixXd::Identity(p, p);     return XtX.ldlt().solve(X.transpose() * y); }

// Конструктор дизайн-матрицы: добавляем константу Eigen::MatrixXd with_intercept(const Eigen::MatrixXd& F) {     Eigen::MatrixXd X(F.rows(), F.cols() + 1);     X.col(0).setOnes();     X.rightCols(F.cols()) = F;     return X; }

QR с перестановками столбцов — компромисс между устойчивостью и скоростью. Для задач регрессии в риск‑движке этого достаточно, а где плохо — включаем ridge с небольшим альфа.

PCA-факторная модель как статистическая альтернатива

Бывает, что не нужны интерпретируемые SMB/HML, а нужна компактная статистическая модель риска: PCA на матрице доходностей, затем несколько главных компонент как факторы. Уместно для доходностей доходных кривых или при тысячах активов, где K << N.

from sklearn.decomposition import PCA  X = excess.fillna(0.0).values  # [T x N] pca = PCA(n_components=10, svd_solver="full", whiten=False, random_state=0) scores = pca.fit_transform(X)         # факторные доходности [T x K] loadings = pca.components_.T          # нагрузки [N x K] spec_var_pca = np.var(X - scores @ loadings.T, axis=0, ddof=1)  Sigma_f_pca = np.cov(scores, rowvar=False, ddof=1) Delta_pca = np.diag(spec_var_pca) Sigma_assets_pca = loadings @ Sigma_f_pca @ loadings.T + Delta_pca

Сколько компонент выбирать? Два критерия: доля объяснённой дисперсии, шагом добавляете компоненты до порога, например 70–90 процентов; либо во валидации по риску минимизируете ошибку предсказания волатильности портфелей. В scikit‑learn explained_variance_ratio_ поможет.

Для стабилизации оценок ковариации факторов и специфики уместно применить Ledoit‑Wolf. В sklearn есть готовый sklearn.covariance.LedoitWolf.

from sklearn.covariance import LedoitWolf  Sigma_f_pca = LedoitWolf().fit(scores).covariance_ # специфическую тоже можно усадить к диагонали средней дисперсии delta_target = np.full_like(spec_var_pca, spec_var_pca.mean()) shrink = 0.2 spec_var_shrunk = (1 - shrink) * spec_var_pca + shrink * delta_target Delta_pca = np.diag(spec_var_shrunk)

Когда регрессия врёт

Мультиколлинеарность факторов. Симптом: нестабильные β и взрывающиеся стандартные ошибки. Лечится регуляризацией (ridge), исключением лишних факторов, PCA‑предобработкой.
t‑статистика «2.0 достаточно?» — уже давно нет. Литература про факторный зоопарк убеждает смотреть на множественную проверку гипотез и поднимать пороги значимости, часто до |t| ≥ 3. Иначе ловите публикационный шум.
Стандартные ошибки. Для тайм‑серий — HAC Newey‑West. В Fama‑MacBeth используйте либо «по времени» с HAC, либо применяйте шанкeновскую поправку на ошибочность β.
GRS‑тест на нулевые альфы портфелей‑тестов. Если отклоняет, модель не объясняет кросс‑секцию.
Данные факторов. Поддерживайте воспроизводимость: фиксируйте конструкторы SMB/HML/UMD. Официальное описание методики Fama/French — не короткая заметка, там нюансы по портфелям 2×3, сортировкам и календарю.

Небольшой пайплайн оценки и применения

Соберём всё в одну функцию, которая: (а) оценивает β по HAC‑OLS, (б) делает Fama‑MacBeth, (в) строит риск‑модель и (г) считает риск заданного портфеля w.

from dataclasses import dataclass  @dataclass class FactorModel:     betas: pd.DataFrame     alphas: pd.Series     lambdas: pd.Series     Sigma_assets: np.ndarray     Sigma_f: np.ndarray     spec_var: pd.Series  def fit_factor_model(excess_ret: pd.DataFrame, factors_df: pd.DataFrame) -> FactorModel:     Fm = factors_df.copy()     X = sm.add_constant(Fm.values)     betas, alphas = {}, {}     for ticker in excess_ret.columns:         y = excess_ret[ticker].values         res = sm.OLS(y, X, missing="drop").fit()         nlags = int(np.floor(4 * (len(y)/100)**(2/9))) or 1         res_hac = res.get_robustcov_results(cov_type="HAC", maxlags=nlags)         alphas[ticker] = res_hac.params[0]         betas[ticker]  = res_hac.params[1:]     B = pd.DataFrame(betas, index=Fm.columns).T     # Fama-MacBeth     lam_t = []     for t, row in excess_ret.iterrows():         y = row[B.index].values         Xcs = np.c_[np.ones(len(B)), B.values]         coef, *_ = np.linalg.lstsq(Xcs, y, rcond=None)         lam_t.append(coef[1:])     lam = np.mean(np.vstack(lam_t), axis=0)     # ковариации     Sigma_f = np.cov(Fm.values, rowvar=False, ddof=1)     resid = {}     for ticker in B.index:         y = excess_ret[ticker].values         y_hat = sm.add_constant(Fm.values) @ np.r_[alphas[ticker], B.loc[ticker].values]         resid[ticker] = y - y_hat     spec_var = pd.Series({k: np.var(v, ddof=1) for k, v in resid.items()})     Delta = np.diag(spec_var.reindex(B.index).values)     Sigma_assets = B.values @ Sigma_f @ B.values.T + Delta     return FactorModel(         betas=B, alphas=pd.Series(alphas), lambdas=pd.Series(lam, index=Fm.columns),         Sigma_assets=Sigma_assets, Sigma_f=Sigma_f, spec_var=spec_var     )  fm = fit_factor_model(excess, F)

Применение к портфелю. Здесь важно помнить о нормировке весов и проверке размерностей.

def portfolio_var(w: np.ndarray, Sigma: np.ndarray) -> float:     w = np.asarray(w).reshape(-1)     assert Sigma.shape[0] == Sigma.shape[1] == w.size     return float(w @ Sigma @ w)  weights = np.ones(len(fm.betas)) / len(fm.betas) sigma_p = np.sqrt(portfolio_var(weights, fm.Sigma_assets))  # месячная волатильность

Частые вопросы про реализацию

Про регуляризацию β. Если факторов много, оценка по каждому активу шумная. Ridge с маленьким α снижает дисперсию оценок без сильной смещённости. Его удобно включать как «fallback» при плохом условии матрицы X.

Про нестабильность факторов. SMB и HML конструируются предсказуемо, но их временные профили меняются. Fama‑French добавили прибыльность и инвестиции в пятифакторную модель — в некоторых выборках HML становится избыточным. Не удивляйтесь, если трёхфакторка и пятифакторка дадут разные λ.

Про momentum. Carhart добавляет UMD и нередко заметно улучшает описание доходностей в акциях — особенно в коротких горизонтах. Но аккуратнее с пересечениями данных и повторным использованем информации.

Про статистические факторы. PCA‑модель часто выигрывает в прогнозе риска, когда нужна именно ковариация, а не экономическая интерпретация. Для предсказательной точности полезно валидировать число компонент out‑of‑sample.

Тесты, без которых дальше нельзя

• HAC Newey‑West на первой стадии и адекватные SE на второй.
• GRS на нулевые альфы портфелей‑тестов.
• Множественная проверка гипотез для новых факторов, пороги повыше. В духе Harvey‑Liu‑Zhu и их обновлений.

Пример минимальной проверки в коде

Склеим простую диагностику: t‑статы λ, доля объяснённой дисперсии и GRS‑проверка альф. GRS придётся реализовать по формулам или взять готовую имплементацию из учебных конспектов.

def fama_macbeth_summary(fm: FactorModel, lam_t_stack: np.ndarray) -> pd.DataFrame:     T = lam_t_stack.shape[0]     lam_mean = lam_t_stack.mean(axis=0)     lam_std  = lam_t_stack.std(axis=0, ddof=1)     t_naive = lam_mean / (lam_std / np.sqrt(T))     return pd.DataFrame({"lambda": lam_mean, "t": t_naive}, index=fm.lambdas.index)  # Заглушка для GRS: верните F-статистику и p-value по формулам Gibbons-Ross-Shanken def grs_test(alphas: np.ndarray, Sigma_e: np.ndarray, Sigma_f: np.ndarray, T: int):     # Реализация опущена для краткости: смотрите учебные заметки или статьи с формулой.     # Важно: используйте конечновыборочную поправку.     return {"F": np.nan, "p": np.nan}

Итог

Линейная регрессия — рабочая лошадь факторных моделей. Но в 2025 году просто OLS — недостаточно. Нужны устойчивые ковариации по Newey‑West, корректные кросс‑секционные ошибки и дисциплина множественных гипотез, особенно если вы пытаетесь «открывать факторы». В risk‑модели неизбежны shrinkage‑оценки и регуляризация. А если интерпретация не критична, то просто та же PCA‑модель часто даст более стабильный прогноз ковариации.

Таким образом, регрессия остаётся базовым инструментом при построении факторных моделей и оценке рисков, но её практическое применение требует аккуратной настройки и понимания ограничений. Чтобы глубже разобраться, как методы анализа данных и машинного обучения применяются именно в финансовой сфере, можно подключиться к бесплатным открытым урокам курса «ML для финансового анализа».

Ближайшие встречи:

27 августа в 20:00 — «Инструменты тестирования торговых стратегий»

4 сентября в 20:00 — «Введение в технический анализ: построение торговой стратегии»

17 сентября в 20:00 — «Работа с торговой площадкой ByBit»

Кроме того, вы можете пройти бесплатное вступительное тестирование. Оно создано для проверки ваших знаний и навыков и не связано с самим курсом.

Также рекомендуем заглянуть в секцию с отзывами — там вы найдёте мнения участников о формате и содержании занятий.

ссылка на оригинал статьи https://habr.com/ru/articles/939504/