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

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

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

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

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

1. Выбор факторов и тест‑активов. В примерах берут Fama‑French: MKT‑RF, SMB, HML, а ещё часто добавляют momentum из Carhart.

2. Тайм‑серийная регрессия для бета‑нагрузок βᵢ. Классика: для каждого актива i оцениваем
   rᵢₜ − r_fₜ = αᵢ + βᵢ⊤ fₜ + εᵢₜ,
   где fₜ это вектор факторных доходностей.

3. Кросс‑секционная регрессия премий λₜ (Fama‑MacBeth): в каждый момент времени t решаем
   r̄ₜ = 1·γ₀ₜ + B·λₜ + uₜ,
   где r̄ₜ это вектор доходностей активов в момент t, B матрица ранее оценённых β. Потом усредняем λₜ и считаем t‑статистики.

4. Построение 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)

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

1. Мультиколлинеарность факторов. Симптом: нестабильные β и взрывающиеся стандартные ошибки. Лечится регуляризацией (ridge), исключением лишних факторов, PCA‑предобработкой.

2. t‑статистика «2.0 достаточно?» — уже давно нет. Литература про факторный зоопарк убеждает смотреть на множественную проверку гипотез и поднимать пороги значимости, часто до |t| ≥ 3. Иначе ловите публикационный шум.

3. Стандартные ошибки. Для тайм‑серий — HAC Newey‑West. В Fama‑MacBeth используйте либо «по времени» с HAC, либо применяйте шанкeновскую поправку на ошибочность β.

4. GRS‑тест на нулевые альфы портфелей‑тестов. Если отклоняет, модель не объясняет кросс‑секцию.

5. Данные факторов. Поддерживайте воспроизводимость: фиксируйте конструкторы 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»

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

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