Прогнозирование и ценообразование с рядом Фурье

  • Варианты ценообразования с рядом Фурье
  • Хестон, Дисперсия-Гамма и Блэк-Шоулз
  • Прогнозирование с помощью рядов Фурье
  • Понимание преобразования Фурье в финансах
  • Преобразование Фурье и его применение в машинном обучении
  • Роль преобразования Фурье в проектировании признаков

Варианты ценообразования с рядом Фурье


Часть 1. Теория и интуиция о том, как использовать ряды Фурье для оценки деривативов.

  • Записная книжка: Github
  • Пожалуйста, ознакомьтесь с Частью 2 и Частью 3, если вы хотите увидеть больше.

Это первая статья, в которой я собираюсь помочь вам получить интуицию о том, как использовать этот новый мощный инструмент для решения полуаналитических проблем в финансах и заменить ваши методы Монте-Карло.

Знакомство

Мы все знаем и любим методы Монте-Карло для численного интегрирования, но что, если я скажу вам, что вы можете использовать мнимые числа и ряды Фурье, чтобы заменить Монте-Карло?
Основным преимуществом является скорость, которая в ценообразовании опционов очень важна. Это очень важно, так как модель Хестона для ценообразования опционов на акции должна быть численно интегрирована, что с методом Монте-Карло займет около 100 мс, а с рядом Фурье — пару мс.

Раздел 1: Но что такое ряд Фурье?

Для любой функции и для интервала a,b мы можем аппроксимировать f(x) как бесконечную сумму косинусов и синусов с L = b-a.

Раздел 2: Применение математических формул к Python

def get_fourier_approx(f, x:np.array, a:float, b:float, N:int):
fa = lambda x, n : f(x) * cos((2*pi*n*x)/(b — a))
fb = lambda x, n : f(x) * sin((2*pi*n*x)/(b — a))

A0 = 1/(b — a) * quad(f, a, b, limit=200)[0]

Cosine_Sine_Sum = np.zeros_like(x)
for n in range(1, N+1):
A = 2/(b — a) * quad(fa, a, b, args=(n), limit=200)[0]
B = 2/(b — a) * quad(fb, a, b, args=(n), limit=200)[0]
Cosine_Sine_Sum += A*cos((2*pi*n*x)/(b — a)) + B*sin((2*pi*n*x)/(b — a))

fx = A0 + Cosine_Sine_Sum
return fxa = -6
b = 6
x = np.linspace(a, b, 1_000)
y = f(x)

fig, (ax1, ax2) = plt.subplots(2, figsize=(20,12))
blue_shades = [‘#0000FF’, ‘#3399FF’, ‘#66B2FF’, ‘#99CCFF’, ‘#CCE5FF’]

avg_residuals = []
Ns = [8, 16, 32, 64, 128]
for i, N in enumerate(Ns):
fx = get_fourier_approx(f=f, x=x, a=a, b=b, N=N)
ax1.plot(x,fx, blue_shades[i], label=f’N = {N}’)
ax2.plot(x,y-fx, blue_shades[i], label=f’N = {N}’)
avg_residuals.append(np.abs(y-fx).mean())

ax1.set_title(‘Fourier Transform of f(x)’)
ax1.plot(x,y,’tab:red’, linestyle=’—‘)
ax2.set_title(‘Residuals’)
plt.tight_layout() ; ax1.legend();ax2.legend() ; plt.show()

pd.Series(avg_residuals, index=Ns, name=’Avg Residual’)

Квадратная функция:

Источник: Блокнот

N Avg. Residual
———————
8 1.311711
16 0.784683
32 0.440387
64 0.268449
128 0.154604

Линейная функция:

Источник: Блокнот

N Avg. Residual
———————
8 0.447389
16 0.264635
32 0.153540
64 0.088745
128 0.052147

Нормальное распределение

  • Масштабируется так y в [0, 12] с:
    — среднее = 100
    — std = 0,1 * sqrt (5) * 100
    — a = 100 -12 * std
    — b = 100 +12 * стандарт
Источник: Блокнот

N Avg. Residual
———————
8 1.092374e-01
16 8.326020e-05
32 6.878539e-14
64 5.721031e-14
128 5.170898e-14

Замечания

  1. Все распределения масштабируются таким образом, что y колеблется от [0,12], поэтому мы можем сравнить величину остатков.
  2. Из графиков и остатков видно, что чем более изогнута функция, тем быстрее ряды Фурье сходятся к правильному значению. Мы используем это свойство как нормальное, и логарифмическая нормаль не нуждается в большом количестве терминов для вычисления с достаточной точностью в нашем приближении.
  3. В начале и в конце данных значительно выше погрешность. Поэтому рекомендуется включать более высокие пределы, чем те, которые предназначены для использования. Например, рассчитайте ±4std, когда вам нужно ±3. Это затрудняет расчет глубоких опционов без денег.

Раздел 3: Логарифмическое нормальное распределение S(T)

S_T следует за простым GBM под Q, где мы можем вывести плотность вероятности S_T, используя следующие уравнения:

Теперь мы можем определить f(S_T) в Python с помощью следующих функций и определить нижний предел как ( 0, S_0*exp(r*T) + 12 * sigma*sqrt(T)*S_0 )S0 = 100
r = 0.05
sigma = 0.1
T = 5.0

Z = lambda St : np.where(St > 0, ((log(St/S0) — (r — .5*sigma)*T)/(sqrt(T)*sigma)), -np.inf)
f = lambda x : norm.pdf(Z(x))

a = S0*exp(r*T) — 12 * sigma*sqrt(T)*S0
b = S0*exp(r*T) + 12 * sigma*sqrt(T)*S0

Источник: Блокнот

N Avg. (scaled) Residual Avg. Residual Execution Time (sec)
—————————————————————————-
8 0.176429 5.880975e-03 0.112720
16 0.004235 1.411566e-04 0.246473
32 0.000030 9.855127e-07 0.624209
64 0.000027 8.918504e-07 1.936948
128 0.000026 8.530034e-07 6.741019

Замечания:

  1. Я включил как масштабированные, так и немасштабированные остатки. Масштабированные остатки соответствуют (неправильно) масштабированной вероятности, такой, что max{y}=12, где (правильный) немасштабированный, max{y}=0,4. Это было сделано для сравнения соответствия логарифмического нормального распределения остальным функциям, изображенным выше.
  2. Мы можем сделать вывод, что логарифмическое нормальное распределение труднее соответствовать, чем нормальное распределение из-за несимметричной формы.
  3. Мы видим, что при 64 членах в немасштабированной версии ожидаемая погрешность при вычислении P(S_T=x) очень мала и составляет менее 0,0001%.
  4. Очень важно центрировать распределение вокруг ожидаемого значения на T, с 12stds симметрично слева и справа. Я сделал версию, где a и b не симметричны, а остатки распределены неравномерно. Интуитивно понятно, что вы бы взяли a = 0+, но это не даст желаемых результатов. *Остатки S_T рассчитываются для значений S_T>0, так как это невозможно, и нас не волнуют значения меньше 0.
Плотность S_T с a=1e-8, проявляя нежелательные свойства при остатках. Источник: Блокнот

Ограничения — Недостатки — Улучшения:

  1. Нет никакой пользы в аппроксимации любой из ранее упомянутых функций в виде ряда Фурье и использовании численного интегрирования в качестве средства для вычисления An и Bn.
  2. Очень важно рассчитать коэффициенты An и Bn аналитически, поэтому единственной числовой частью является вычисление ряда.
  3. Выгода заключается в другом. Это когда f(x) не имеет явной формы и требуется численное интегрирование, и мы можем решить задачу полуаналитически с характеристической функцией и рядами Фурье.
  4. Если бы мы оценили опцион со стандартным BS в Python с помощью scipy.norm, это заняло бы около 0,06 миллисекунды. Так что использовать Фурье для стандартной БС не имеет смысла.
    Однако, если мы аналитически решим интегралы A0, An, Bn и воспользуемся версией с комплексными числами, мы получим около 0,6 миллисекунды, что сопоставимо. Мы будем использовать это в наших интересах в модели Хестона, которая является отраслевым стандартом для ценообразования таких опций, в части 3.

Часть 2

Содержание:

  • Косинусное разложение Фурье
  • Функция плотности S(T)
  • Формула косинуса Эйлера
  • Характеристическая функция — характеристическая функция геометрического броуновского движения.
  • Формулы ценообразования опционов
  • Сходимость по N членам и важность выбора хороших пределов
  • Анализ погрешностей над аналитическим решением

«Нет ничего невозможного для того, кто постарается»
Александр Македонский

Косинусное разложение Фурье

Подобно синусно-косинусному разложению, мы можем сделать то же самое только с cos-частью уравнения:

Переход от косинусно-синусной трансформации к косинусной трансформации

Из следующей функции, которую мы также сделали в части 1, вы можете видеть, что косинус способен очень хорошо аппроксимировать те же функцииdef get_fourier_approx(f, x:np.array, a:float, b:float, N:int):
fa = lambda x, n : f(x) * cos((2*pi*n*x)/(b — a))
fb = lambda x, n : f(x) * sin((2*pi*n*x)/(b — a))

A0 = 1/(b — a) * quad(f, a, b, limit=200)[0]

Cosine_Sine_Sum = np.zeros_like(x)
for n in range(1, N+1):
A = 2/(b — a) * quad(fa, a, b, args=(n), limit=200)[0]
B = 2/(b — a) * quad(fb, a, b, args=(n), limit=200)[0]
Cosine_Sine_Sum += A*cos((2*pi*n*x)/(b — a)) + B*sin((2*pi*n*x)/(b — a))

fx = A0 + Cosine_Sine_Sum
return fx

mean = 100
std = .1 *sqrt(5)*100
f = lambda x : (1/(std*sqrt(2*pi))) * exp(-(x-mean)**2/(2*std**2)) #*12/0.0175

a = mean — 12 * std
b = mean + 12 * std

Квадратная функция — Источник:Notebook
Линейная функция — Источник: Записная книжка
Нормальное распределение — Источник: Notebook

Avg Residuals
N Square f Line f Normal Distribution
8 1.311711 0.447389 1.593045e-04
16 0.784683 0.264635 1.214211e-07
32 0.440387 0.153540 7.100142e-18
64 0.268449 0.088745 3.628299e-14
128 0.154604 0.052147 2.992953e-12

Замечания

  • Выводы те же, что и в части 1, что нормальное распределение очень хорошо подходит, и нам нужно взять адекватный диапазон значений [a, b]
  • Так как используется только косинусная часть. Это упрощает многие следующие шаги, которые мы должны сделать, и ускоряет процесс, поскольку нам приходится рассчитывать меньше вещей.

Функция плотности S(T)

Чтобы начать обретать интуицию в задаче, которую мы должны решить, хорошо начать с создания функции и ее преобразования Фурье, которое при некоторых начальных значениях, S0, T, r, sigma, вычисляет вероятность S_T = X.

Метод расчета функции плотности S(T)

Важно отметить, что интервал [-12 сигма, 12 сигма] покрывает 99. плюс еще 33 9-ки! Кто-то может сказать, что это экстремально, но важно использовать высокую сигму, чтобы функция плотности охватывала глубокие варианты. (Я вернусь к этому позже)

Функция плотности S(T) с S0=100, T=5Y, r=5%, сигма=10% — Источник: Notebook

N Avg Residual
8 1.897794e-02
16 7.159815e-03
32 5.616416e-04
64 2.326956e-06
128 1.384069e-10

Варианты ценообразования с рядами Фурье

Шаг 1: Формула Эйлера для преобразования cos(x)

Мы собираемся начать наше путешествие по выводу функции для оценки опциона с рядами Фурье, используя формулу Эйлера для преобразования cos(x). Это поможет со следующими шагами, которые мы должны сделать.

Шаг 2: Характеристическая функция случайной величины

В теории вероятностей характеристическая функция — это комплекснозначная функция, определенная для случайной величины. Он фиксирует моменты и свойства распределения случайной величины, позволяя анализировать различные статистические свойства с использованием таких методов, как преобразования Фурье. Например, для случайной величины, следующей нормальному распределению или распределению Пуассона:

Характеристическая функция случайной величины

Шаг 3: Аналитическое решение интеграла с использованием характеристической функции

В нашем случае использования мы можем преобразовать функцию плотности, т.е. log(S_T). Продолжая преобразование f(x) из шага 1 An в сочетании с шагом 2 и зная, что f(x) — плотность случайной величины с известной характеристической функцией:

Знание того, что мы можем преобразовать интеграл в характеристическую функцию, позволяет нам явно вычислить интеграл, что я и делал численно с четверкой из сципиона ранее. Математическое ожидание предполагается для достаточно большого интервала [a, b], поскольку f(z) является функцией плотности случайной величины.

Шаг 4: Связываем все это вместе, чтобы явно рассчитать коэффициенты An

Где En — часть интеграла за пределами [a,b], который очень мал, если для массы вероятности f(x) мы выберем a,b, где f(x) близко к нулю (±12sigma, для нормали). Аналогичная концепция может быть применена, если f(x) является выигрышем опциона, а выигрыш очень мал для a,bdef Fourier_cosine_char_function(char_function, x:np.array, a, b, N:int):
A = lambda n : 2/(b-a) * np.real( char_function((n*pi)/(b-a))*exp((-1j*n*pi*a)/(b-a)) )
fx = .5*A(0)
for n in range(1, N+1):
fx += A(n) * cos(n*pi * (x-a)/(b-a))

return fx

S0 = 100
K = 100
r = 0.05
sigma = 0.1
T = 5.0

# For normal distribution
mu = 100
sigma = 25
char_function = lambda u: exp(1j*u*mu -.5*sigma**2*u**2)
f = lambda x: norm.pdf(x, loc=mu, scale=sigma)

a = 100 — 12 * sigma
b = 100 + 12 * sigma

Плотность нормального распределения с использованием явного решения с характеристической функцией — Источник: Notebook

N Avg. Residual Execution Time
8 1.182498e-03 0.000059
16 1.588727e-04 0.000092
32 1.160981e-07 0.000180
64 2.983220e-18 0.000358
128 2.986019e-18 0.000767

Замечания

  • Как видите, у нас все очень хорошо. Время выполнения невелико, и мы получаем очень быструю конвергенцию всего с 32 терминами
  • Использование характеристической функции повышает точность приближения рядов Фурье и упрощает расчеты.
  • Использование преобразования мнимых чисел для cos позволяет преобразовать интеграл в коэффициентах в характеристическую функцию функции плотности.

Ценообразование опционов с использованием рядов Фурье

Шаг 1: Резюме

Для функции f, где она представляет собой функцию плотности случайной величины z, мы имеем косинусное разложение Фурье:

Функция плотности случайной величины z с известной характеристической функцией

Шаг 2: Характерная функция плотности будущих цен на активы

Когда мы оцениваем активы, у нас есть функция плотности логарифма стоимости актива, которая принимает параметры: t, T, S(t), r, sigma и выводит вероятность S(T). В стохастическом исчислении легко вычислить характеристическую функцию случайной величины на основе шага 2 предыдущего раздела. Мы собираемся использовать этот новый инструмент для оценки производной F с определенным выигрышем F(T, S_T):

Формула для расчета цены дериватива с известным выигрышем при T

Шаги, которые я сделал выше:

  • Дисконтируем значение ожидаемого выигрыша дериватива на T. (строка 1)
  • Заменим функцию плотности преобразованием рядов Фурье, которое мы вывели выше. (строка 2)
  • Разделите термины, содержащие z, и упакуйте их в коэффициент gn

Шаг 3: Аналитическое решение g(n)

Поскольку интеграл в gn относится к значениям при T, мы можем аналитически вычислить интеграл, если функция выплат является явной. В этом случае мы попытаемся оценить европейский опцион пут с максимальной выплатой {K-S(T), 0}. Следовательно:

Аналитическое решение для коэффициентов Gn

Чтобы получить показанные здесь решения, вы можете использовать онлайн-калькулятор определенного интеграла, и для g0 вы собираетесь использовать предел gn по мере приближения n к 0, чтобы получить решение, поскольку подстановка n = 0 означает, что g0 не может быть определено.

Шаг 4: Характеристическая функция GBM

Теперь, когда у нас есть наше «оружие», мы пропускаем только пули. В этом случае предполагается, что случайная величина z следует геометрическому броуновскому движению при Q. С характерной функцией «пули»:

Характеристическая функция геометрического броуновского движения

def Fourier_BS_Put(S0, K, T, sigma, r, N):
_sigma, _T = max(.1, sigma), max(1, T)
a = log(S0) + r * _T — 12 * _sigma * np.sqrt(_T) # Upper limit -12 std from the expected value
b = log(S0) + r * _T + 12 * _sigma * np.sqrt(_T) # Lower limit -12 std from the expected value

h = lambda n : (n*pi) / (b-a)
g = lambda n : (exp(a) — (K/h(n))*sin(h(n)*(a — log(K))) — K*cos(h(n)*(a — log(K)))) / (1 + h(n)**2)
g0 = K*(log(K) — a — 1) + exp(a)

m = log(S0) + (r — .5*sigma**2)*(T)
gbm_char= lambda u : exp(1j * u * m -0.5*(u**2 * sigma**2)*T)

F = g0
for n in range(1, N+1):
h_n = h(n)
F += 2*bs_char(h_n) * exp(-1j*a*h_n) * g(n)

F = exp((-r*T))/(b-a) * np.real(F)
return F

Цена опциона по разным ценам исполнения. S0=100, T=5Y, r=5%, сигма=10% — Источник: Notebook

N avg_residuals Ex. Time
8 5.075343e-01 0.000185
16 2.556765e-02 0.000179
32 5.810656e-06 0.000265
64 2.319065e-14 0.000461
128 2.319066e-14 0.000787

Замечания

  • Этот метод оказался очень быстрым и эффективным способом расчета стоимости европейского опциона пут
  • Единственные вещи, которые мы контролируем, чтобы повлиять на ошибку в приближении, — это N, количество членов и интервал [a, b] значений S(T), на которые мы смотрим.

Выбор правильных пределов [a, b]

На приведенном ниже графике показаны 2 очень важные вещи, которые необходимо понимать при использовании рядов Фурье.
1. Что лимиты [a=S_T(-z_score), b=S_T(z_score)] должны покрывать все возможные страйк-цены
2. По мере увеличения сигмы увеличивается и погрешность аппроксимации

Иллюстрация правильного выбора пределов распределения и того, как ошибка в аппроксимации увеличивается с большими сигмами — Источник: Notebook

На первом подграфике вы можете видеть, что если бы у нас была цена исполнения при S(T) = 600, мы не смогли бы правильно ее оценить, поскольку мы моделировали значения только от S(-12) до S(12). Следовательно На втором графике видно, что значение такой страйк-цены было бы очень неправильным!
На третьем графике вы можете увидеть сигмы, охватывающие широкий диапазон страйк-цен. Однако из-за формы выплаты сложнее точно аппроксимировать функцию с тем же количеством слагаемых. Следовательно, по мере увеличения сигмы может потребоваться увеличить количество слагаемых N.

Средняя погрешность первого подучастка составляет: 7.676e-12
Средняя погрешность первого подучастка составляет: 2.830e-02
Средняя погрешность первого подучастка составляет: 1.405e-03
*для N = 128

Погрешность аппроксимации по разному числу слагаемых N

Ошибка аппроксимации аналитического решения европейского опциона пут — Источник: Notebook

Я провел простой эксперимент, чтобы увидеть, как уменьшается погрешность на разное количество слагаемых. Вы можете видеть, что у нас экспоненциальная сходимость, в то время как вычислительные затраты линейны. Мечта каждого ученого-вычислительщика. Это означает, что когда мы заменяем характеристическую функцию моделью Хестона, нам не нужно платить штраф за численное интегрирование по методу Монте-Карло, в котором ошибка уменьшается с квадратным корнем из N.

Анализ погрешностей по всем возможным входным значениям для европейского пут

Для акций стоимостью 100 евро. Я смоделировал условия ценообразования Put, где:
К = [0, 500]
T = [0, 30Y]
r = [-3%, 20%]
сигма = [0,5%, 100%]

Визуализация средней погрешности при различных значениях параметров — Источник: Notebook

Замечания

  • Похоже, что для опционов с экспирацией T = 0 аппроксимация имеет повышенную погрешность, но это связано с известной выплатой, которая представляет собой линейную функцию, которую я показал в начале статьи, которую сложнее аппроксимировать, чем кривую функцию. Для борьбы с этим можно использовать большее количество терминов.
  • Сигма имеет линейную зависимость с погрешностью в аппроксимации, поскольку она изменяет форму функции на что-то, что труднее аппроксимировать. (Как я показал ранее) Опять же, для борьбы с этим можно использовать большее количество терминов.

Общая средняя погрешность: 0,00289 EUR по всем заданным диапазонам значений K, T, r и сигма, которые я определил ранее.

Заключение

Шаг за шагом я объяснил, как перейти от общего преобразования Фурье к пользовательским функциям для расчета выигрыша европейского опциона пут.
Однако, если вы хотите оценить деривативы, где выплата не зависит от фиксированных интервалов, вам придется использовать более сложные численные методы. Например, конечная разность или метод Монте-Карло.

Несмотря на то, что до сих пор это был «бесполезный» инструмент, мы заложили основу для ценообразования с использованием модели Хестона, чего мы не можем сделать аналитически. Надеемся, что вы сможете использовать те же принципы для ценообразования с вариантами

Я надеюсь, что вы нашли эту статью полезной, я надеюсь увидеть вас в части 3, где мы увидим, как использовать модель Хестона с рядами Фурье и как откалибровать ее в цепочке опционов.

Часть 3

Содержание:

  • Резюме: Преобразование Фурье, характеристическая функция
  • Использованный материал: модель Хестона, характеристическая функция модели Хестона, пределы интегрирования (кумулянты)
  • Методика расчета IV, Risk Free Rate Model (Nelson Siegel Svensson)
  • Сбор данных о цепочке опционов из yfinance
  • 3D-графики цепочки опционов: (Цена, IV) против (Время) против (Страйк)
  • Калибровка параметров модели Хестона. мин IV (прогнозируемый) vs IV (наблюдаемый)
  • Анализ ошибок: взвешенное среднеквадратичное значение
  • Поверхность Хестона IV
  • Обсуждение: Ограничения этого метода. Почему опционы Deep In/Out of money трудно правильно оценить. («Ловушка Фурье»)

«Поскольку законы математики относятся к реальности, они не являются достоверными, а поскольку они достоверны, они не относятся к реальности».
Альберт Эйнштейн

Преобразование Фурье, характеристическая функция случайной величины

Давайте быстро вернемся к математическому инструменту, который мы используем для оценки европейского опциона пут. В части 1 я углубился в преобразование Фурье, технику, которая преобразует функцию в серию косинусных волн. В части 2 я обсуждал, как использовать характеристическую функцию распределения случайной величины для точной аппроксимации ее PDF/CDF. Это позволяет нам эффективно оценивать европейский опцион пут.

Резюме: Преобразование Фурье ценообразования опциона с характеристической функцией

Чего не хватает в нашем инструментарии, так это характерной функции распределения и пределов интегрирования (кумулянтов), которые были рассмотрены в части 2. Там я продемонстрировал, как оценить европейский опцион пут с помощью геометрического броуновского движения (GBM), и сравнил его с его аналитическим аналогом. Результаты показали, что этот метод отличается высокой точностью.

Характеристическая функция геометрического броуновского движения

Модель Хестона

Ниже приведены уравнения, которые модель Хестона использует для описания скорости изменения S(t) с помощью модели стохастической волатильности. Я предполагаю, что вы уже знакомы с моделью Хестона. Вместо того, чтобы перефразировать его основы, давайте сосредоточимся на некоторых его ключевых свойствах, которые делают преобразование Фурье (FT) ценным инструментом для этой модели.

Уравнения модели Хестона для моделирования скорости изменения S_t и sigma_t.

Если мы используем метод Эйлера-Монте-Карло (MC) для аппроксимации dS, мы столкнемся с двумя проблемами:

  • Эффективность времени: MC занимает значительно больше времени, чем FT. При калибровке по цепочке опционов, содержащей около 6000 европейских опционов, каждый набор параметров, протестированных в MC, занимал около 10 минут. Напротив, FT заняла менее секунды. Это делает MC вычислительно неосуществимым, когда вам нужно протестировать сотни комбинаций параметров.
  • Отрицательные значения в Vt: Термин ζ×sqrt{dt×Vt}×Z может давать отрицательные значения для Vt. Это проблематично, потому что квадратный корень из отрицательного числа не определен в области действительных чисел. Хотя вы можете изменить модель для использования |Vt| и очень маленький dt, этот обходной путь делает невозможным сравнение результатов MC и FT.

Существуют альтернативные подходы для MC, такие как добавление дополнительных членов Ito или явное вычисление интеграла Vt.

Источник: Милан Мразек и Ян Поспишил [1]

Характерная функция модели Хестона

Я настоятельно рекомендую взглянуть на «Ловушку Маленького Хестона»[2]. Этот источник предлагает подробное объяснение правильного уравнения и динамики характеристической функции модели Хестона. Понимание этой характеристической функции имеет решающее значение для достижения устойчивых приближений, независимо от параметров, с которыми вы работаете.

В этом проекте я опирался на характеристическую функцию, представленную в моем курсе «Стохастическое исчисление». Уравнение выглядит следующим образом:

Характеристическая функция модели Хестона

def heston_char(u):
t0 = 0.0 ; q = 0.0
m = log(S0) + (r — q)*(T-t0)
D = sqrt((rho*zeta*1j*u — kappa)**2 + zeta**2*(1j*u + u**2))
C = (kappa — rho*zeta*1j*u — D) / (kappa — rho*zeta*1j*u + D)
beta = ((kappa — rho*zeta*1j*u — D)*(1-exp(-D*(T-t0)))) / (zeta**2*(1-C*exp(-D*(T-t0))))
alpha = ((kappa*theta)/(zeta**2))*((kappa — rho*zeta*1j*u — D)*(T-t0) — 2*log((1-C*exp(-D*(T-t0))/(1-C))))
return exp(1j*u*m + alpha + beta*v0)

? Внимание: будьте осторожны с α и β значениями в уравнении. Это не то же самое, что пределы интеграции. (Не спрашивайте меня, откуда я знаю.)

Лимиты интеграции, кумулянты

Кумулянты, как определено в Википедии, по сути, являются экспоненциальными центральными моментами распределения. Характеристическая функция и кумулянты неразрывно связаны, что позволяет нам использовать второй кумулянт для расчета дисперсии. Это, в свою очередь, помогает нам определить пределы интеграции [a,b].

В части 2 я обсудил первый и второй кумулянты для GBM и то, как они связаны с ограничениями интеграции. Точно так же модель Хестона имеет кумулянтно-порождающую функцию, из которой выводится c2.

Первый и второй кумулянты GBM и пределы интеграции

Точно так же модель Хестона имеет кумулянтно-порождающую функцию, из которой выводится c2. Источник: «Обобщенный подход с преобразованием Фурье к мерам риска» [3]

Первый и второй кумулянты модели Хестона и пределы интегрирования

Определение этих пределов может быть таким же сложным, как и сама основная модель. Heston PDF/CDF заметно перекошен и демонстрирует более высокий эксцесс по сравнению с GBM. Хотя ранее мы использовали z = 12, более высокое значение zz, такое как z = 24, дает здесь лучшие результаты.

? Ловушка Фурье: значения Z (пределы интегрирования) и N (количество членов) взаимосвязаныУвеличение Z требует более высокого N, но у этого есть предел, который я подробно расскажу позже как «Ловушка Фурье».

Код Python для ценообразования европейского варианта

Ниже приведен полный код Python, который позволяет оценить европейский опцион с использованием преобразованной модели Хестона Фурье. Пара ключевых моментов, на которые следует обратить внимание:

  • Мы оцениваем опцион пут, и формула C=P+S0−K×exp(−r×T) используется для оценки колл-опциона.
  • Эта функция была делом любви, поэтому, если вы найдете ее полезной, подумайте о том, чтобы выпить в мою честь и, возможно, дать ссылку на этот блог.

? Приветствую Numba: Использование JIT-компилятора (JIT), такого как Numba, может ускорить выполнение кода в 164 раза! Мы говорим о 26 мс (Numpy) против 159 мкс (Numba). Итак, давайте поднимем бокал за Numba за то, что он сделал этот процесс молниеносным в Python.import numpy as np
from numpy import sqrt, exp, pi, cos, sin, log, abs
from numba import njit

@njit
def Fourier_Heston_Put(S0, K, T, r,
# Heston Model Paramters
kappa, # Speed of the mean reversion
theta, # Long term mean
rho, # correlation between 2 random variables
zeta, # Volatility of volatility
v0, # Initial volatility
opt_type,
N = 1_012,
z = 24
):

def heston_char(u):
t0 = 0.0 ; q = 0.0
m = log(S0) + (r — q)*(T-t0)
D = sqrt((rho*zeta*1j*u — kappa)**2 + zeta**2*(1j*u + u**2))
C = (kappa — rho*zeta*1j*u — D) / (kappa — rho*zeta*1j*u + D)
beta = ((kappa — rho*zeta*1j*u — D)*(1-exp(-D*(T-t0)))) / (zeta**2*(1-C*exp(-D*(T-t0))))
alpha = ((kappa*theta)/(zeta**2))*((kappa — rho*zeta*1j*u — D)*(T-t0) — 2*log((1-C*exp(-D*(T-t0))/(1-C))))
return exp(1j*u*m + alpha + beta*v0)

# # Parameters for the Function to make sure the approximations are correct.
c1 = log(S0) + r*T — .5*theta*T
c2 = theta/(8*kappa**3)*(-zeta**2*exp(-2*kappa*T) + 4*zeta*exp(-kappa*T)*(zeta-2*kappa*rho)
+ 2*kappa*T*(4*kappa**2 + zeta**2 — 4*kappa*zeta*rho) + zeta*(8*kappa*rho — 3*zeta))
a = c1 — z*sqrt(abs(c2))
b = c1 + z*sqrt(abs(c2))

h = lambda n : (n*pi) / (b-a)
g_n = lambda n : (exp(a) — (K/h(n))*sin(h(n)*(a — log(K))) — K*cos(h(n)*(a — log(K)))) / (1 + h(n)**2)
g0 = K*(log(K) — a — 1) + exp(a)

F = g0
for n in range(1, N+1):
h_n = h(n)
F += 2*heston_char(h_n) * exp(-1j*a*h_n) * g_n(n)

F = exp(-r*T)/(b-a) * np.real(F)
F = F if opt_type == ‘p’ else F + S0 — K*exp(-r*T)
return F if F > 0 else 0

print(‘Speed Analysis of Fourier_Heston_Put: Per Option’)
%timeit Fourier_Heston_Put(S0,K, T, r, kappa, theta, rho, zeta, v0, ‘p’, N = 1_012, z = 24)

«»» Results (Macbook Air M1):
Speed Analysis of Fourier_Heston_Put: Per Option
159 µs ± 8.51 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)
«»»

Подразумеваемая волатильность (IV) и безрисковая модель ставок (модель Нельсона Сигеля Свенссона)

Расчет подразумеваемой волатильности (IV)

Когда вы имеете дело с глубокими опционами «в деньгах» или «вне денег», расчет IV становится числовой игрой. Эта задача может быть столь же ресурсоемкой, как преобразование Фурье (FT). Как правило, используется подход, основанный на методе Ньютона, который, хотя и точен, может занять много времени из-за его итеративного характера.

? Инструмент выбора: По этой причине я выбрал библиотеку py_vollib_vectorized.vectorized_implied_volatility Он использует Numba для более быстрых вычислений и невероятно удобен для пользователя.import py_vollib_vectorized
price = 0.10 ; S = 95 ; K = 100 ; t = .2 ; r = .2 ; flag = ‘c’

def implied_volatility(price, S, K, t, r, flag):
return py_vollib_vectorized.vectorized_implied_volatility(
price, S, K, t, r, flag, q=0.0, on_error=’ignore’, model=’black_scholes_merton’,return_as=’numpy’)

print(‘Speed Analysis of the implied volatility Function: Per Option’)
%timeit implied_volatility(price, S, K, t, r, flag)

«»» Result (Mackbook Air M1):
Speed Analysis of the implied volatility Function: Per Option
49 µs ± 2.79 µs per loop (mean ± std. dev. of 7 runs, 10,000 loops each)
«»»

Кривая безрисковой ставки (RFR)

Для определения RFR я использовал казначейские облигации США, доступ к которым можно получить по этой ссылке.

Метод Нельсона Сигеля Свенссона

Чтобы построить кривую доходности, я использовал метод Нельсона Сигеля Свенссона. Уравнение для этого метода выглядит следующим образом:

Метод Нельсона Сигеля Свенссона

? Реализация: Процесс прост. Во-первых, соберите доходность из Казначейства США. Затем сопоставьте эти выходы с помощью модели Нельсона Сигеля Свенссона с библиотекой Python nelson_siegel_svensso(ссылка).from nelson_siegel_svensson.calibrate import calibrate_nss_ols
# Collected at : 06/01/2023
yield_maturities = np.array([1/12, 2/12, 3/12, 4/12, 6/12, 1, 2, 3, 5, 7, 10, 20, 30])
yields = np.array([5.30,5.39,5.50,5.50,5.44,5.11,4.33,3.98,3.70,3.66,3.61,3.98,3.84])

#NSS model calibrate
curve_fit, status = calibrate_nss_ols(yield_maturities,yields)
#Results: NelsonSiegelSvenssonCurve(beta0=1.2595121714700128, beta1=4.367186541206554, beta2=-1457.4165045532893, beta3=1458.9304400823748, tau1=6.83997027106895, tau2=6.877834685907713)

Результирующая подгонка на основе предоставленных данных.- Источник: Записная книжка

Сбор данных о цепочке опционов из yfinance

Фокус на ^SPX (опционы S&P 500)

Мы будем использовать платформу yfinance, чтобы получить доступ к полной цепочке опционов для ^SPX (опционы S&P 500). Но прежде чем погрузиться в анализ, крайне важно очистить и отформатировать данные.

  1. Получите даты истечения срока действия и последнюю цену закрытия S&P 500. Это закладывает основу для остальной части анализа.

SPX = yf.Ticker(«^SPX»)

# Step 1: get last clossing price
S0 = SPX.history(perdiod=’3mo’)[‘Close’][-1] ; print(f’S0={S0}’)

# Step 2: get Expiration Dates
expiration_dates = SPX.options

2. Цикл по датам истечения срока действия: Перебирайте все даты истечения срока действия, чтобы собрать данные как о путах, так и о коллах, что дает нам полное представление о цепочке опционов.

  • Количество дней, прошедших с момента последней сделки: Рассчитайте количество дней, прошедших с момента последней торговли каждым опционом.
  • Вес опции: Вычислите вес каждого варианта, который необходим для калибровки параметров Хестона.
  • Цена опциона: Предположим, что цена каждого опциона является серединой между ценами покупки и продажи. Обратите внимание, что это может быть неточно для неликвидных вариантов, но они будут иметь меньший вес в нашем анализе из-за их распространения.
  • Срок погашения: Рассчитайте дату погашения для каждого опциона на основе даты истечения срока его действия.

from tqdm import tqdm
from datetime import datetime
import yfinance as yf
from tqdm import tqdm

df = []
for expiration_date in tqdm(expiration_dates, desc=’Collecting Option Data’):
# Collect all Calls
_df = SPX.option_chain(expiration_date).calls
_df[‘DaysSinceLastTraded’] = (pd.to_datetime(_df[‘lastTradeDate’]).dt.tz_localize(None) — pd.to_datetime(datetime.now())).dt.days
_df = _df[_df[‘DaysSinceLastTraded’] > -14]
_df[‘Weight’] = 1 / (_df[‘bid’] — _df[‘ask’])**2
_df[‘price’] = (_df[‘bid’] + _df[‘ask’])/2
_df[‘maturity’] = (pd.to_datetime(expiration_date) — pd.to_datetime(datetime.now())).days / 365
_df[‘Type’] = ‘c’
_df[‘S’] = S0
df.append(_df)

# Collect all Puts
_df = SPX.option_chain(expiration_date).puts
_df[‘DaysSinceLastTraded’] = (pd.to_datetime(_df[‘lastTradeDate’]).dt.tz_localize(None) — pd.to_datetime(datetime.now())).dt.days
_df = _df[_df[‘DaysSinceLastTraded’] > -14]
_df[‘Weight’] = 1 / (_df[‘bid’] — _df[‘ask’])**2
_df[‘price’] = (_df[‘bid’] + _df[‘ask’])/2
_df[‘maturity’] = (pd.to_datetime(expiration_date) — pd.to_datetime(datetime.now())).days / 365
_df[‘Type’] = ‘p’
_df[‘S’] = S0
df.append(_df)

df = pd.concat(df)
df = df[df[‘maturity’] > 0]

3. Отфильтруйте опционы по дате их последней сделки, чтобы свести к минимуму проблемы с подгонкой, связанные с «нерелевантными» опционами. В идеале, наличие точных деталей, таких как время, когда опцион был продан, и позиция S&P 500 в этот момент, помогло бы рассчитать подразумеваемую волатильность (IV). Однако такой уровень детализации, как правило, недоступен. Можно было бы изучить альтернативные методы аппроксимации цены опциона, помимо простого определения средней точки между спросом и предложением.

Распределение дат экспирации и количества доступных опционов по сравнению с DaysSinceLastTraded

4. Примените модель Нельсона-Сигела-Свенссона для оценки безрисковой ставки, соответствующей сроку погашения каждого опциона. Эта модель предлагает надежный способ учета срочной структуры процентных ставок.volSurfaceLong[‘rate’] = volSurfaceLong[‘maturity’].apply(curve_fit) / 100
volSurfaceLong = volSurfaceLong.reset_index()

3D-графики цепочки опционов: цена, IV, время и забастовка

Для европейских опционов, доступных на S&P 500, я использовал Plotly для визуализации цепочки опционов, поскольку он отлично справляется с рендерингом 3D-графиков.

Визуализированные варианты, предлагаемые для SPX. Срок погашения, страйк и цена опциона

Мы будем использовать py_vollib_vectorized для расчета подразумеваемой волатильности (IV) для этих опционов.

Визуализированные варианты, предлагаемые для SPX. Зрелость, страйк и подразумеваемая волатильность

Ключевые замечания:

Из сюжетов выделяются следующие моменты:

  1. Ухмылка или улыбка: Ярко выраженная ухмылка или улыбка очевидны, особенно в IV сюжете.
  2. Симметрия в IV: Для опционов с одинаковым страйком и зрелостью IV коллов и путов удивительно близки.
  3. Низкий IV Близкий к истечением срока действия: Опционы, которые близки к экспирации и находятся в деньгах, демонстрируют самый низкий IV.
  4. Перекошенный IV для коротких экспираций: IV имеет тенденцию быть как более высоким, так и более искаженным для вариантов с более коротким сроком действия.

Рекомендации по моделированию:

Эти наблюдения должны быть эффективно инкапсулированы любой моделью, которую мы решим использовать.

Примечание к IV расчету: Учитывая, что наши методы сбора данных и расчета IV не безупречны, существует вероятность неточностей. Например, если в py_vollib_vectorized введены неправильные входные данные, по умолчанию он равен нулю. Учитывая, что IV особенно чувствителен к опционам, срок действия которых приближается к истечению, крайне важно точно ввести S0 и цену транзакционного опциона.

Калибровка параметров модели Хестона

Теперь, когда у нас есть инструмент для ценообразования европейских опционов на основе конкретных параметров и входных данных, следующим шагом является определение оптимальных параметров модели Хестона, которые наилучшим образом соответствуют рыночным данным. Давайте рассмотрим проблемы и нюансы, связанные с этой калибровкой:

  • Метрика оптимизации: цена опциона или подразумеваемая волатильность (% в годовом исчислении) — критическая. Подразумеваемая волатильность (IV) рассматривает все опционы одинаково, в отличие от цены опциона, которая придает больший вес опционам с глубокими деньгами или вне денег.
  • Метод оптимизацииоптимизатор последовательного программирования методом наименьших квадратов (SLSQP), имитация отжига (SA) и т.д. — высокий. Ваш выбор здесь влияет как на качество посадки, так и на время сходки. Я остановил свой выбор на SLSQP из-за его эффективности и надежности.
  • Функция ошибки: MAE, MSE, MAPE (средняя абсолютная процентная ошибка) — высокая. MSE может показаться заманчивым из-за его чувствительности к выбросам и более быстрой конвергенции, но это не всегда лучший выбор.
  • Вес: 1/(спред²),1/(|спред|), 1/(√спред) — средний. Выбор взвешивания является дискреционным, хотя, по моему опыту, опционы ближе к экспирации и при деньгах, как правило, имеют меньшие спреды.
Уравнение: Минимизируйте взвешенную среднеквадратичную ошибку (MSE) IV между наблюдаемыми и прогнозируемыми значениями.
  • Начальные параметры Хестона и ограничения пространства поиска: Обратитесь к тем, которые используются в проекте — Высокие. Исходные параметры должны быть «разумными», чтобы обеспечить эффективную конвергенцию. Вот почему другие предпочитают такие методы, как SA, но, начиная с сильной начальной точки, SLSQP становится более выгодным.
  • Допуск (tol, ftol): Я установил его на уровне 1×10−51×10−5, что достаточно для наших нужд, учитывая процент внутривенного вливания в годовом исчислении, например, 30,05%. Если улучшение меньше, алгоритм завершает работу — обычно до достижения максимального количества итераций.
  • Обработка ошибок: При работе с опционами, близкими к экспирации, или с высокими ценами исполнения преобразование Фурье может испытывать трудности. Вы можете либо отфильтровать эти параметры, либо установить их IV на ноль. Я выбрал последнее, чтобы поддерживать надежный набор параметров.

Используя как преобразование Фурье (FT), так и Numba, калибровка полной цепочки опционов из 4,200 вызовов и путов занимает менее 6 минут на одноядерном ноутбуке. Результаты? Минимальное взвешенное RMSE(IV) составляет 3,4009, а MAE(IV) — 2,7232 (взвешенное) и 1,3806 (невзвешенное). Это средняя погрешность волатильности в годовом исчислении, составляющая всего 1%.

Пример результирующей подгонкиResults:
v0=0.0228 ; kappa=0.0807 ; theta=0.0363 ; zeta=1.1760 ; rho=-0.3021
>>Zeros : 1

# ——————————— Error for Prices ———————————
MAE of Price : 13.6089
WMAE of Price : 1.2354

MAPE of Price : 54.966 %
WMAPE of Price : 77.121 %
Max APE of Price : 4961.760 %

RMSE of Price : 31.5974
WRMSE of Price : 3.9844

# ——————————— Error for IV ———————————

MAE of IV : 2.7232
WMAE of IV : 1.3806

MAPE of IV : 15.250 %
WMAPE of IV : 6.704 %
Max APE of IV : 230.815 %

RMSE of IV : 4.2373
-> WRMSE of IV: 3.4009

Результирующая посадка

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

Что мы ожидаем от алгоритма оптимизации

  1. Сходитесь к надежным параметрам.
  2. Выполняйте эффективно, используя предыдущие попытки.
  3. Избегайте застревания, исследуя различные области пространства параметров.

Код — реализация

from scipy.optimize import minimize
# Define variables to be used in optimization
_K = volSurface[‘strike’].to_numpy()
_C = volSurface[‘price’].to_numpy()
_type = volSurface[‘Type’].to_numpy(dtype=str)
_T = volSurface[‘maturity’].to_numpy()
_r = volSurface[‘rate’].to_numpy()
_S = np.full_like(_K, S0)
_IV = volSurface[‘IV’].to_numpy()
_Weight = volSurface[‘Weight’].to_numpy()

# Parameters to search for in the simualation: ———————————
params = {
«v0»: {«x0»: 0.002874, «lbub»: [1e-3, 1.2]},
«kappa»: {«x0»: 1.6891, «lbub»: [1e-3, 10]},
«theta»: {«x0»: 0.0190, «lbub»: [1e-3, 1.2]},
«zeta»: {«x0»: 3.7472, «lbub»: [1e-2, 4]},
«rho»: {«x0»: -0.2731, «lbub»: [-1, 1]},
}
x0 = [param[«x0»] for key, param in params.items()]
bnds = [param[«lbub»] for key, param in params.items()]

def SqErr(x):
v0, kappa, theta, zeta, rho = [param for param in x]

# Prices Using Heston Model
Price_Heston = get_resutls_array_Heston(volSurface, v0=v0, kappa=kappa, theta=theta, zeta=zeta, rho=rho, N=1_012, z=24)

# Inverse approximation of implied Volatilities
IV_Heston = implied_volatility(price=Price_Heston, S=_S, K=_K, t=_T, r=_r, flag=_type)

# Safty Feature: Exclude the IV caluclations which are not definable.
diff = IV_Heston — _IV
idx = np.isnan(diff) | np.isinf(diff)
diff[idx] = 0 — _IV[idx]
IV_Heston[idx] = 0

# RMSE
rmse = sqrt(np.mean((diff*100)**2 * _Weight))

zeros = int(np.where(IV_Heston==0, 1,0).sum())
wmae = np.mean(np.abs(diff*100) * _Weight)

print(f»>>v0={v0:7.4}; kappa={kappa:5.4f}; theta={theta:5.4f}; zeta={zeta:5.4f}; rho={rho:6.4f} | WMAE(IV): {wmae:5e} | Nulls: {(idx).sum()}/{idx.shape[0]} | Zeros: {zeros}/{idx.shape[0]} | WRMSE(IV): {rmse:5e}»)
return rmse

# Perform scipy minimize to find the best parameters in the Heston Model which match the observed values
result = minimize(SqErr, x0, tol = 1e-5, method=’SLSQP’, options={‘maxiter’: 80 , ‘ftol’: 1e-5, ‘disp’:True}, bounds=bnds, jac=’3-point’)
v0, kappa, theta, zeta, rho = [param for param in result.x]

«»» Example output:
>>v0=0.01459; kappa=1.0330; theta=0.0268; zeta=3.5362; rho=-0.2757 | WMAE(IV): 1.606438e+00 | Nulls: 0/4283 | Zeros: 6/4283 | WRMSE(IV): 3.789558e+00
>>v0=0.01459; kappa=1.0330; theta=0.0268; zeta=3.5362; rho=-0.2757 | WMAE(IV): 1.606397e+00 | Nulls: 0/4283 | Zeros: 6/4283 | WRMSE(IV): 3.789579e+00
>>v0=0.01459; kappa=1.0330; theta=0.0268; zeta=3.5362; rho=-0.2757 | WMAE(IV): 1.606400e+00 | Nulls: 0/4283 | Zeros: 6/4283 | WRMSE(IV): 3.789582e+00
>>v0=0.01459; kappa=1.0330; theta=0.0268; zeta=3.5362; rho=-0.2757 | WMAE(IV): 1.606403e+00 | Nulls: 0/4283 | Zeros: 6/4283 | WRMSE(IV): 3.789586e+00

Current function value: 3.4043308974510307
Iterations: 41
Function evaluations: 486
Gradient evaluations: 41
Time: 5:52mins

Визуализация посадки

Прогнозируемые цены больше расходятся для опционов с более длительным сроком действия, в первую очередь потому, что алгоритм фокусируется на подгонке опционов ближе к экспирации из-за их меньших спредов. Та же логика применяется при рассмотрении IV, что делает его логически подходящей метрикой для модели Хестона.

Удобство посадки — Цены на опционы
Пригодность — вариант IV

Поверхность Хестона IV

Конечная цель здесь состоит в том, чтобы построить 3D-график поверхности для подразумеваемой волатильности Хестона (IV). Сгенерировав 2D-массив, содержащий данные из 15 000 вариантов, мы можем визуализировать характеристики этой поверхности IV.

В предыдущем обсуждении наблюдаемой поверхности IV были выделены следующие свойства:

  1. Наличие ухмылки или улыбки — эта особенность очевидна в нашей модели Heston.
  2. Самый низкий IV для опционов, которые находятся как в деньгах, так и близки к истечению срока действия — это справедливо и для модели Хестона.
  3. Увеличенное и искаженное внутривенное вливание для краткосрочных экспираций — опять же, модель подтверждает это.
Поверхностный график подразумеваемой волатильности Хестона (IV)

С параметрами:
Хестон: v0=0,001 ; каппа=1,6483 ; тета=0,0261 ; Zeta=4.0000 ; rho=-0.0832
S(0)=4
405.71 долларов США

Устранение численных неустойчивостей

Одна из проблем, которая возникает, заключается в том, что модель не может точно рассчитать значения опционов при определенных параметрах. Это особенно очевидно в следующем графике, в котором основное внимание уделяется опционам со сроком погашения менее одного года.

Численные неустойчивости к преобразованной модели Хестона Фурье

«Ловушка Фурье»: понимание подводных камней ценообразования опционов с FT

Чтобы пролить свет на эти тонкости, я вернулся к предыдущему коду, чтобы построить функцию плотности вероятности (PDF) модели Хестона. Параметры по-прежнему согласуются с графиком поверхности подразумеваемой волатильности (IV), о котором говорилось ранее. PDF-файл соответствует параметру, срок действия которого истекает через один день, т. е. T = 1/365.

График 1: Хестон PDF — Логарифмическая плотность — N=10 000
График 2: Хестон PDF — Логарифмическая плотность — N = 1,000

Очевидно, что график с меньшим количеством слагаемых (N=1000) демонстрирует большую нестабильность, чем его аналог с большим количеством слагаемых (N=10000). Итак, что же дает? Возвращаясь к части 1, мы вспомним, что по мере приближения N к бесконечности ряды Фурье сходятся к исходной функции. Однако предостережение заключается в том, что до этого момента мы работаем с приближениями. Поскольку в нашем приближении используются синусоидальные волны, исходную функцию, близкую к нулю, становится трудно аппроксимировать, особенно для экстремальных z-баллов.def Fourier_cosine_char_function(char_function, x:np.array, a, b, N:int):
A = lambda n : 2/(b-a) * np.real(char_function((n*pi)/(b-a))*exp((-1j*n*pi*a)/(b-a)) )
fx = .5*A(0)
for n in np.arange(1, N+1):
fx += A(n) * cos(n*pi * (x-a)/(b-a))
return fx

Ловушка Фурье

Можно заподозрить, что эти несоответствия как-то связаны с особенностями Heston PDF. Однако это предположение не подтверждается.

Обычный PDF-файл с использованием преобразования Фурье и характеристической функции — N=128
Обычный PDF-файл с использованием преобразования Фурье и характеристической функции — N=100 000

Очевидно, что независимо от количества используемых терминов (N) приближение достигает точки, когда оно перестает сходиться с исходной функцией. Этот подводный камень я назвал «ловушкой Фурье». Я предполагаю, что подобное явление происходит и в ценообразовании опционов пут.

Нормальная CDF с использованием преобразования Фурье и характеристической функции — N=100 000
Нормальный CDF с использованием преобразования Фурье — остатки Abs — N = 100 000

Заключение: ловушка Фурье

Приведенные выше графики объясняют, почему числовая нестабильность возникает в левой части расчета. Синие значения отрицательны, что невозможно в CDF. Точно так же в модели Хестона область, которая отсутствует, связана с тем, что в Хестоне-CDF есть отрицательные значения, что приводит к тому, что цена производной равна 0.

Решение проблемы

  • Увеличение количества сроков (N): Онлайн-источники обычно экономят на терминах при расчете опционов, сосредотачиваясь в основном на «удобных» сроках погашения и опционах, близких к деньгам. Учитывая, что добавление членов в ряды Фурье является вычислительно дешевым, нет ничего плохого в использовании 1,024 членов. Однако после этого момента, выгоды, по-видимому, выходят на плато. Например, увеличение количества членов с 1 024 до 10 000 показывает незначительные улучшения в уменьшении численной нестабильности, но увеличивает время вычислений в 10 раз.
  • Выбор оптимальных кумулянтов (c1, c2): Хотя во многих статьях не учитываются интегральные пределы, статья [3] (2018) является исключением. Многим методологиям, как в статье [4] (2008), кажется, не хватает строгости. Во-первых, кумулянты отличаются от кумулянтов в [3]. Во-вторых, предлагаемый диапазон ±12 * std не подходит для коротких сроков погашения или опционов с глубокими вложениями / вне денег. Наконец, они не тщательно анализируют различные параметры Хестона. После обширного тестирования использование кумулянтов из статьи [3] в сочетании с z = ±24 дало наименьшую среднюю квадратичную ошибку (MSE).
  • Установка соответствующего диапазона (z): Учитывая высокий эксцесс и асимметрию в функции плотности вероятности Хестона (PDF), особенно когда rho значительно меньше нуля, крайне важно использовать достаточно широкий диапазон. Значение z=±24 кажется оптимальным, исходя из самых низких результатов MSE

Давайте углубимся в некоторые показательные сюжеты:

Heston PDF — N=1,024 — T=1 неделя
Heston PDF — N=10,024 — T=1 неделя

Взгляд на ловушку Фурье в модели Хестона

  1. Больший z помогает в ценообразовании опционов с глубокими деньгами / вне их.
  2. Чем выше значение z, тем менее стабильны результаты.
  3. Увеличение z требует более высокого N для стабильных результатов.
  4. Большее z увеличивает значение плотности в хвостах, что, возможно, указывает на большую точность.
  5. Левый хвост имеет тенденцию быть более восприимчивым к разрушению, особенно когда rho отрицательный (например, rho = -0,5).
  6. Численная стабильность изменяется в зависимости от конкретных параметров Хестона:
    — Более высокое значение v0 по отношению к тета(θ) дает более стабильные результаты.
    — Более высокая каппа (κ) улучшает стабильность.
    — Повышенная дзета(ζ) приводит к нестабильности.
    — Rho(ρ) не оказывает заметного влияния на стабильность результата.

Дальнейшая работа

  • Задача 1: Разработать функцию для оптимизации переменных z zz и NN на основе заданных параметров модели Хестона. Цель состоит в том, чтобы свести к минимуму численную нестабильность по различным опционным контрактам.
  • Цель 2: Глубокое погружение в лежащую в основе математическую структуру, чтобы определить потенциальные модификации, которые могли бы решить проблему ловушки Фурье. Непосредственной областью исследований является одновременное использование функций синуса и косинуса, поскольку их взаимодействие может сводить на нет друг друга, потенциально давая более стабильные, близкие к нулю значения.

Источники

  • [1] «Калибровка и моделирование модели Хестона», ссылка
  • [2] «Ловушка Маленького Хестонассылка
  • [3] «Обобщенный подход к мерам риска с преобразованием Фурье», ссылка
  • [4] «НОВЫЙ МЕТОД ЦЕНООБРАЗОВАНИЯ ДЛЯ ЕВРОПЕЙСКИХ ОПЦИОНОВ, ОСНОВАННЫЙ НА РАЗЛОЖЕНИЯХ РЯДОВ ФУРЬЕ-КОСИНУСА», ссылка (2008)

Хестон, Дисперсия-Гамма и Блэк-Шоулз

В этой статье я покажу, как можно оценивать опционы с помощью быстрого преобразования Фурье (БПФ). Концепции и часть кода, которые я здесь представлю, основаны на курсе «Финансовая инженерия и риск-менеджмент» Колумбийского университета. Он предлагается через Coursera, и я очень рекомендую его.

БПФ — это математический метод, который эффективно вычисляет дискретное преобразование Фурье и был адаптирован для ценообразования опционов благодаря своим вычислительным преимуществам. Используя БПФ в оценке опционов, вы можете быстро и точно определить цены различных производных финансовых инструментов.

Во-первых, я покажу общую концепцию и то, как выполнять вычисления. Затем я покажу 3 модели ценообразования и то, как вычислять цены опционов на основе их характерных функций. Мы будем работать с 3 моделями: Black-Scholes, Heston и Variance Gamma.

Метод БПФ может быть применен также к капитализации, минимальному уровню и свопам, которые могут быть выражены в виде простых опционов пут и колл, что позволяет быстро оценить эти инструменты.

Основы

Несмотря на то, что опцион можно оценить путем интеграции, это не является эффективным методом. Лучшим методом является ценообразование на БПФ.

Конечно, мы могли бы написать наши явные функции БПФ для определения цены опций, однако модуль NumPy np.fft.fft намного быстрее, поэтому мы будем использовать его.

Существует ряд правил, которым мы должны следовать, чтобы наши расчеты были правильными.
1) log Strike — довольно просто, вам просто нужно написать, например, Strike = 100 как log(100)
2) Для того, чтобы БПФ было эффективным, вы всегда должны использовать 2^n (векторы, например, 2³ = 8)
3) Логарифмический процесс определения цены акций
4) Определенные параметры: eta, alpha, n и betaS0 = 100
K = 80
k = np.log(K)
r = 0.055
q = 0.03

Ниже приведены примеры параметров# Parameters for FFT
n = 12
N = 2**n
#step-size
eta = 0.25
# damping factor
alpha = 1

lambda_ = (2*np.pi/N)/eta;
beta = np.log(K)

Это основы. Теперь давайте рассмотрим, с какими моделями будет работать и реализовывать ценообразование БПФ.

Блэк-Шоулз

Модель Блэка-Шоулза — это фундаментальная математическая модель, используемая в финансах для расчета теоретической цены опционов европейского типа.

  • Эффективные рынки: Модель предполагает, что финансовые рынки эффективны и на них нет возможностей для арбитража. Он предполагает непрерывную торговлю, отсутствие транзакционных издержек и то, что вся информация отражается в текущей рыночной цене.
  • Логарифмически нормальное распределение: Модель предполагает, что цена базового актива следует логарифмически нормальному распределению. Это предположение в сочетании с непрерывной торговлей приводит к геометрическому броуновскому движению цены акций.
  • Постоянная волатильность: предполагает, что волатильность базового актива постоянна в течение всего срока действия опциона. Это предположение является одним из основных ограничений модели Блэка-Шоулза, поскольку оно противоречит эмпирическим наблюдениям, согласно которым волатильность имеет тенденцию изменяться с течением времени.
  • Безрисковая процентная ставка: Безрисковая процентная ставка считается известной и постоянной на протяжении всего срока действия опциона.

Модель Блэка-Шоулза является основополагающим инструментом в ценообразовании опционов, обеспечивая теоретическую основу для оценки справедливой цены европейских опционов при определенных допущениях. Несмотря на свою простоту и различные допущения, он заложил основу для дальнейшего развития теории финансовых деривативов и опционного ценообразования. Тем не менее, его ограничения в учете изменяющейся волатильности и другой динамики рынка привели к разработке более сложных моделей, таких как модель Хестона и модель Variance Gamma, упомянутые ниже.

Следующие уравнения представляют процесс БС и его характеристическую функцию, необходимую для расчетов БПФ:

Хестон

Модель Хестона — это широко используемая стохастическая модель волатильности в финансах, которая устраняет некоторые ограничения модели Блэка-Шоулза путем введения стохастической волатильности. Модель направлена на то, чтобы отразить улыбку волатильности, наблюдаемую на рынках опционов, и тенденцию волатильности колебаться с течением времени.

  • Стохастическая волатильность: В отличие от модели Блэка-Шоулза, где волатильность предполагается постоянной, модель Хестона включает в волатильность стохастичность. Он предполагает, что сама дисперсия следует стохастическому процессу, обычно моделируемому как процесс диффузии квадратного корня с возвратом к среднему значению (процесс CIR).
  • Возврат к среднему значению: Стохастическая модель волатильности Хестона включает в себя поведение возврата к среднему значению, при котором волатильность имеет тенденцию двигаться к долгосрочному среднему уровню. Это свойство обращения к среднему значению помогает уловить наблюдаемое поведение волатильности на финансовых рынках.
  • Улыбка волатильности: Одной из основных мотиваций для модели Хестона является ее способность создавать улыбку волатильности, эмпирическое явление, при котором подразумеваемая волатильность изменяется в зависимости от цены исполнения и срока погашения на рынках опционов.
  • Корреляция между ценой актива и волатильностью: Модель Хестона допускает параметр корреляции (обычно обозначаемый ρ) между процессом цены актива и процессом волатильности. Эта корреляция отражает взаимосвязь между ценами на активы и соответствующей волатильностью. Если корреляция равна нулю, это означает независимость между двумя процессами.

Введение стохастической волатильности в модели Хестона устраняет некоторые недостатки модели Блэка-Шоулза, обеспечивая более динамичное представление поведения рынка. Его способность улавливать улыбку волатильности и учитывать тенденцию волатильности колебаться с течением времени делает его популярным выбором в приложениях для ценообразования опционов и управления рисками на финансовых рынках.

Дисперсия-Гамма

Модель Variance Gamma (VG) является расширением модели Блэка-Шоулза, которая включает в себя как скачки цен на активы, так и стохастическую волатильность

  • Включение скачков: Модель VG вводит компоненты скачков в ценовую динамику, позволяя более реалистично вести себя на финансовых рынках. Эти скачки помогают фиксировать внезапные и крупные движения, которые время от времени происходят в ценах на активы.
  • Гибкость в асимметрии и эксцессе: VG обеспечивает большую гибкость в моделировании асимметрии (асимметрии) и эксцесса (жирного хвоста) по сравнению с более простыми моделями, такими как Блэк-Шоулз.
  • Стохастическая волатильность: Подобно модели Хестона, VG включает в себя стохастическую волатильность, позволяя самой волатильности стохастически изменяться с течением времени. Эта функция обеспечивает лучшее представление об эмпирическом поведении волатильности на финансовых рынках.

Если параметры асимметрии и эксцесса равны нулю, модель VG сходится к модели Блэка-Шоулза. Кроме того, модель VG сводится к модели скачкообразной диффузии Мертона, когда член стохастической волатильности равен нулю.

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

Характеристическая функция VG логарифма цены акции:

Реализация кода

После рассмотрения моделей мы будем работать с реализацией кода. Первая функция «compute_characteristic_function» как раз и является прямым выполнением приведенных выше формул.

Функция принимает входные параметры, специфичные для каждой модели, и генерирует соответствующую характеристическую функцию. Переменная phi содержит результирующую характеристическую функцию, которая затем возвращается.def compute_characteristic_function(u, params, S0, r, q, T, model):
if model == ‘BS’:
volatility = params[0]
drift = np.log(S0) + (r — q — volatility ** 2 / 2) * T
diffusion = volatility * np.sqrt(T)
phi = np.exp(1j * drift * u — (diffusion * u) ** 2 / 2)
elif model == ‘Heston’:
kappa = params[0]
theta = params[1]
sigma = params[2]
rho = params[3]
v0 = params[4]

temp1 = (kappa — 1j * rho * sigma * u)
g = np.sqrt((sigma ** 2) * (u ** 2 + 1j * u) + temp1 ** 2)
power_term = 2 * kappa * theta / (sigma ** 2)
numerator = (kappa * theta * T * temp1) / (sigma ** 2) + 1j * u * T * r + 1j * u * np.log(S0)
log_denominator = power_term * np.log(np.cosh(g * T / 2) + (temp1 / g) * np.sinh(g * T / 2))
temp2 = ((u * u + 1j * u) * v0) / (g / np.tanh(g * T / 2) + temp1)
log_phi = numerator — log_denominator — temp2
phi = np.exp(log_phi)
elif model == ‘VG’:
sigma = params[0]
nu = params[1]
theta = params[2]

if nu == 0:
mu = np.log(S0) + (r — q — theta — 0.5 * sigma ** 2) * T
phi = np.exp(1j * u * mu) * np.exp((1j * theta * u — 0.5 * sigma ** 2 * u ** 2) * T)
else:
mu = np.log(S0) + (r — q + np.log(1 — theta * nu — 0.5 * sigma ** 2 * nu) / nu) * T
phi = np.exp(1j * u * mu) * ((1 — 1j * nu * theta * u + 0.5 * nu * sigma ** 2 * u ** 2) ** (-T / nu))
return phi

Функция «genericFFT» по сути использует БПФ для эффективного вычисления цен опционов путем преобразования значений характеристических функций и страйк-цен в подходящее пространство для расчета цен, что обеспечивает более быструю и точную оценку опционов.def calculate_fft_values(params, S0, K, r, q, T, alpha, eta, n, model):
«»»
S0: Initial asset price.
K: Strike price of the option.
r: Risk-free interest rate.
q: Continuous dividend yield.
T: Time to maturity.
alpha: Constant used in the Carr-Madan method for damping.
eta: Step size for numerical integration in FFT.
n: Exponential factor to determine the number of steps in the FFT.
«»»
N = 2 ** n
delta = (2 * np.pi / N) / eta
beta = np.log(K)
km_values = np.zeros(N)
x_values = np.zeros(N)
# discount factor
discount_factor = np.exp(-r * T)

nuJ = np.arange(N) * eta
psi_nuJ = generic_characteristic_function(nuJ — (alpha + 1) * 1j, params, S0, r, q, T, model) / (
(alpha + 1j * nuJ) * (alpha + 1 + 1j * nuJ))

km_values = beta + delta * np.arange(N)
w_values = eta * np.ones(N)
w_values[0] = eta / 2

x_values = np.exp(-1j * beta * nuJ) * discount_factor * psi_nuJ * w_values
y_values = np.fft.fft(x_values)

cT_km_values = np.zeros(N)
multiplier = np.exp(-alpha * km_values) / np.pi

cT_km_values = multiplier * np.real(y_values)

return km_values, cT_km_values

А теперь код для ценообразования путов. Он инициализирует матричный put_matrix для хранения результатов. Каждая строка представляет собой уникальную комбинацию параметров, а столбцы представляют eta, n, alpha и цену пут.def calculate_put_prices(params, S0, K, r, q, T, model, alpha_values, eta_values, n_values):
num_prices = len(eta_values) * len(n_values) * len(alpha_values)
# Columns correspond to eta, n, alpha, and put price
put_matrix = np.zeros([num_prices, 4])
i = 0
for eta in eta_values:
for n in n_values:
for alpha in alpha_values:
# Pricing via FFT
put = 0
k_values, option_values = genericFFT(params, S0, K, r, q, T, alpha, eta, n, model)
put = option_values[0] # Considering only one strike
put_matrix[i] = np.array([eta, n, alpha, put])
i += 1
return put_matrix

Теперь давайте воспользуемся кодом и оценим некоторые опции! Начнем с Блэка Шоулзаmod = ‘BS’
sig = 0.3
params = [sig]def print_put_prices(params, S0, K, r, q, T, model, alpha_values, eta_values, n_values):
put_matrix = calculate_put_prices(params, S0, K, r, q, T, model, alpha_values, eta_values, n_values)
num_prices = put_matrix.shape[0]

print(‘Model = ‘ + model)
print(‘eta\tN\talpha\tput’)

for i in range(num_prices):
print(‘%.2f\t2^%i\t%.2f\t%.4f’ % (put_matrix[i, 0], put_matrix[i, 1], put_matrix[i, 2], put_matrix[i, 3]))

Вывод:
— — — — — — — — — — —
Модель = BS
eta N alpha put
0,10 2⁶ -1,01 89,7431
0,10 2⁶ -1,25 2,7396
0,10 2⁶ -1,50 2,7569
0,10 2⁶ -1,75 2,7701
0,10 2⁶ -2,00 2,7789
0,10 2⁶ -5,00 2,6727
0,10 2¹⁰ -1,01 89,7316
0,10 2¹⁰ -1,25 2,7080
0.10 2¹⁰ -1.50 2.7080
0,10 2¹⁰ -1,75 2,7080
0,10 2¹⁰ -2,00 2,7080
0,10 2¹⁰ -5,00 2,7080
0,25 2⁶ -1,01 269,0367
0,25 2⁶ -1,25 2,8504
0,25 2⁶ -1,50 2,7083
0,25 2⁶ -1,75 2,7080
0,25 2⁶ -2,00 2,7080
0,25 2⁶ -5,00 2,7080
0,25 2¹⁰ -1,01 269,0367
0,25 2¹⁰ -1,25 2,8504
0,25 2¹⁰ -1,50 2,7083
0,25 2¹⁰ -1,75 2,7080
0,25 2¹⁰ -2,00 2,7080
0,25 2¹⁰ -5,00 2,7080

Хестон:mod = ‘Heston’
kappa = 2.
theta = 0.05
lda = 0.3
rho = -0.7v0 = 0.04
params = [kappa, theta, lda, rho, v0]
heston_matrix = price_all_puts(params, S0, K, r, q, T, mod, alpha_vec, eta_vec, n_vec)

print(‘Model = ‘ + mod)
print(‘eta\tN\talpha\tput’)
for i in range(num_prices):
print(‘%.2f\t2^%i\t%.2f\t%.4f’ % (heston_matrix[i,0], heston_matrix[i,1], heston_matrix[i,2], heston_matrix[i,3]))

Вывод:
— — — — — — — — — — —
Модель = Heston
eta N alpha put
0,10 2⁶ -1,01 88,0744
0,10 2⁶ -1,25 1,1102
0,10 2⁶ -1,50 1,1666
0,10 2⁶ -1,75 1,2167
0,10 2⁶ -2,00 1,2605
0,10 2⁶ -5,00 1,3979
0,10 2¹⁰ -1,01 88,3648
0,10 2¹⁰ -1,25 1,3412
0,10 2¹⁰ -1,50 1,3412
0,10 2¹⁰ -1,75 1,3412
0,10 2¹⁰ -2,00 1,3412
0,10 2¹⁰ -5,00 1,3412
0,25 2⁶ -1,01 267,6738
0,25 2⁶ -1,25 1,4873
0,25 2⁶ -1,50 1,3450
0,25 2⁶ -1,75 1,3445
0,25 2⁶ -2,00 1,3442
0,25 2⁶ -5,00 1,3415
0,25 2¹⁰ -1,01 267,6698
0,25 2¹⁰ -1,25 1,4835
0,25 2¹⁰ -1,50 1,3414
0,25 2¹⁰ -1,75 1,3412
0,25 2¹⁰ -2,00 1,3412
0,25 2¹⁰ -5,00 1,3412

и дисперсия-гамма:# model-specific parameters
mod = ‘VG’
sigma = 0.3
nu = 0.5
theta = -0.4
params = [sigma, nu, theta]vg_matrix = price_all_puts(params, S0, K, r, q, T, mod, alpha_vec, eta_vec, n_vec)
print(‘Model = ‘ + mod)
print(‘eta\tN\talpha\tput’)
for i in range(num_prices):
print(‘%.2f\t2^%i\t%.2f\t%.4f’ % (vg_matrix[i,0], vg_matrix[i,1], vg_matrix[i,2], vg_matrix[i,3]))

Вывод:
— — — — — — — — — — —
Модель = VG
eta N alpha put
0,10 2⁶ -1,01 92,2103
0,10 2⁶ -1,25 5,2047
0,10 2⁶ -1,50 5,2230
0,10 2⁶ -1,75 5,2401
0,10 2⁶ -2,00 5,2556
0,10 2⁶ -5,00 0,0030
0,10 2¹⁰ -1,01 92,3373
0,10 2¹⁰ -1,25 5,3137
0,10 2¹⁰ -1,50 5,3137
0,10 2¹⁰ -1,75 5,3137
0,10 2¹⁰ -2,00 5,3137
0,10 2¹⁰ -5,00 0,0000
0,25 2⁶ -1,01 271,6399
0,25 2⁶ -1,25 5,4539
0,25 2⁶ -1,50 5,3121
0,25 2⁶ -1,75 5,3122
0,25 2⁶ -2,00 5,3124
0,25 2⁶ -5,00 0,0019
0,25 2¹⁰ -1,01 271,6423
0,25 2¹⁰ -1,25 5,4560
0,25 2¹⁰ -1,50 5,3139
0,25 2¹⁰ -1,75 5,3137
0,25 2¹⁰ -2,00 5,3137
0,25 2¹⁰ -5,00 0,0020

Интерпретация результатов

  • eta: представляет собой размер шага для численного интегрирования в методе БПФ. Меньшие значения обычно дают более точные результаты, но требуют больше вычислительных ресурсов.
  • n: Указывает экспоненциальный коэффициент, используемый при определении числа шагов в БПФ. Большее n означает более точные вычисления, но повышенную вычислительную сложность.
  • alpha: Это константа, используемая в методе Карра-Мадана для демпфирования. Это влияет на сглаживание кривой цены опциона, сгенерированной с помощью БПФ.
  • put: Обозначает вычисленную цену пут-опциона.

Наблюдая за строками для различных комбинаций eta, n и alpha, вы можете увидеть, как изменение этих параметров влияет на результирующие цены опционов пут.

Результаты не следует интерпретировать вне контекста! Целью сравнения значений является сравнение их с аналитическими решениями и определение правильного выбора параметров для дальнейшего ценообразования.

Заключение

Таким образом, представленный анализ предлагает всестороннюю перспективу влияния различных параметров на цены пут-опционов в рамках модели VG. Этот анализ подчеркивает чувствительность модели к изменяющимся входным параметрам, что имеет решающее значение для обеспечения точности и стабильности ценообразования.

Изучение ценообразования опционов на основе БПФ дает базовое понимание, подчеркивая ключевую роль чувствительности к параметрам. Эти знания помогут специалистам, занимающимся количественными финансами, эффективно оценивать деривативы. Тем не менее, это также подчеркивает необходимость точного выбора параметров для обеспечения точности ценообразования в практике финансовой оценки.

Ссылки:

  • Введение в математику производных финансовых инструментов — Али Хирса, Салих А. Нефтчи
  • Специализация «Финансовая инженерия и риск-менеджмент»; Вычислительные методы в ценообразовании и калибровке моделей

Прогнозирование с помощью рядов Фурье

Основная цель — показать, как ряды Фурье могут быть подогнаны для создания прогнозов временных рядов для высокосезонных данных, как это делает очень популярный прогнозист Prophet.

Prophet, разработанный инженерами Meta, является «аналитиком в цикле» для данных временных рядов. Он был разработан с учетом бизнес-сценария использования Facebook, где их данные, как правило, демонстрируют ярко выраженные сезонные свойства, в отличие от данных авторегрессии — подумайте о том, как веб-трафик меняется по часам в день, день в неделю или неделю в году, в отличие от цен на акции, где сезонность очень мала, а будущие цены могут быть сильно коррелированы с ценой в предыдущую единицу времени.

Также очень важна роль «аналитика в петле». Создатели ожидают, что пользователи будут обладать глубокими знаниями предметной области (в отличие от опыта прогнозирования), чтобы настраивать параметры прогнозиста. Интересное исследование, признанное инженерами Meta в своем блоге, показало, что другие прогнозисты (такие как AutoArima) могут быть значительно лучше по сравнению с параметрами Prophet по умолчанию (потенциально), развенчивая идею о том, что Prophet является лучшим «готовым» инструментом для таймсерий.

Помня об этом, нам казалось разумным, что мы могли бы воссоздать (неидентичную) версию Prophet и попутно изучить моделирование временных рядов.

Ряды Фурье

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

При указании типов сезонности (суточный, недельный, годовой и т.д.) в прогнозе Prophet пользователь также может определить порядок Фурье для каждого типа сезонности, это относится к количеству рядов Фурье, которые суммируются для каждого компонента сезонности.

Пример со страницы в Википедии. По мере увеличения количества терминов посадка улучшается.

Глядя на код Профета, мы видим, что для каждого порядка Фурье создаются синусоида и косинус. Для первого триместра частота этих волн равна одной за указанные периоды для данной сезонности (например, 1/7 для недельных). И по мере увеличения порядка Фурье частота увеличивается пропорционально (например, 1/7, 2/7, 3/7).@staticmethod
def fourier_series(
dates: pd.Series,
period: Union[int, float],
series_order: int,
) -> NDArray[np.float_]:
«»»Provides Fourier series components with the specified frequency
and order.

Parameters
———-
dates: pd.Series containing timestamps.
period: Number of days of the period.
series_order: Number of components.

Returns
——-
Matrix with seasonality features.
«»»
if not (series_order >= 1):
raise ValueError(«series_order must be >= 1»)

# convert to days since epoch
t = dates.to_numpy(dtype=np.int64) // NANOSECONDS_TO_SECONDS / (3600 * 24.)

x_T = t * np.pi * 2
fourier_components = np.empty((dates.shape[0], 2 * series_order))
for i in range(series_order):
c = x_T * (i + 1) / period
fourier_components[:, 2 * i] = np.sin(c)
fourier_components[:, (2 * i) + 1] = np.cos(c)
return fourier_components

Подгонка синусоидальных волн

Prophet использует Stan, вероятностный фреймворк программирования для статистического вывода. Тем не менее, для этой реализации я хотел попробовать градиентный спуск.

Для любого алгоритма градиентного спуска необходима функция стоимости, чтобы обеспечить обратную связь для обновлений на каждой итерации. Среднеквадратичная ошибка является одной из наиболее часто используемых стоимостных функций. Производная первого порядка вычисляется путем умножения на член экспоненты (2) и уменьшения экспоненты на 1.

ошибка = y_pred — ymean_squared_error = (1 / n_samples) * (погрешность) ** 2mse_derivative = (2 / n_samples) * (ошибка)

Это обеспечивает градиент затрат в любой заданной точке и, следовательно, в какую сторону двигаться для следующей итерации и насколько. Вот пример аппроксимации весов линейной регрессии — член смещения можно рассматривать как часть регрессоров X:
X = …
y = …
learning_rate = .001
n_iterations = 10_000

n_samples, n_features = X.shape
params = np.zeros(n_features, dtype=np.float64)

for _ in range(n_iterations):
# create new set of preds
y_pred = X @ params
# calculate error term
error = y_pred — y
# calculate cost for each weight
dw = X.T @ ((2 / n_samples) * (error))
# update weights
params -= learning_rate * dw

В отличие от линейной регрессии, где для каждого признака вычисляется один параметр, для синусоидальных волн необходимо подобрать три параметра (амплитуда, фаза и частота). Но в любом случае принципы остаются прежними — для каждой итерации каждый параметр сначала увеличивается или уменьшается на небольшую величину по отдельности (чтобы узнать направление и величину изменения), а затем все они обновляются вместе для следующей итерации.

Синусоида может быть выражена следующим образом:

амплитуда * sin(фаза + частота * x)

Для каждой амплитуды, фазы и частоты нам нужно знать производную синусоиды по отношению к изменению этого параметра. Вы можете сделать это, либо вспомнив тригонометрию, либо, как в моем случае, обратившись к Wolfram Alpha:

Амплитуда -> d/da(a sin(p + f x)) = sin(f x + p)Фаза -> d/dp(a sin(p + f x)) = a cos(f x + p)Частота -> d/df(A sin(p + f x)) = A x cos(f x + p)

Перед подгонкой необходимо инициировать параметры волны в соответствии с количеством присутствующих сезонностей и количеством заданных ордеров Фурье. Как и в случае с Пророком, для каждого порядка Фурье инициируется волна и волна смещения — помните, что косинусоидальная волна — это синусоида, смещенная на половину цикла.seasonality_terms = {
7.: weekly_seasonality_terms,
30.43: monthly_seasonality_terms,
91.31: quarterly_seasonality_terms,
365.25: yearly_seasonality_terms
}

# no. wave -> total no. fourier terms * 2 (sine + cosine)
n_waves = 2 * sum(seasonality_terms.values())
amplitudes = np.ones(n_waves, dtype=np.float64)
phases = np.zeros(n_waves, dtype=np.float64)
frequencies = np.zeros(n_waves, dtype=np.float64)
count = 0
for periods, terms in seasonality_terms.items():
for j in range(terms):
f = (j + 1) / periods
frequencies[count] = f
frequencies[count + 1] = f
phases[count] = 0 # sine
phases[count + 1] = np.pi # cosine
count += 2

Эти параметры могут быть подогнаны градиентным спуском:y = …
X = …
t = np.arrange(y.size) * 2 * np.pi
n_samples = y.size

learning_rate = .001
n_iterations = 10_000

for _ in range(n_iterations):

# create predictions
y_pred = np.zeros_like(y, dtype=np.float64)
for w in range(n_waves):
y_pred += amplitudes[w] * np.sin(phases[w] + frequencies[w] * t)

# calculate error term
error = y_pred — y

# calculate cost for each weight
cost_derivative = (2 / n_samples) * (error)

# update amp and phase for each wave
for w in range(n_waves):
da = 1 * np.sin(phases[w] + frequencies[w] * t)
dp = amplitudes[w] * np.cos(phases[w] + frequencies[w] * t)

amplitudes[w] -= da.sum() * learning_rate
phases[w] -= dp.sum() * learning_rate

Вы можете заметить, что в этом примере обновляются только фазы и амплитуды, в то время как частоты выводятся непосредственно из заданного количества периодов. Это ловушка, в которую я попал, оказывается, градиентный спуск будет сходиться только для выпуклых функций, а вогнутые функции могут не приводить к стабильным результатам. Производная синусоиды по частоте является вогнутой функцией, и поэтому ее трудно совместить с этим процессом. Вероятно, именно поэтому многие алгоритмы, такие как Prophet, предопределяют периоды для заданных сезонностей (или пользователь может указать свои собственные).

Процесс аппроксимации рядов Фурье также можно комбинировать с аппроксимацией регрессоров в предыдущем примере, а также с трендом и смещением.

Пример из реальной жизни

В этом примере будет использоваться общее количество пассажиров поездов и автобусов от Транспортного управления Чигаго. Набор данных ежедневный и охватывает период с января 2001 года по октябрь 2018 года. Этот набор данных был выбран потому, что кажется разумным ожидать еженедельной и годовой сезонности.

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

Учитывая, что праздники могут влиять на пассажиропоток, праздники Python также реализованы для поиска государственных праздников США.

Полный код FourierForecaster можно найти в этом репозитории GitHub. Это позволяет подгонять тренд и смещение, а также любые указанные регрессоры (например, праздники). Numba также реализован для повышения производительности.

Компоненты сюжета

После подгонки можно вызвать .plot_components() для визуализации смещения, тренда, сезонности, регрессоров (праздников) и шума:

Оценка

Модель обучается на данных за два года и строятся прогнозы на третий год, затем на следующей итерации обучение и прогнозирование переносятся на один год вперед. Это повторяется, чтобы прогнозы составлялись на каждый год с 2003 по 2018 год. Например:

Итерация 1: Обучение 2001 и 2002 гг. Прогнозы: 2003 год.

Итерация 2: Обучение 2002 и 2003 гг. Прогнозы: 2004 год.

Итерация 3: Обучение 2003 и 2004 гг. Прогнозы: 2005 год.

Итерация 16: Обучение 2016 и 2017 гг. Прогнозы: 2018 год.

Метрикой оценки будет MAPE (не обязательно лучшая, но легко интерпретируемая).from fourier_forecast.fourier_forecast import FourierForecast
from datetime import date
import pandas as pd
import numpy as np

df = …
hols = …

train_years = 2
first_year, last_year = 2003, 2018
ff_results = np.zeros(last_year — first_year + 1, dtype=np.float64)
for i, test_year in enumerate(range(first_year, last_year + 1)):

train_dates = (date(test_year — train_years, 1, 1), date(test_year — 1, 12, 31))
test_dates = (date(test_year, 1, 1), date(test_year, 12, 31))

train_index = df[(df[‘ds’] >= train_dates[0]) & (df[‘ds’] <= train_dates[1])].index
test_index = df[(df[‘ds’] >= test_dates[0]) & (df[‘ds’] <= test_dates[1])].index

train_df = df[df.index.isin(train_index)]
test_df = df[df.index.isin(test_index)]
train_hols = hols[train_index]
test_hols = hols[test_index]

ff = FourierForecast()
y = np.log(train_df[‘y’].values)
ff.fit(train_df[‘ds’].values, y, regressors=train_hols)

preds = ff.predict(ds=test_df[‘ds’].values, regressors=test_hols)
preds = np.exp(preds)

mape = np.abs(preds / test_df[‘y’].values — 1).mean()
ff_results[i] = mape

Порядок Фурье по умолчанию для недельной и годовой сезонности равен 3 и 10 соответственно по той причине, что они также являются значениями по умолчанию, используемыми Prophet, что делает сравнение более похожим.

Сравнение с Prophet

Prophet также включает в себя множество параметров, таких как точки изменений, гибкие тренды, мультипликативная сезонность, предыдущие масштабы и т. д. Многие из них могут быть выведены автоматически или настроены аналитиком.

Поскольку Prophet может обрабатывать мультипликативную сезонность внутри компании, логарифмическое преобразование данных не требуется. Праздники можно добавить так же легко. Кроме этих двух, никакие другие параметры настраиваться не будут.from prophet import Prophet
from datetime import date
import pandas as pd

df = …

train_years = 2
first_year, last_year = 2003, 2018
prophet_results = np.zeros(last_year — first_year + 1, dtype=np.float64)
for i, test_year in enumerate(range(first_year, last_year + 1)):

train_dates = (date(test_year — train_years, 1, 1), date(test_year — 1, 12, 31))
test_dates = (date(test_year, 1, 1), date(test_year, 12, 31))

train_index = df[(df[‘ds’] >= train_dates[0]) & (df[‘ds’] <= train_dates[1])].index
test_index = df[(df[‘ds’] >= test_dates[0]) & (df[‘ds’] <= test_dates[1])].index

train_df = df[df.index.isin(train_index)]
test_df = df[df.index.isin(test_index)]

m = Prophet(seasonality_mode=’multiplicative’)
m.add_country_holidays(country_name=’US’)
m.fit(train_df)

future_df = m.make_future_dataframe(periods=test_index.size).tail(test_index.size)
pred_df = m.predict(future_df)

mape = np.abs(preds / test_df[‘y’].values — 1).mean()
prophet_results[i] = mape

Сравнение с AutoArima

Учитывая ссылку в начале статьи, для полноты картины можно взглянуть на реализацию AutoArima в StatsForecast. Похоже, что есть возможность установить только одну сезонную длину, поэтому она установлена на 7, учитывая сильные недельные циклы. Остальные параметры остаются нетронутыми, учитывая заявления о том, что он является лучшим «нестандартным» прогнозистом.from statsforecast.models import AutoARIMA
from datetime import date
import pandas as pd
import numpy as np

df = …
hols = …

train_years = 2
first_year, last_year = 2003, 2018
aa_results = np.zeros(last_year — first_year + 1, dtype=np.float64)
for i, test_year in enumerate(range(first_year, last_year + 1)):

train_dates = (date(test_year — train_years, 1, 1), date(test_year — 1, 12, 31))
test_dates = (date(test_year, 1, 1), date(test_year, 12, 31))

train_index = df[(df[‘ds’] >= train_dates[0]) & (df[‘ds’] <= train_dates[1])].index
test_index = df[(df[‘ds’] >= test_dates[0]) & (df[‘ds’] <= test_dates[1])].index

train_df = df[df.index.isin(train_index)]
test_df = df[df.index.isin(test_index)]
train_hols = hols[train_index]
test_hols = hols[test_index]

aa = AutoARIMA(season_length=7)
aa.fit(train_df[‘y’].values, train_hols.astype(np.float64))
preds = aa.predict((test_dates[1] — test_dates[0]).days + 1, test_hols.astype(np.float64))[‘mean’]

mape = np.abs(preds / test_df[‘y’].values — 1).mean()
aa_results[i] = mape

Результаты

Ни один набор данных не может дать полную картину производительности и возможностей прогнозистов, но для наглядных целей результаты довольно интересны:FourierForecast scores: [0.06757175 0.06847599 0.07174342 0.06946523 0.07235714 0.08831864
0.09115687 0.06550169 0.08012167 0.06707774 0.08542078 0.07360266
0.07506197 0.07625026 0.07252445 0.06323971]
Prophet scores: [0.06361311 0.08916864 0.07080911 0.0656165 0.0788751 0.11362799
0.06367285 0.06227564 0.13398795 0.07868203 0.06452716 0.07208917
0.08931415 0.08308844 0.07541953 0.05568152]
AutoArima: [0.07567813 0.07967367 0.08599838 0.0809436 0.08907974 0.1112463
0.08570553 0.09364843 0.13209986 0.08488383 0.08312174 0.08522938
0.14317758 0.1041176 0.08812333 0.07335536]

FourierForecast average MAPE: 7.4%
Prophet average MAPE: 7.9%
AutoArima average MAPE: 9.4%

В этом конкретном наборе данных и оценочной метрике FourierForecaster показал лучшие результаты, чем Prophet и AutoArima. Тем не менее, возможно, у Prophet есть больше возможностей для улучшения, если параметры прогнозиста будут настроены дальше. И вполне вероятно, что в сценариях, где присутствует более сильная авторегрессионная составляющая, AutoArima может затмить обе модели в других ситуациях.

Ресурсы:

Преобразование Фурье в финансах

В динамичном мире финансов инновационные математические инструменты, такие как преобразование Фурье, все чаще используются для улучшения торговых стратегий. Этот трансформационный подход предлагает уникальную перспективу, особенно в хеджировании и арбитражной торговле.

Фото ThisisEngineering RAEng с сайта Unsplash

Понимание преобразования Фурье в финансах

Преобразование Фурье — это математическое понятие, которое разлагает функцию (в данном случае финансовые данные) на составляющие ее частоты. В финансах это выражается в анализе рыночных данных для выявления основных закономерностей и циклов.

Применение в стратегиях хеджирования

Хеджирование, стратегия управления рисками, используемая для компенсации потенциальных убытков, может значительно выиграть от преобразования Фурье.

Шаг 1: Сбор и подготовка данных

Соберите исторические данные об активе и связанных с ним финансовых показателях.

Шаг 2: Примените преобразование Фурье

Используя библиотеки Python, такие как NumPy или SciPy, примените преобразование Фурье к данным.import numpy as np

# Sample data
prices = np.array([data_points]) # Replace with actual data
fft_result = np.fft.fft(prices)

Шаг 3: Частотный анализ

Изучите выходные данные, чтобы определить доминирующие частоты, связанные с волатильностью или тенденциями рынка.

Шаг 4: Разработка стратегии

Разработайте стратегию хеджирования на основе выявленных закономерностей для снижения риска.

Преобразование Фурье в арбитражных возможностях

Арбитраж включает в себя использование разницы в ценах на различных рынках или в разных формах. Преобразование Фурье может выявлять возможности неправильной оценки, анализируя частотную область ценовых рядов.

Определение возможностей для арбитража

Проанализируйте частотные составляющие двух похожих активов. Значительные отклонения в их частотных паттернах могут указывать на арбитражные возможности.

Пример кода Python

import numpy as np

# Data of two similar assets
prices_asset1 = np.array([data_points_asset1])
prices_asset2 = np.array([data_points_asset2])

# Fourier Transform
fft_asset1 = np.fft.fft(prices_asset1)
fft_asset2 = np.fft.fft(prices_asset2)

# Frequency comparison logic here

Математические идеи и реализации Python

Реальная сила преобразования Фурье в финансах заключается в его способности преобразовывать сложные данные во временной области (например, цены на акции) в более простую информацию в частотной области. Это преобразование позволяет трейдерам выявлять закономерности, которые не так легко увидеть во временной области.

Реализация Python

Python предлагает доступную платформу для реализации этих концепций. Для вычислений и визуализации можно использовать такие библиотеки, как numpy и matplotlib.import matplotlib.pyplot as plt

# Plotting Frequency Components
plt.plot(np.abs(fft_result))
plt.title(‘Frequency Components’)
plt.xlabel(‘Frequency’)
plt.ylabel(‘Amplitude’)
plt.show()

Дальнейший путь

Применение преобразования Фурье в финансах все еще развивается. Его способность выявлять скрытые закономерности в рыночных данных открывает новые горизонты в торговых стратегиях, особенно в хеджировании и арбитраже. Однако это требует глубокого понимания как финансовых концепций, так и математических методов. По мере того, как финансовый мир становится все более управляемым данными, такие методы, как преобразование Фурье, будут играть ключевую роль в формировании будущих торговых стратегий.

Преобразование Фурье и его применение в машинном обучении

Знакомство

Преобразование Фурье — это математический метод, который играет ключевую роль в различных научных и инженерных областях, применяясь в различных областях от обработки сигналов до квантовой механики. В последние годы он обрел новое значение в области машинного обучения. В этом эссе исследуются основы преобразования Фурье и его растущее значение в приложениях машинного обучения.

Гармонизируя мир данных, преобразование Фурье раскрывает скрытые симфонии внутри, организуя идеи, которые резонируют во всей области машинного обучения.

Общие сведения о преобразовании Фурье

Преобразование Фурье, названное в честь французского математика и физика Жана-Батиста Жозефа Фурье, представляет собой математическую операцию, которая разлагает сигнал на составляющие его частотные составляющие. Он позволяет анализировать частотное содержание сигнала и представляет его в частотной области. Это преобразование особенно полезно при работе со сложными сигналами, так как упрощает анализ лежащих в их основе паттернов.

Непрерывное преобразование Фурье (CFT) и дискретное преобразование Фурье (DFT) являются двумя распространенными вариантами. CFT используется для непрерывных сигналов, в то время как DFT применяется для дискретных сигналов, что делает его более актуальным для цифровых данных и задач машинного обучения. Быстрое преобразование Фурье (БПФ) является эффективным алгоритмом для вычисления ТПФ, что еще больше облегчает его широкое использование в различных приложениях.

Применение в обработке сигналов

Одним из наиболее традиционных и хорошо зарекомендовавших себя применений преобразования Фурье является обработка сигналов. Он используется в таких задачах, как обработка звука, анализ изображений и сжатие данных. Например, при обработке звука преобразование Фурье помогает идентифицировать различные частоты, присутствующие в аудиосигнале, позволяя выполнять такие задачи, как распознавание речи, классификация музыки и шумоподавление.

При анализе изображений преобразование Фурье можно использовать для извлечения информации о текстуре и узоре из изображений. Преобразование изображений в частотную область упрощает распознавание краев, форм и других визуальных особенностей. Это важно для таких задач, как распознавание изображений, обнаружение объектов и сжатие изображений.

Приложения в машинном обучении

Машинное обучение, область, которая фокусируется на разработке алгоритмов, способных обучаться на основе данных, все чаще использует преобразование Фурье. Его применение в этой области разнообразно и эффективно:

  1. Анализ временных рядов: В таких областях, как финансы, здравоохранение и прогнозирование погоды, данные временных рядов имеются в изобилии. Преобразование Фурье может помочь извлечь релевантные признаки из данных временных рядов, анализируя их частотные компоненты. Это имеет решающее значение для таких задач, как обнаружение аномалий, анализ тенденций и прогнозирование.
  2. Обработка естественного языка: Текстовые данные, представленные в виде последовательности слов, могут рассматриваться как дискретный сигнал. Применяя преобразование Фурье, становится возможным анализировать текстовые данные в частотной области, что находит применение в классификации текстов, анализе тональности и тематическом моделировании.
  3. Проектирование признаков: Проектирование признаков — важнейший шаг в конвейере машинного обучения. Преобразуя данные в частотную область, можно извлечь ценные признаки, которые может быть трудно зафиксировать во временной области. Это может привести к созданию более надежных и точных моделей машинного обучения.
  4. Сверточные нейронные сети (СНС): CNN — это популярная архитектура глубокого обучения для анализа изображений. Преобразование Фурье может быть использовано для создания сверточных фильтров, которые специализируются на обнаружении определенных частотных составляющих в изображениях. Это может повысить производительность СНС в таких задачах, как классификация изображений и распознавание объектов.
  5. Аугментация данных: Аугментация данных — это метод, используемый для увеличения размера обучающего набора данных. При обработке изображений преобразование Фурье может быть использовано для создания дополненных данных путем изменения частотных составляющих изображений. Это помогает улучшить обобщение и надежность моделей машинного обучения.

Код

Чтобы выполнить преобразование Фурье для анализа временных рядов в Python, вы можете использовать библиотеки numpy и matplotlib. Я предоставлю вам полный пример кода Python с использованием примера набора данных и созданием соответствующих графиков. Во-первых, вам нужно установить необходимые библиотеки, если вы еще этого не сделали:pip install numpy matplotlib

Ниже приведен код Python с образцом набора данных и графиками для анализа временных рядов с использованием преобразования Фурье:import numpy as np
import matplotlib.pyplot as plt

# Generate a sample time series dataset
# You can replace this with your own time series data
# Ensure that the data is in a NumPy array or a list
time = np.arange(0, 10, 0.01) # Time values from 0 to 10 with a step of 0.01
signal = 2 * np.sin(2 * np.pi * 1 * time) + 1 * np.sin(2 * np.pi * 2 * time)

# Plot the original time series
plt.figure(figsize=(10, 4))
plt.subplot(2, 1, 1)
plt.plot(time, signal)
plt.title(‘Original Time Series’)
plt.xlabel(‘Time’)
plt.ylabel(‘Amplitude’)

# Perform the Fourier Transform
fourier_transform = np.fft.fft(signal)
frequencies = np.fft.fftfreq(len(signal), 0.01) # Frequency values (assuming a sampling interval of 0.01)

# Plot the magnitude of the Fourier Transform
plt.subplot(2, 1, 2)
plt.plot(frequencies, np.abs(fourier_transform))
plt.title(‘Fourier Transform’)
plt.xlabel(‘Frequency (Hz)’)
plt.ylabel(‘Magnitude’)
plt.xlim(0, 5) # Limit the x-axis to show frequencies up to 5 Hz

plt.tight_layout()
plt.show()

В этом коде:

  1. Мы создаем образец набора данных временных рядов с двумя синусоидальными компонентами. Вы должны заменить их своими собственными данными временных рядов.
  2. Мы используем np.fft.fft для выполнения преобразования Фурье для данных временных рядов.
  3. Мы вычисляем соответствующие частоты с помощью np.fft.fftfreq.
  4. Мы создадим два подграфика: один для исходного временного ряда, а другой для величины преобразования Фурье.
  5. Наконец, мы выводим графики с помощью plt.show().

При работе с реальными данными обязательно замените образец набора данных собственными данными временных рядов. Приведенный выше код даст вам график исходного временного ряда и величины преобразования Фурье. Вы можете настроить параметры печати и метки в соответствии с вашими конкретными потребностями.

Заключение

Преобразование Фурье, первоначально разработанное для обработки сигналов, нашло свое применение в самом сердце машинного обучения. Его способность анализировать и извлекать информацию о частоте из данных неоценима в различных приложениях машинного обучения, включая анализ временных рядов, обработку естественного языка, конструирование признаков и улучшение моделей глубокого обучения. По мере того, как машинное обучение продолжает развиваться и расширять свои горизонты, преобразование Фурье остается мощным инструментом, позволяющим исследователям и практикам получать более глубокие знания и разрабатывать более эффективные алгоритмы. Его универсальность и адаптивность делают его незаменимым компонентом современного инструментария машинного обучения.

Роль преобразования Фурье в проектировании признаков

Знакомство

Преобразование Фурье является мощным инструментом в конструировании признаков, широко используемым в различных областях, таких как обработка сигналов, анализ изображений и наука о данных. Его значение заключается в преобразовании сигналов во времени или пространстве в частотную область, предлагая другую перспективу для анализа и обработки данных. В этом эссе мы рассмотрим концепции преобразования Фурье и его применение в проектировании признаков.

В потоке данных преобразование Фурье является маяком, раскрывающим скрытые частоты озарения.

Общие сведения о преобразовании Фурье

По своей сути преобразование Фурье разлагает функцию времени (сигнал) на составляющие его частоты. Это очень важно, потому что во многих практических сценариях анализ частотных составляющих сигнала может быть более проницательным, чем изучение сигнала в его исходной временной области. Например, при обработке сигналов он может выявлять скрытые периодичности или доминирующие частоты, которые не очевидны во временной области.

Применение в проектировании признаков

Конструирование признаков — это процесс использования знаний предметной области для извлечения признаков из необработанных данных. Эти особенности обеспечивают эффективную работу алгоритмов машинного обучения. Преобразование Фурье особенно полезно в этом контексте по нескольким причинам:

  1. Шумоподавление: В реальных данных шум неизбежен. Преобразование Фурье может помочь изолировать шумовую составляющую от сигнала, облегчая уменьшение или устранение шума, тем самым повышая качество данных.
  2. Извлечение признаков: В данных временных рядов определенные закономерности могут быть различимы только в частотной области. Преобразуя данные, можно обнаружить новые признаки, которые могут иметь решающее значение для прогнозного моделирования.
  3. Сжатие данных: Преобразование Фурье может значительно уменьшить размерность данных без потери критически важной информации. Это особенно полезно при обработке изображений и аудиосигналов, где уменьшение размера данных может привести к более эффективному хранению и обработке.
  4. Гомогенизация вариаций шкалы времени: В наборах данных, где присутствуют вариации шкалы времени (например, два похожих сигнала, но один растянут во времени), преобразование Фурье может помочь нормализовать эти сигналы до общего масштаба.

Практические примеры

На практике преобразование Фурье применяется различными способами. Например, при обработке аудиосигнала он используется для идентификации различных музыкальных нот в композиции. При обработке изображений он помогает в распознавании краев и фильтрации изображений. В финансах преобразования Фурье используются для анализа данных временных рядов для выявления тенденций и циклического поведения в данных фондового рынка.

Код

Чтобы продемонстрировать применение преобразования Фурье в конструировании признаков с помощью Python, мы создадим синтетический набор данных, применим к нему преобразование Фурье, а затем визуализируем результаты. Наш синтетический набор данных будет состоять из комбинации синусоидальных сигналов, что является распространенным сценарием при обработке сигналов.

Пошаговый процесс:

  1. Generate Synthetic Data (Генерация синтетических данных): Создание данных временных рядов, состоящих из синусоидальных сигналов.
  2. Применение преобразования Фурье: Используйте быстрое преобразование Фурье (БПФ) для преобразования данных временных рядов в частотную область.
  3. Визуализация результатов: построение графика как исходных данных временных рядов, так и их представления в частотной области.

Код на Python:import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import fft

# Time array
t = np.linspace(0, 1, 400, endpoint=False)

# Generating a synthetic signal: combination of two sine waves
signal = np.sin(2 * np.pi * 5 * t) + np.sin(2 * np.pi * 10 * t)

# Plotting the synthetic signal
plt.figure(figsize=(12, 6))
plt.plot(t, signal)
plt.title(«Synthetic Time Series Data»)
plt.xlabel(«Time»)
plt.ylabel(«Amplitude»)
plt.show()

# Applying Fourier Transform using FFT
signal_fft = fft(signal)

# Getting frequencies
n = len(signal)
freq = np.fft.fftfreq(n, d=(t[1]-t[0]))

# Plotting the Fourier Transform (Spectrum)
plt.figure(figsize=(12, 6))
plt.plot(freq, np.abs(signal_fft))
plt.title(«Fourier Transform of the Synthetic Signal»)
plt.xlabel(«Frequency (Hz)»)
plt.ylabel(«Amplitude»)
plt.xlim([0, 50]) # Limiting to show only positive frequencies
plt.show()

Этот код будет создавать графики как для сигнала во временной области, так и для его представления в частотной области. На первом графике показан синтетический сигнал, представляющий собой комбинацию двух синусоидальных волн. Второй график, преобразование Фурье, показывает частотные составляющие этого сигнала, показывая пики на частотах, соответствующих добавленным синусоидальным волнам (5 Гц и 10 Гц).

Вы можете запустить этот код в среде Python, чтобы визуализировать, как преобразование Фурье может быть использовано в конструировании признаков для извлечения признаков из данных временных рядов. Помните, что это упрощенный пример для иллюстрации. В реальных сценариях интерпретация и применение преобразования Фурье может потребовать более сложной обработки данных и соображений, специфичных для конкретной предметной области.

Ограничения и рекомендации

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

Заключение

Преобразование Фурье — это универсальный и мощный инструмент в области проектирования признаков, предлагающий другую призму для просмотра и обработки данных. Его способность преобразовывать данные в частотную область открывает множество возможностей для шумоподавления, извлечения признаков, сжатия данных и многого другого. Понимание и применение этого метода может значительно повысить качество и эффективность анализа данных и моделей машинного обучения. Тем не менее, важно осознавать его ограничения и необходимость знаний в предметной области для эффективной интерпретации результатов.

Источник

Источник

Источник

Источник

Источник

Источник