Когда я впервые столкнулся с байесовскими моделями пространственно-временных рядов в контексте разработки торговых алгоритмов, меня поразила их элегантность и мощь. В отличие от традиционных методов прогнозирования временных рядов, которые часто игнорируют неопределенность или делают слишком жесткие предположения о структуре данных, байесовские state-space модели предоставляют гибкий фреймворк для работы с динамическими системами в условиях неопределенности.
Данная технология стала краеугольным камнем современного количественного анализа, особенно в хедж-фонд индустрии, где точность прогнозов и оценка рисков определяют успех стратегий. Понимание байесовских моделей открывает совершенно новые горизонты в анализе временных рядов.
Философия байесовского подхода к временным рядам
Традиционные методы анализа временных рядов, такие как ARIMA или экспоненциальное сглаживание, работают с фиксированными параметрами и предполагают стационарность процессов. В реальном мире финансовых рынков такие предположения редко выполняются. Волатильность меняется во времени, корреляции между активами эволюционируют, а структурные сдвиги происходят регулярно.
Байесовский подход радикально меняет парадигму анализа. Вместо поиска «истинных» значений параметров, мы работаем с распределениями вероятностей, которые отражают нашу неопределенность. Каждое новое наблюдение обновляет наши убеждения согласно теореме Байеса, создавая адаптивную систему обучения.
State-space models добавляют еще один уровень гибкости. Наблюдаемые данные рассматриваются как зашумленные проекции скрытого состояния системы, которое эволюционирует во времени согласно определенной динамике. Это позволяет естественным образом моделировать латентные факторы, такие как изменения в рыночном режиме или скрытые драйверы волатильности.
Мне нравится аналогия с навигацией корабля в тумане. Капитан не может видеть истинное положение корабля (скрытое состояние), но получает зашумленные сигналы от различных приборов (наблюдения). Байесовская фильтрация позволяет оптимально объединять эту информацию для получения наилучшей оценки текущего положения и прогноза будущего движения.
Математический фундамент
Байесовская state-space модель состоит из двух основных уравнений. Уравнение состояния описывает эволюцию скрытого состояния системы:
x_t = f(x_{t-1}, θ) + w_t
- где x_t — вектор состояния в момент времени t,
- f() — функция перехода состояния,
- θ — параметры модели,
- w_t — шум процесса.
Уравнение наблюдения связывает скрытое состояние с наблюдаемыми данными:
y_t = h(x_t, φ) + v_t
- где y_t — вектор наблюдений,
- h() — функция наблюдения,
- φ — параметры наблюдения,
- v_t — шум наблюдений.
Ключевое преимущество этого представления заключается в его универсальности. Функции f() и h() могут быть линейными или нелинейными, параметры могут меняться во времени, а распределения шумов — не обязательно гауссовскими. Это открывает возможности для моделирования сложных динамических процессов, характерных для финансовых временных рядов.
Архитектура и компоненты модели
Состояние системы и его эволюция
Центральным элементом байесовской state-space модели является вектор состояния, который содержит всю релевантную информацию о системе в данный момент времени. В контексте финансовых приложений состояние может включать различные компоненты: текущий уровень цены, тренд, сезонные компоненты, волатильность, корреляции между активами.
Проектирование вектора состояния требует глубокого понимания моделируемого процесса. Слишком простое состояние не сможет захватить важные динамические особенности, в то время как излишне сложное приведет к переобучению и вычислительным проблемам. В моей работе с валютными парами я обычно начинаю с трехкомпонентного состояния: уровень, тренд и волатильность, постепенно добавляя дополнительные элементы по мере необходимости.
Эволюция состояния моделируется через матрицу перехода или более сложные нелинейные функции. Для многих финансовых приложений хорошо работают модели случайного блуждания для уровня, AR(1) процессы для тренда и GARCH-подобные динамики для волатильности. Однако байесовский подход позволяет гибко адаптировать эти компоненты под специфику конкретных данных.
Механизм наблюдений и измерений
Уравнение наблюдения определяет, как скрытое состояние проявляется в наблюдаемых данных. В простейшем случае наблюдения могут быть прямыми измерениями состояния с добавлением шума. Однако в реальных приложениях связь часто более сложная.
Например, при моделировании котировок акций наблюдаемые цены могут зависеть не только от истинной стоимости актива, но и от микроструктурных эффектов, таких как bid-ask спред, влияние крупных ордеров, или временные дисбалансы спроса и предложения. Байесовская модель позволяет явно учесть эти факторы в уравнении наблюдения.
Особенно интересны случаи с пропущенными данными или нерегулярными наблюдениями. В отличие от классических методов, которые требуют предварительной интерполяции, байесовские state-space модели естественным образом работают с неполными данными, обновляя апостериорные распределения только при наличии новых наблюдений.
Априорные распределения и их выбор
Выбор априорных распределений — это искусство, требующее баланса между экспертным знанием и гибкостью модели. Слишком информативные априорные распределения могут ограничить способность модели адаптироваться к данным, в то время как неинформативные могут привести к неустойчивости и медленной сходимости.
В финансовых приложениях я часто использую слабо информативные априорные распределения, которые отражают базовые экономические принципы, не будучи слишком ограничивающими. Например, для параметров волатильности естественно использовать log-normal или inverse-gamma распределения, которые обеспечивают положительность, но позволяют широкий диапазон значений.
Иерархические априорные распределения особенно полезны при работе с панельными данными или портфелями активов. Они позволяют моделировать как индивидуальные особенности каждого актива, так и общие рыночные факторы, создавая естественный механизм регуляризации через частичное объединение (partial pooling).
Алгоритмы фильтрации и сглаживания
Фильтр Калмана и его байесовские расширения
Фильтр Калмана (Kalman Filter) представляет собой точное решение байесовской фильтрации для линейных гауссовских state-space моделей. Его элегантность заключается в том, что он дает замкнутые формулы для обновления апостериорных распределений, делая вычисления быстрыми и численно стабильными.
Алгоритм состоит из двух этапов: предсказания и обновления. На этапе предсказания используется модель эволюции состояния для получения априорного распределения на следующем временном шаге. Этап обновления объединяет это априорное распределение с новым наблюдением согласно теореме Байеса.
Однако линейность и гауссовость — серьезные ограничения для многих практических приложений. Расширенный фильтр Калмана (EKF) адресует проблему нелинейности через линеаризацию вокруг текущей оценки состояния. Unscented фильтр Калмана (UKF) использует более сложную sigma-point трансформацию для лучшего захвата нелинейных эффектов.
В моей практике UKF часто показывает лучшие результаты при моделировании опционных цен или сложных деривативов, где нелинейности играют критическую роль. Вычислительная сложность при этом конечно возрастает, однако улучшение точности обычно оправдывает дополнительные затраты.
Particle фильтрация для сложных моделей
Когда предположения о гауссовости или даже параметрическая форма распределений неприемлемы, particle фильтры предоставляют гибкое решение. Они представляют апостериорные распределения через набор взвешенных частиц (samples), каждая из которых представляет возможную траекторию скрытого состояния.
Алгоритм particle фильтра состоит из трех основных шагов: sampling, weighting, и resampling. На каждом временном шаге генерируются новые частицы согласно модели эволюции состояния, затем их веса обновляются на основе правдоподобия наблюдений. Этап resampling предотвращает вырождение фильтра, когда большинство частиц имеют пренебрежимо малые веса.
Основная проблема particle фильтров — это «проклятие размерности» (curse of dimensionality). Количество частиц, необходимое для точного представления распределения, экспоненциально растет с размерностью пространства состояний. Однако для умеренно размерных задач (до 10-15 переменных состояния) современные вычислительные ресурсы позволяют получать хорошие результаты.
Я часто использую particle фильтры для моделирования режимных переключений в рыночных данных, где распределения состояний естественно многомодальны. Например, при анализе волатильности опционов различные рыночные режимы (низкая/высокая волатильность, trending/ranging рынки) могут быть представлены отдельными кластерами частиц.
Сглаживание и ретроспективный анализ
Фильтрация дает оптимальные оценки состояния на основе информации, доступной до текущего момента времени. Однако для многих аналитических задач полезно использовать всю доступную информацию, включая будущие наблюдения. Это задача сглаживания (smoothing).
Алгоритм RTS (Rauch-Tung-Striebel) сглаживания для линейных гауссовских моделей работает в обратном направлении времени, начиная с последнего наблюдения и используя уже вычисленные фильтрованные оценки. Для каждого временного шага сглаженная оценка объединяет фильтрованную оценку с информацией из будущих наблюдений.
Сглаживание особенно важно при работе с структурными изменениями в данных. Например, при обнаружении моментов изменения рыночного режима ретроспективный анализ позволяет более точно идентифицировать точки переключения и лучше понять динамику переходных периодов.
В практических приложениях я использую сглаживание для построения исторических оценок скрытых факторов риска, анализа эффективности торговых стратегий в различных рыночных условиях, и калибровки параметров модели на исторических данных.
Практическая реализация на Python
Базовая структура байесовской state-space модели
Для демонстрации концепций рассмотрим реализацию простой байесовской модели локального уровня с изменяющейся во времени волатильностью. Эта модель особенно полезна для анализа финансовых временных рядов, где волатильность кластеризуется во времени.
# !pip install filterpy
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import norm, invgamma
from scipy.optimize import minimize
import yfinance as yf
from filterpy.kalman import KalmanFilter
from sklearn.metrics import mean_squared_error
import warnings
warnings.filterwarnings('ignore')
class BayesianStateSpaceModel:
"""
Байесовская модель пространственно-временного ряда с изменяющейся волатильностью
"""
def __init__(self, n_states=2, n_obs=1):
"""
Инициализация модели
Parameters:
-----------
n_states : int, количество переменных состояния
n_obs : int, количество наблюдаемых переменных
"""
self.n_states = n_states
self.n_obs = n_obs
self.fitted = False
# Параметры модели
self.params = {
'sigma_level': 1.0, # Шум уровня
'sigma_slope': 0.1, # Шум тренда
'sigma_obs': 1.0, # Шум наблюдений
'phi': 0.95, # AR коэффициент для тренда
}
# Априорные распределения
self.priors = {
'sigma_level': {'alpha': 2, 'beta': 1},
'sigma_slope': {'alpha': 2, 'beta': 0.1},
'sigma_obs': {'alpha': 2, 'beta': 1},
'phi': {'mean': 0.9, 'std': 0.1}
}
def _setup_matrices(self):
"""Настройка матриц состояния"""
# Матрица перехода состояния (level + slope model)
self.F = np.array([
[1, 1],
[0, self.params['phi']]
])
# Матрица наблюдения
self.H = np.array([[1, 0]])
# Ковариационная матрица шума процесса
self.Q = np.array([
[self.params['sigma_level']**2, 0],
[0, self.params['sigma_slope']**2]
])
# Ковариационная матрица шума наблюдений
self.R = np.array([[self.params['sigma_obs']**2]])
def log_likelihood(self, y, params=None):
"""
Вычисление логарифма правдоподобия
"""
if params is not None:
self.params.update(params)
self._setup_matrices()
# Инициализация фильтра Калмана
kf = KalmanFilter(dim_x=self.n_states, dim_z=self.n_obs)
kf.F = self.F.copy()
kf.H = self.H.copy()
kf.Q = self.Q.copy()
kf.R = self.R.copy()
# Начальные условия
kf.x = np.array([y[0], 0.0])
kf.P = np.eye(self.n_states) * 100
log_lik = 0
self.filtered_states = []
self.predicted_states = []
for i, obs in enumerate(y):
# Предсказание
kf.predict()
self.predicted_states.append(kf.x.copy())
# Обновление
kf.update(np.array([obs]))
self.filtered_states.append(kf.x.copy())
# Вычисление логарифма правдоподобия
innovation = obs - np.dot(kf.H, kf.x_prior)
S = np.dot(kf.H, np.dot(kf.P_prior, kf.H.T)) + kf.R
log_lik += -0.5 * (np.log(2 * np.pi * S) + innovation**2 / S)
return log_lik
def log_prior(self):
"""Вычисление логарифма априорной плотности"""
log_prior = 0
# Log-prior для sigma параметров (inverse gamma)
for param in ['sigma_level', 'sigma_slope', 'sigma_obs']:
alpha = self.priors[param]['alpha']
beta = self.priors[param]['beta']
log_prior += invgamma.logpdf(self.params[param], alpha, scale=beta)
# Log-prior для phi (truncated normal)
if 0 < self.params['phi'] < 1:
log_prior += norm.logpdf(self.params['phi'],
self.priors['phi']['mean'],
self.priors['phi']['std'])
else:
log_prior = -np.inf
return log_prior
def log_posterior(self, y, params=None):
"""Логарифм апостериорной плотности"""
if params is not None:
temp_params = self.params.copy()
self.params.update(params)
log_lik = self.log_likelihood(y)
log_prior = self.log_prior()
if params is not None:
self.params = temp_params
return log_lik + log_priorЭта реализация демонстрирует ключевые компоненты байесовской state-space модели. Класс инкапсулирует матрицы состояния, параметры модели и априорные распределения. Метод log_likelihood использует фильтр Калмана для вычисления правдоподобия данных, а log_posterior объединяет правдоподобие с априорными убеждениями.
Особенность этой реализации — гибкая параметризация через словарь params, что позволяет легко экспериментировать с различными спецификациями модели. Априорные распределения заданы через параметры соответствующих семейств распределений, что обеспечивает прозрачность байесовского вывода.
Использование filterpy библиотеки для фильтра Калмана обеспечивает численную стабильность и производительность, особенно важные при работе с реальными финансовыми данными. Метод автоматически отслеживает отфильтрованные и предсказанные состояния, что полезно для последующего анализа и диагностики модели.
Продвинутые методы байесовского вывода
Для более сложных моделей, где фильтр Калмана неприменим, необходимы продвинутые методы байесовского вывода. Рассмотрим реализацию Metropolis-Hastings алгоритма для сэмплирования из апостериорного распределения:
import pymc as pm
import numpy as np
import matplotlib.pyplot as plt
import arviz as az
from scipy import stats
class AdvancedBayesianSSM:
"""
Продвинутая байесовская state-space модель с MCMC выводом
"""
def __init__(self):
self.model = None
self.trace = None
self.fitted = False
def build_model(self, y, stochastic_volatility=True):
"""
Построение PyMC модели
Parameters:
-----------
y : array-like, наблюдаемые данные
stochastic_volatility : bool, использовать стохастическую волатильность
"""
n_obs = len(y)
with pm.Model() as model:
# Гиперпараметры
sigma_level = pm.HalfCauchy('sigma_level', beta=2.5)
sigma_slope = pm.HalfCauchy('sigma_slope', beta=1.0)
if stochastic_volatility:
# Стохастическая волатильность
sigma_vol = pm.HalfCauchy('sigma_vol', beta=1.0)
phi_vol = pm.Beta('phi_vol', alpha=9, beta=1) # ~ 0.9
# Логарифм волатильности как AR(1) процесс
log_vol_init = pm.Normal('log_vol_init', 0, 1)
log_vol_innovations = pm.Normal('log_vol_innovations', 0, sigma_vol, shape=n_obs-1)
# Построение AR(1) процесса
log_vol = pm.Deterministic('log_vol',
pm.math.concatenate([
[log_vol_init],
log_vol_init * phi_vol + pm.math.cumsum(log_vol_innovations * pm.math.sqrt(1 - phi_vol**2))
]))
vol = pm.Deterministic('vol', pm.math.exp(log_vol))
else:
vol = pm.HalfCauchy('vol', beta=2.5)
# Состояние: уровень и наклон как случайные блуждания
level_init = pm.Normal('level_init', y[0], 10)
level_innovations = pm.Normal('level_innovations', 0, sigma_level, shape=n_obs-1)
level = pm.Deterministic('level',
pm.math.concatenate([
[level_init],
level_init + pm.math.cumsum(level_innovations)
]))
slope_init = pm.Normal('slope_init', 0, 1)
slope_innovations = pm.Normal('slope_innovations', 0, sigma_slope, shape=n_obs-1)
slope = pm.Deterministic('slope',
pm.math.concatenate([
[slope_init],
slope_init + pm.math.cumsum(slope_innovations)
]))
# Латентный процесс
mu = level + slope
# Наблюдения
if stochastic_volatility:
obs = pm.Normal('obs', mu=mu, sigma=vol, observed=y)
else:
obs = pm.Normal('obs', mu=mu, sigma=vol, observed=y)
self.model = model
return model
def fit(self, y, draws=500, tune=200, chains=2, cores=None):
"""
Фиттинг модели через MCMC
"""
if self.model is None:
self.build_model(y, stochastic_volatility=False) # Упрощенная модель без стохастической волатильности
print(f"Запуск MCMC: {draws} draws, {tune} tune, {chains} chains")
with self.model:
# MCMC сэмплирование с прогресс-баром
self.trace = pm.sample(
draws=draws,
tune=tune,
chains=chains,
cores=cores,
return_inferencedata=True,
target_accept=0.85, # Компромисс между скоростью и качеством
random_seed=42,
progressbar=True,
idata_kwargs={"log_likelihood": False}, # Отключаем log_likelihood для экономии памяти
init="auto" # Автоматическая инициализация
)
print("MCMC завершен!")
# Проверка качества сходимости
try:
rhat_stats = az.rhat(self.trace)
# Извлекаем максимальное значение R-hat из всех переменных
max_rhat = 1.0
for var_name in rhat_stats.data_vars:
var_rhat = float(rhat_stats[var_name].max().values)
max_rhat = max(max_rhat, var_rhat)
if max_rhat > 1.1:
print(f"⚠️ Предупреждение: max R-hat = {max_rhat:.3f} > 1.1 (плохая сходимость)")
print("Рекомендуется увеличить количество итераций или изменить параметры")
else:
print(f"✅ max R-hat = {max_rhat:.3f} (хорошая сходимость)")
except Exception as e:
print(f"⚠️ Не удалось вычислить R-hat: {e}")
print("Но модель обучена, можно продолжать")
self.fitted = True
return self.trace
def predict(self, steps_ahead=10, samples=1000):
"""
Генерация прогнозов
"""
if not self.fitted:
raise ValueError("Модель должна быть обучена перед прогнозированием")
# Извлечение последних состояний
last_level = self.trace.posterior['level'][:, :, -1].values.flatten()
last_slope = self.trace.posterior['slope'][:, :, -1].values.flatten()
# Параметры для прогнозирования
sigma_level = self.trace.posterior['sigma_level'].values.flatten()
sigma_slope = self.trace.posterior['sigma_slope'].values.flatten()
# Обработка волатильности (может быть скаляром или вектором)
if 'vol' in self.trace.posterior:
vol_data = self.trace.posterior['vol'].values
if vol_data.ndim > 2: # Если это временной ряд
vol = vol_data[:, :, -1].flatten()
else: # Если это скаляр
vol = vol_data.flatten()
else:
vol = np.ones_like(last_level) # fallback
# Случайная выборка из апостериорного распределения
sample_idx = np.random.choice(len(last_level), samples, replace=True)
predictions = []
for i in range(samples):
idx = sample_idx[i]
level = last_level[idx]
slope = last_slope[idx]
s_level = sigma_level[idx]
s_slope = sigma_slope[idx]
v = vol[idx] if hasattr(vol, '__len__') and len(vol) > idx else (vol[0] if hasattr(vol, '__len__') else vol)
forecast = []
for step in range(steps_ahead):
# Эволюция состояния
level += np.random.normal(0, s_level)
slope += np.random.normal(0, s_slope)
# Прогноз наблюдения
pred = level + slope + np.random.normal(0, v)
forecast.append(pred)
predictions.append(forecast)
predictions = np.array(predictions)
# Статистики прогнозов
pred_mean = np.mean(predictions, axis=0)
pred_std = np.std(predictions, axis=0)
pred_quantiles = np.percentile(predictions, [5, 25, 50, 75, 95], axis=0)
return {
'mean': pred_mean,
'std': pred_std,
'quantiles': pred_quantiles,
'samples': predictions
}
def plot_diagnostics(self):
"""Диагностические графики"""
if not self.fitted:
raise ValueError("Модель должна быть обучена")
# Trace plots
var_names = ['sigma_level', 'sigma_slope']
if 'vol' in self.trace.posterior:
var_names.append('vol')
if 'phi_vol' in self.trace.posterior:
var_names.append('phi_vol')
az.plot_trace(self.trace, var_names=var_names)
plt.tight_layout()
plt.show()
# Posterior plots
az.plot_posterior(self.trace, var_names=var_names)
plt.tight_layout()
plt.show()
# Autocorrelation
az.plot_autocorr(self.trace, var_names=['sigma_level', 'sigma_slope'])
plt.tight_layout()
plt.show()
def quick_fit(self, y, draws=100, tune=50):
"""
Быстрый фит для тестирования (низкое качество, но быстро)
"""
if self.model is None:
self.build_model(y, stochastic_volatility=False) # Упрощенная модель
print(f"Быстрый режим: {draws} draws, {tune} tune")
with self.model:
self.trace = pm.sample(
draws=draws,
tune=tune,
chains=1, # Только одна цепь
cores=1,
return_inferencedata=True,
target_accept=0.7,
random_seed=42,
progressbar=True,
discard_tuned_samples=True
)
print("Быстрый фит завершен!")
self.fitted = True
return self.trace
def summary(self):
"""Краткая сводка результатов"""
if not self.fitted:
raise ValueError("Модель должна быть обучена")
return az.summary(self.trace)
# Пример использования
if __name__ == "__main__":
# Генерация тестовых данных
np.random.seed(42)
n = 50
true_level = np.cumsum(np.random.normal(0, 0.1, n))
true_slope = np.cumsum(np.random.normal(0, 0.05, n))
noise = np.random.normal(0, 0.2, n)
y = true_level + true_slope + noise
print("Начинаем обучение модели...")
print(f"Размер данных: {len(y)} точек")
# Создание и обучение модели
model = AdvancedBayesianSSM()
# Выберите один из вариантов:
# 1. Быстрый тест:
# model.quick_fit(y)
# 2. Стандартный режим:
trace = model.fit(y, draws=500, tune=200, chains=2)
print("Обучение завершено!")
# Прогнозирование
predictions = model.predict(steps_ahead=10, samples=200)
# Построение графиков
plt.figure(figsize=(12, 8))
# Исходные данные
plt.subplot(2, 1, 1)
plt.plot(y, 'b-', label='Наблюдения')
plt.plot(true_level + true_slope, 'r--', alpha=0.7, label='Истинный тренд')
plt.legend()
plt.title('Исходные данные')
# Прогнозы
plt.subplot(2, 1, 2)
future_x = np.arange(len(y), len(y) + len(predictions['mean']))
plt.plot(y, 'b-', label='Наблюдения')
plt.plot(future_x, predictions['mean'], 'g-', label='Прогноз (среднее)')
plt.fill_between(future_x,
predictions['quantiles'][0],
predictions['quantiles'][4],
alpha=0.3, color='green', label='90% доверительный интервал')
plt.legend()
plt.title('Прогнозы')
plt.tight_layout()
plt.show()
# Диагностика (опционально)
# model.plot_diagnostics()
# Сводка
try:
print("\n=== СВОДКА РЕЗУЛЬТАТОВ ===")
summary_df = model.summary()
print(summary_df)
# Дополнительная диагностика
print(f"\nКоличество параметров: {len(summary_df)}")
# Проверка эффективного размера выборки
try:
ess_stats = az.ess(model.trace)
min_ess = float('inf')
for var_name in ess_stats.data_vars:
var_ess = float(ess_stats[var_name].min().values)
min_ess = min(min_ess, var_ess)
if min_ess == float('inf'):
min_ess = 0
print(f"Минимальный ESS: {min_ess:.0f}")
if min_ess < 100:
print("⚠️ ESS < 100: рекомендуется больше итераций")
else:
print("✅ ESS достаточный")
except Exception as e:
print(f"Не удалось вычислить ESS: {e}")
except Exception as e:
print(f"Ошибка при выводе сводки: {e}")
print("Но модель обучена успешно!")
# Раскомментируйте для мгновенного теста:
# quick_demo()Начинаем обучение модели...
Размер данных: 50 точек
Запуск MCMC: 500 draws, 200 tune, 2 chains
Progress Draws Divergences Step size Grad evals Sampling Speed Elapsed Remaining
─────────────────────────────────────────────────────────────────────────────────────────────────────────────────
━━━━━━━━━━━━━━━━━━━━━━━━━ 700 2 0.02 255 32.13 draws/s 0:00:21 0:00:00
━━━━━━━━━━━━━━━━━━━━━━━━━ 700 0 0.02 255 16.26 draws/s 0:00:42 0:00:00
MCMC завершен!
⚠️ Предупреждение: max R-hat = 1.437 > 1.1 (плохая сходимость)
Рекомендуется увеличить количество итераций или изменить параметры
Обучение завершено!
Рис. 1: Визуализация исходных данных и байесовского прогноза (зеленая линия с 90% доверительным интервалом) будущих значений временного ряда
=== СВОДКА РЕЗУЛЬТАТОВ ===
mean sd hdi_3% hdi_97% mcse_mean mcse_sd \
level[0] -0.054 0.991 -1.805 1.872 0.039 0.041
level[1] -0.017 0.990 -1.782 1.913 0.039 0.040
level[2] 0.030 0.992 -1.894 1.827 0.039 0.039
level[3] 0.086 0.989 -1.758 1.915 0.039 0.040
level[4] 0.165 0.991 -1.610 2.147 0.040 0.039
... ... ... ... ... ... ...
slope_innovations[45] -0.039 0.092 -0.254 0.093 0.020 0.020
slope_innovations[46] 0.002 0.078 -0.188 0.149 0.002 0.019
slope_innovations[47] 0.028 0.086 -0.119 0.206 0.013 0.018
slope_innovations[48] -0.006 0.075 -0.171 0.130 0.002 0.017
vol 0.204 0.032 0.150 0.269 0.002 0.001
ess_bulk ess_tail r_hat
level[0] 659.0 575.0 1.00
level[1] 662.0 589.0 1.00
level[2] 651.0 565.0 1.00
level[3] 639.0 609.0 1.00
level[4] 616.0 566.0 1.00
... ... ... ...
slope_innovations[45] 36.0 75.0 1.17
slope_innovations[46] 1369.0 133.0 1.17
slope_innovations[47] 144.0 197.0 1.17
slope_innovations[48] 1730.0 184.0 1.16
vol 277.0 267.0 1.01
[203 rows x 9 columns]
Количество параметров: 203
Минимальный ESS: 4
⚠️ ESS < 100: рекомендуется больше итерацийЭта продвинутая реализация использует PyMC для полного байесовского вывода через MCMC. Ключевые преимущества этого подхода включают:
- автоматическое получение неопределенности в параметрах;
- возможность работы с произвольными априорными распределениями;
- естественную интеграцию стохастической волатильности.
Метод build_model демонстрирует гибкость PyMC в спецификации сложных иерархических моделей. Реализация случайных блужданий для состояний через кумулятивные суммы обеспечивает правильную факторизацию апостериорного распределения. Опциональная стохастическая волатильность моделируется через AR(1) процесс в логарифмическом пространстве, что гарантирует положительность волатильности.
Метод прогнозирования predict генерирует полные апостериорные распределения прогнозов, включая квантили неопределенности. Это критически важно для риск-менеджмента и принятия торговых решений. Сэмплирование из апостериорного распределения параметров обеспечивает правильное распространение неопределенности в прогнозы.
Применение к реальным финансовым данным
Продемонстрируем применение разработанных моделей к анализу реальных рыночных данных:
import yfinance as yf
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import pymc as pm
import arviz as az
# Загрузка данных
def load_market_data(symbol='SPY', period='2y'):
"""Загрузка рыночных данных"""
ticker = yf.Ticker(symbol)
data = ticker.history(period=period)
# Логарифмические доходности
returns = np.log(data['Close'] / data['Close'].shift(1)).dropna()
# Кумулятивные логарифмические цены
log_prices = np.log(data['Close'])
return {
'prices': data['Close'],
'log_prices': log_prices,
'returns': returns,
'volume': data['Volume'],
'dates': data.index
}
# Простая модель для вычисления правдоподобия
class SimpleBayesianModel:
def __init__(self):
pass
def log_likelihood(self, data):
"""Простое вычисление логарифма правдоподобия"""
# Предполагаем нормальное распределение
mu = np.mean(data)
sigma = np.std(data)
# Логарифм правдоподобия для нормального распределения
log_lik = -0.5 * len(data) * np.log(2 * np.pi * sigma**2)
log_lik -= 0.5 * np.sum((data - mu)**2) / sigma**2
return log_lik
# Продвинутая модель из предыдущего кода
class AdvancedBayesianSSM:
def __init__(self):
self.model = None
self.trace = None
self.fitted = False
def build_model(self, y, stochastic_volatility=False):
n_obs = len(y)
with pm.Model() as model:
sigma_level = pm.HalfNormal('sigma_level', sigma=0.01)
sigma_slope = pm.HalfNormal('sigma_slope', sigma=0.005)
vol = pm.HalfNormal('vol', sigma=0.02) # Константная волатильность
# Центрированная параметризация для лучшей сходимости
level_raw = pm.Normal('level_raw', 0, 1, shape=n_obs)
slope_raw = pm.Normal('slope_raw', 0, 1, shape=n_obs)
# Некоррелированные инновации
level_innovations = sigma_level * level_raw[1:]
slope_innovations = sigma_slope * slope_raw[1:]
# Построение состояний
level = pm.Deterministic('level',
pm.math.concatenate([
[y[0] + sigma_level * level_raw[0]], # Начальный уровень
y[0] + sigma_level * level_raw[0] + pm.math.cumsum(level_innovations)
]))
slope = pm.Deterministic('slope',
pm.math.concatenate([
[sigma_slope * slope_raw[0]], # Начальный наклон
sigma_slope * slope_raw[0] + pm.math.cumsum(slope_innovations)
]))
# Латентный процесс (упрощенный)
mu = level
# Наблюдения
obs = pm.Normal('obs', mu=mu, sigma=vol, observed=y)
self.model = model
return model
def fit(self, y, draws=1000, tune=500, chains=4, cores=None):
"""Фиттинг модели"""
if self.model is None:
self.build_model(y, stochastic_volatility=False)
print(f"Запуск MCMC: {draws} draws, {tune} tune, {chains} chains")
with self.model:
self.trace = pm.sample(
draws=draws,
tune=tune,
chains=chains,
cores=cores,
return_inferencedata=True,
target_accept=0.95,
random_seed=42,
progressbar=True,
idata_kwargs={"log_likelihood": True},
init="adapt_diag" # Адаптивная диагональная инициализация
)
print("MCMC завершен!")
# Проверка дивергенций
try:
divergences = self.trace.sample_stats.diverging.sum().values
total_divergences = int(divergences) if hasattr(divergences, 'item') else int(divergences.sum())
if total_divergences > 0:
print(f"⚠️ {total_divergences} дивергенций обнаружено")
else:
print("✅ Дивергенций не обнаружено")
except:
print("Не удалось проверить дивергенции")
self.fitted = True
return self.trace
def predict(self, steps_ahead=10, samples=500):
"""Генерация прогнозов"""
if not self.fitted:
raise ValueError("Модель должна быть обучена")
# Извлечение последних состояний
last_level = self.trace.posterior['level'][:, :, -1].values.flatten()
last_slope = self.trace.posterior['slope'][:, :, -1].values.flatten()
# Параметры
sigma_level = self.trace.posterior['sigma_level'].values.flatten()
sigma_slope = self.trace.posterior['sigma_slope'].values.flatten()
vol_data = self.trace.posterior['vol'].values
if vol_data.ndim > 2:
vol = vol_data[:, :, -1].flatten()
else:
vol = vol_data.flatten()
# Сэмплирование
sample_idx = np.random.choice(len(last_level), samples, replace=True)
predictions = []
for i in range(samples):
idx = sample_idx[i]
level = last_level[idx]
slope = last_slope[idx]
s_level = sigma_level[idx]
s_slope = sigma_slope[idx]
v = vol[idx] if hasattr(vol, '__len__') and len(vol) > idx else vol[0] if hasattr(vol, '__len__') else vol
forecast = []
for step in range(steps_ahead):
level += np.random.normal(0, s_level)
slope += np.random.normal(0, s_slope)
pred = level + slope + np.random.normal(0, v)
forecast.append(pred)
predictions.append(forecast)
predictions = np.array(predictions)
return {
'mean': np.mean(predictions, axis=0),
'std': np.std(predictions, axis=0),
'quantiles': np.percentile(predictions, [5, 25, 50, 75, 95], axis=0),
'samples': predictions
}
def plot_diagnostics(self):
"""Диагностические графики"""
if not self.fitted:
raise ValueError("Модель должна быть обучена")
var_names = ['sigma_level', 'sigma_slope', 'vol']
az.plot_trace(self.trace, var_names=var_names)
plt.tight_layout()
plt.show()
# Пример использования
print("Загрузка данных SPY...")
market_data = load_market_data('SPY', '5y')
log_prices = market_data['log_prices'].values
returns = market_data['returns'].values
print(f"Загружено {len(returns)} наблюдений")
# Фиттинг простой модели
simple_model = SimpleBayesianModel()
log_lik = simple_model.log_likelihood(returns)
print(f"Логарифм правдоподобия: {log_lik:.2f}")
# Фиттинг продвинутой модели
returns_subset = returns[-250:] # Последний год торгов
print(f"\nОбучение на {len(returns_subset)} последних наблюдениях...")
advanced_model = AdvancedBayesianSSM()
# Выберите режим обучения:
# Режим 1: Полное обучение
trace = advanced_model.fit(returns_subset, draws=1000, tune=500, chains=4)
# Режим 2: Быстрое обучение
# trace = advanced_model.fit(returns_subset, draws=300, tune=200, chains=2)
# Генерация прогнозов
print("Генерация прогнозов...")
predictions = advanced_model.predict(steps_ahead=30, samples=500)
# Визуализация результатов
fig, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(12, 10))
# График 1: Исходные цены
dates = market_data['dates']
ax1.plot(dates, market_data['prices'], 'k-', alpha=0.7, label='Цена SPY')
ax1.set_title('Цены SPY', fontsize=12)
ax1.set_ylabel('Цена ($)')
ax1.legend()
ax1.grid(True, alpha=0.3)
# График 2: Доходности
ax2.plot(dates[1:], returns, 'darkgray', alpha=0.6, label='Доходности')
ax2.axhline(y=0, color='red', linestyle='--', alpha=0.5)
ax2.set_title('Доходности SPY', fontsize=12)
ax2.set_ylabel('Доходность')
ax2.legend()
ax2.grid(True, alpha=0.3)
# График 3: Прогнозы
lookback = 60
future_steps = len(predictions['mean'])
# Исторические данные для контекста
ax3.plot(range(-lookback, 0), returns_subset[-lookback:],
'k-', alpha=0.7, label='Исторические данные')
# Прогнозы
future_x = range(0, future_steps)
pred_mean = predictions['mean']
pred_quantiles = predictions['quantiles']
ax3.plot(future_x, pred_mean, 'b-', linewidth=2, label='Прогноз (среднее)')
ax3.fill_between(future_x,
pred_quantiles[0], pred_quantiles[4],
alpha=0.2, color='blue', label='90% доверительный интервал')
ax3.fill_between(future_x,
pred_quantiles[1], pred_quantiles[3],
alpha=0.3, color='blue', label='50% доверительный интервал')
ax3.axvline(x=0, color='red', linestyle='--', alpha=0.5, label='Начало прогноза')
ax3.set_title(f'Прогноз доходностей на {future_steps} дней вперед', fontsize=12)
ax3.set_ylabel('Прогнозируемая доходность')
ax3.set_xlabel('Дни')
ax3.legend()
ax3.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# Диагностика модели
print("\nПостроение диагностических графиков...")
advanced_model.plot_diagnostics()
# Вычисление метрик качества
def calculate_metrics(trace):
"""Вычисление метрик качества модели"""
try:
# WAIC (Widely Applicable Information Criterion)
waic_result = az.waic(trace)
waic_value = float(waic_result.elpd_waic)
# Effective sample size
ess_stats = az.ess(trace)
min_ess = float('inf')
for var_name in ess_stats.data_vars:
var_ess = float(ess_stats[var_name].min().values)
min_ess = min(min_ess, var_ess)
# R-hat статистика
rhat_stats = az.rhat(trace)
max_rhat = 1.0
for var_name in rhat_stats.data_vars:
var_rhat = float(rhat_stats[var_name].max().values)
max_rhat = max(max_rhat, var_rhat)
return {
'WAIC': waic_value,
'min_ESS': min_ess,
'max_R_hat': max_rhat,
'waic_full': waic_result
}
except Exception as e:
print(f"Ошибка при вычислении метрик: {e}")
return None
print("\nВычисление метрик качества...")
metrics = calculate_metrics(trace)
if metrics:
print("Метрики качества модели:")
print(f"WAIC ELPD: {metrics['WAIC']:.2f}")
print(f"Минимальный ESS: {metrics['min_ESS']:.0f}")
print(f"Максимальный R-hat: {metrics['max_R_hat']:.3f}")
# Дополнительная информация о WAIC
waic_full = metrics['waic_full']
print(f"WAIC (полный): {waic_full}")
# Интерпретация
if metrics['max_R_hat'] < 1.1:
print("✅ Хорошая сходимость (R-hat < 1.1)") else: print("⚠️ Плохая сходимость (R-hat >= 1.1)")
if metrics['min_ESS'] >= 100:
print("✅ Достаточный эффективный размер выборки")
else:
print("⚠️ Недостаточный эффективный размер выборки")
else:
print("Не удалось вычислить метрики качества")
print("\nАнализ завершен!")Загрузка данных SPY...
Загружено 1255 наблюдений
Логарифм правдоподобия: 3867.84
Обучение на 250 последних наблюдениях...
Запуск MCMC: 1000 draws, 500 tune, 4 chains
Progress Draws Divergences Step size Grad evals Sampling Speed Elapsed Remaining
─────────────────────────────────────────────────────────────────────────────────────────────────────────────────
━━━━━━━━━━━━━━━━━━━━━━━━━ 1500 0 0.11 31 136.83 draws/s 0:00:10 0:00:00
━━━━━━━━━━━━━━━━━━━━━━━━━ 1500 5 0.16 31 67.51 draws/s 0:00:22 0:00:00
━━━━━━━━━━━━━━━━━━━━━━━━━ 1500 0 0.11 31 44.10 draws/s 0:00:33 0:00:00
━━━━━━━━━━━━━━━━━━━━━━━━━ 1500 0 0.13 31 33.34 draws/s 0:00:44 0:00:00
MCMC завершен!
⚠️ 5 дивергенций обнаружено
Генерация прогнозов...
Рис. 2: Графики цен закрытия SPY, доходностей SPY, прогноза доходностей на 30 дней вперед
Рис. 3: Диагностические графики модели (sigma_level, sigma_slope, vol)
Диагностические графики показывают сходимость MCMC цепей для параметров модели:
- trace plots демонстрируют стабильность сэмплирования (цепи должны «смешиваться» без трендов);
- posterior plots показывают апостериорные распределения параметров волатильности уровня (sigma_level), наклона (sigma_slope) и наблюдений (vol), где узкие распределения указывают на хорошую идентифицируемость параметров.
Вычисление метрик качества...
Метрики качества модели:
WAIC ELPD: 725.80
Минимальный ESS: 1642
Максимальный R-hat: 1.006
WAIC (полный): Computed from 4000 posterior samples and 250 observations log-likelihood matrix.
Estimate SE
elpd_waic 725.80 38.44
p_waic 13.14 -
✅ Хорошая сходимость (R-hat < 1.1)
✅ Достаточный эффективный размер выборки
Анализ завершен!Данный пример демонстрирует полный цикл применения байесовской state-space модели к реальным финансовым данным. Код загружает данные SPY ETF, применяет модель к логарифмическим доходностям и генерирует прогнозы с полными доверительными интервалами.
Особое внимание уделено визуализации неопределенности прогнозов через квантили апостериорного распределения. Это критически важно для практических приложений, где понимание рисков не менее важно самих прогнозов. Диагностические метрики (DIC, WAIC, ESS, R-hat) позволяют оценить качество MCMC сэмплирования и сравнивать различные спецификации модели.
Реализация обеспечивает баланс между теоретической строгостью и практической применимостью. Использование современных библиотек (PyMC, ArviZ) гарантирует численную стабильность и соответствие лучшим практикам байесовского анализа.
Практические соображения и ограничения
Вычислительная сложность и масштабируемость
Основной вызов при работе с байесовскими state-space моделями — их вычислительная сложность. MCMC методы, хотя и теоретически обоснованы, могут требовать значительного времени для сходимости, особенно в высокомерных задачах. Particle фильтры страдают от проклятия размерности, требуя экспоненциально растущего числа частиц.
Вариационные методы предоставляют привлекательную альтернативу, заменяя точный байесовский вывод оптимизационной задачей. Вариационная нижняя граница (ELBO) может быть оптимизирована стохастическими градиентными методами, обеспечивая масштабируемость к большим датасетам. Однако качество аппроксимации зависит от выбора вариационного семейства распределений.
Для практических приложений критически важна параллелизация вычислений. Современные реализации (Stan, PyMC, Edward) эффективно используют многоядерные процессоры и GPU для ускорения MCMC и вариационного вывода. Распределенные вычисления позволяют масштабировать модель к портфелям с тысячами активов.
Лично я предпочитаю гибридный подход: использование быстрых приближенных методов для первичного анализа и скрининга моделей, с последующим применением точных байесовских методов к наиболее перспективным спецификациям.
Выбор модели и сравнение альтернатив
Байесовский фреймворк предоставляет принципиальные методы сравнения моделей через маргинальные правдоподобия и информационные критерии. Однако на практике вычисление маргинальных правдоподобий часто численно нестабильно, особенно для сложных моделей.
Информационные критерии (DIC, WAIC, LOO-CV) предоставляют практические альтернативы, балансируя качество подгонки с сложностью модели. Методы кросс-валидации, адаптированные для временных рядов, позволяют оценивать производительность на out-of-sample выборках без нарушения временной структуры данных.
Важно подойти с умом к спецификации априорных распределений. Слишком информативные априорные могут доминировать над данными, в то время как неинформативные могут приводить к неидентифицируемости. Вот почему анализ чувствительности к выбору априорных распределений должен стать частью любого байесовского анализа.
В практических приложениях я рекомендую начинать с простых моделей и постепенно увеличивать сложность, тщательно документируя влияние каждого дополнительного компонента на качество прогнозов и интерпретируемость результатов.
Валидация и тестирование моделей
Валидация байесовских state-space моделей требует комплексного подхода, включающего как статистические тесты, так и экономические критерии. Следует использовать:
- Апостериорные проверки прогнозируемости (Posterior predictive checks). Они позволяют оценить способность модели воспроизводить ключевые статистические свойства данных;
- Анализ остатков (Residual analysis) остается важным инструментом диагностики. Стандартизированные остатки должны демонстрировать отсутствие автокорреляции, гомоскедастичность и соответствие предполагаемому распределению;
- Q-q графики и тесты нормальности помогут выявить несоответствия в спецификации шумовых компонентов.
Для финансовых приложений критически важна оценка экономической значимости результатов. Следует учитывать, что статистически значимые улучшения в метриках качества прогнозов не всегда транслируются в экономическую прибыль после учета транзакционных издержек и рыночных ограничений.
Бэктестинг на out-of-sample данных должен имитировать реальные условия торговли, включая ограничения ликвидности, проскальзывания и издержки, связанные с исполнением ордеров. Анализ скользящих окон (Rolling window analysis) позволит оценить стабильность модели во времени и выявить периоды, когда модель работает плохо.
Заключение
Байесовские модели пространственно-временных рядов представляют собой мощный и гибкий инструментарий для анализа динамических систем в условиях неопределенности. Их способность естественным образом инкорпорировать априорную информацию, квантифицировать неопределенность и адаптироваться к изменяющимся условиям делает их особенно ценными для финансовых приложений.
Ключевые преимущества этого подхода включают принципиальную обработку неопределенности, гибкость в моделировании сложных зависимостей, и возможность инкорпорирования экспертного знания через априорные распределения. Современные вычислительные методы делают эти модели практически применимыми даже для высокомерных задач.
Вместе с тем, успешное применение Bayesian State-Space Time Series Models требует глубокого понимания как теоретических основ, так и практических аспектов реализации. Выбор подходящей архитектуры модели, спецификация априорных распределений, и валидация результатов остаются нетривиальными задачами, требующими экспертизы и опыта.
В контексте современных финансовых рынков, характеризующихся высокой волатильностью и структурными изменениями, байесовские state-space модели предоставляют интересные методы для построения адаптивных торговых систем и риск-моделей. Их способность количественно оценивать неопределенность делает их незаменимыми для принятия обоснованных инвестиционных решений в условиях ограниченной информации.