Оценка L1-нормы полусверток гауссовских плотностей

[_hum_]

Определим для любых абсолютно-интегрируемых функций $f_1$ u $f_2$ операцию "полусвертки" $[f_1\, @ f_2](t) = \int_{0}^{\infty}f_1(s)f_2(t-s)ds$ (отличие от свертки - нижний предел интегрирования). Договоримся также считать эту операцию левоассоциативной и обозначать $\underset{\text{n}}{\underbrace{f\,@\,\cdots\,@\,f\,}}$ через $f^{@ n}$.
Пусть $\varphi_u(t) = \frac{1}{\sqrt{2\pi}}e^{-\frac{(x-u)^2}{2}}$. Для случая, когда $u < 0$ интересует скорость убывания по $n$ величин:

$$\int_{0}^{\infty}\varphi_u^{@ n}$$

(а в идеале, величин, когда нижний предел интегрирования - произвольное неотрицательное число).

С одной стороны можно было воспользоваться тем, что полусвертка неотрицательных функций не больше свертки оных, то есть, $\varphi_u^{@ n} \leq \varphi_u^{* n} = \varphi_{n u,\sqrt{n}}$, но кажется, что получится слишком грубая оценка. Для нахождения более тонкой думал было обратиться к Young's convolution inequality. Но как-то с наскоку не удалось им воспользоваться. Вот и хотелось спросить у знатоков, имеет ли смысл вообще в эту сторону копать?
Спасибо.

[DeBill]

_hum_
Пусть $p_n$ - та самая величина.
1. Рассмотрим случайное блуждание на прямой, когда величина каждого шага - случайная, с плотностью $\varphi =\varphi_u$. Тогда $p_n$ - (вроде бы) вероятность того, что после $n$ шагов точка осталась на луче $[0,\infty)$. Может, что то из теории таких блужданий нам поможет?
2. (опять же, если я не ошибся при замене переменной) Имеем:
$p_n=\int\limits_{U_n}^{}\varphi(x_1)...\varphi(x_n)dx$,
где $U_n=\{x_1>0,x_1+x_2>0,...,x_1+...+x_n>0\}$.
Рассмотрим еще интегралы
$q_n=\int\limits_{V_n}^{}\varphi(x_1)...\varphi(x_n)dx$,
где $V_n=\{x_1+...+x_n>0\}$. Эти интегралы стандартной заменой сводятся к функции Лапласа (или: $q_n$ есть вероятность того, что сумма $n$ независимых, с плотностью $\varphi$, больше 0; т.к., сумма нормальных - нормальна, то ответ выражается через $\Phi$), и асимптотику $q_n$ можно получить из известной асимптотики для нее. Осталось установить связь асимптотик для $p_n$ и $q_n$. И тут будут только некие правдоподобные рассуденя.
Линии уровня $L_C$ подынтегральной функции $f(x)=\varphi(x_1)...\varphi(x_n)$ -сферы с центром в точке $(u,u,...,u)$. В частности, сфера $L_c$ для $c=f(0)$ касается плоскости $S=\partial V_n=\{x_1+...+x_n=0\}$. Поэтому, для малых $\varepsilon>0$ кусок сферы $L_{c+\varepsilon}$, попавший в $V_n$, похож на малый диск в $S$. Аналогично, для малых $\varepsilon>0$ кусочек сферы $L_{c+\varepsilon}$, попавший в $U_n$, похож на часть этого малого диска. Заметим, что (см. ниже) этот кусочек в $n$ раз меньше куска (почти). Поэтому, в асимптотике $p_n$ по сравнению с $q_n$ появится дополнительный множитель $\frac{1}{n}$ ...Всё.
3. Пусть $B$ - единичный шар в $\mathbb{R}^n$, и $S=\{x_1+...+x_n=0\}$ - гиперплоскость. Какую часть от "площади" диска $b=B\cap S$ составляют точки, для которых $x_1>0,x_1+x_2>0,...,x_1+...+x_{n-1}>0$? Ответ: $\frac{1}{n}$
Я это отношение сосчитал ручками для $n=2,3$ и 4; получился красивый ответ. "Эта "ж"- неспроста", решил я. И правда, получилось.
Выложу эту задачку в Олимпиадном разделе.

[_hum_]

DeBill в сообщении #1545986 писал(а):

1. Рассмотрим случайное блуждание на прямой, когда величина каждого шага - случайная, с плотностью $\varphi =\varphi_u$. Тогда $p_n$ - (вроде бы) вероятность того, что после $n$ шагов точка осталась на луче $[0,\infty)$. Может, что то из теории таких блужданий нам поможет?

почему-то мне не удается найти подобного рода блуждания - всюду либо с дискретным шагом, либо с непрерывным временем :(

насчет связи интегралов $p_n$ и $q_n$ мысль не уловил - какую роль тут играет "малое $\varepsilon$", и почему

DeBill в сообщении #1545986 писал(а):

ля малых $\varepsilon>0$ кусок сферы $L_{c+\varepsilon}$, попавший в $V_n$, похож на малый диск в $S$

?

Вот рисунок границ соответствующих областей для $n=3$ (оранжевая - гиперплоскость):

[ipgmvq]

DeBill в сообщении #1545986 писал(а):

3. Пусть $B$ - единичный шар в $\mathbb{R}^n$, и $S=\{x_1+...+x_n=0\}$ - гиперплоскость. Какую часть от "площади" диска $b=B\cap S$ составляют точки, для которых $x_1>0,x_1+x_2>0,...,x_1+...+x_{n-1}>0$? Ответ: $\frac{1}{n}$
Я это отношение сосчитал ручками для $n=2,3$ и 4; получился красивый ответ. "Эта "ж"- неспроста", решил я. И правда, получилось.
Выложу эту задачку в Олимпиадном разделе.

У меня эмпирически $\frac{1}{n}$ не получается

Где-то ошибка или я задачу неправильно понял?...

Используется синтаксис Python

import numpy as np
rng = np.random.default_rng()
повторов = 1000000
for измерений in range(2,18):
    выборка = rng.uniform(size=(повторов,измерений)) * 2.0 - 1.0
    фильтр = np.power(выборка, 2).sum(axis=1) <= 1.0
    выборка = выборка[фильтр,:]
    нормы = выборка.cumsum(axis=1)[:,:-1]
    условие = np.all(нормы > 0.0, axis=1)
    print("Измерений: %i, выборка: %i, доля: %.3f%%" % (измерений, выборка.shape[0], 100.0 * условие.mean()))
 

[DeBill]

_hum_ в сообщении #1546201 писал(а):

насчет связи интегралов $p_n$ и $q_n$

Перейдем в интегралах к сферической системе координат (с центром в точке $(u,u,...,u)$). Тогда
$q_n=(2\pi)^{-\frac{n}{2}}\int\limits_{-u\sqrt{n}}^{\infty}e^{-\frac{1}{2}r^2}\sigma_n(r)dr$,
где $\sigma_n(r) $ - "площадь" части сферы радиуса $r$, попавшей в $V_n$. Главный член асимптотики будем считать методом Лапласа (см, напр., Зорич, т.2). Из формул для асимптотического разложения (придется их использовать, поскоку коэф-т $\sigma_n$ при экспоненте зануляется в начале отрезка интегрирования) видим: асимптотика определяется показателем экспоненты, и главным членом асимптотики для $\sigma_n$. То же верно для $p_n$, с функцией $\widetilde{\sigma_n}$.
Но, по п.3., $\lim\limits_{r\to -u\sqrt{n}+0}\frac{\sigma_n(r)}{\widetilde{\sigma_n(r)}}=n$ . Значит, таким же будет и отношение главных членов асимптотики для $q_n$ и $p_n$.
Ну, что то вроде этого....

[_hum_]

DeBill
кажется, начинаю понимать ход Ваших мыслей. Но тогда вопрос, если п.3 справедлив для любых шаров, то, получается, должно выполняться равенство $q_n = p_n n$, ведь $q_n=(2\pi)^{-\frac{n}{2}}\int e^{-\frac{1}{2}r^2}\sigma_n(r)dr = (2\pi)^{-\frac{n}{2}}\int e^{-\frac{1}{2}r^2}n \widetilde{\sigma_n}(r)dr = n p_n$. Нет?

ipgmvq
к сожалению, я не знаком с numpy, но по коду у Вас как-то подозрительно присутствует rng.uniform, намекающее на равномерное распределение вместо нормального.

[ipgmvq]

_hum_ в сообщении #1546232 писал(а):

к сожалению, я не знаком с numpy, но по коду у Вас как-то подозрительно присутствует rng.uniform, намекающее на равномерное распределение вместо нормального.

Да, это потому, что в пункте 3 был единичный шар.
Но я неправильно посчитал в том посте. После пересчета в симуляции действительно получается очень близко к $\frac{1}{n}$.
Тут распечатка результатов.
И новый код:

Используется синтаксис Python

import numpy as np
rng = np.random.default_rng()
повторов = 1000000
for измерений in range(2,15):
    выборка = rng.uniform(size=(повторов,измерений-1)) * 2.0 - 1.0
    доп_измерение = (-выборка.sum(axis=1)).reshape((повторов, 1))
    выборка = np.hstack((выборка, доп_измерение))
    фильтр_L2 = np.power(выборка, 2).sum(axis=1) <= 1.0
    выборка = выборка[фильтр_L2,:]
    случайный_сдвиг_измерений = rng.integers(измерений, size=(повторов,))
    нормы_L1 = выборка.cumsum(axis=1)[:,:-1]
    условие = np.all(нормы_L1 > 0.0, axis=1)
    print("Измерений: %i,\tвыборка: %i,\tдоля: %.3f%%" % (измерений, выборка.shape[0], 100.0 * условие.mean()))
    выборка = np.roll(выборка, случайный_сдвиг_измерений, axis=1)
    нормы_L1 = выборка.cumsum(axis=1)[:,:-1]
    условие = np.all(нормы_L1 > 0.0, axis=1)
    print("Измерений: %i,\tвыборка: %i,\tдоля: %.3f%%, доп. случайность" % (измерений, выборка.shape[0], 100.0 * условие.mean()))
 

Задача на глаз похожа на то, как будто Вы хотите что-то оптимизировать на плотности распределения лог-доходностей портфеля out-of-the-money колл-опционов (без учета их стоимости и по одному страйку на акцию) для очень высокодивидентных акций.

[DeBill]

_hum_ в сообщении #1546232 писал(а):

если п.3 справедлив для любых шаров

Это вроде так.

_hum_ в сообщении #1546232 писал(а):

получается, должно выполняться равенство $q_n = p_n n$

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

[ipgmvq]

_hum_
Сделал симуляцию по Вашей спецификации (надеюсь без ошибок на сей раз).
Финальный интеграл от n судя по ней не зависит и при фиксированной $u \leqslant 0$ является константой.

Вот к примеру зависимость финального интеграла от n при $u = 0$:

Оценка интеграла с увеличением n идёт в разнос из-за прогрессивного уменьшения эффективного размера выборки.

Поэтому я взял $n = 2$ и получил зависимость финального интеграла от $u$:


Мой код симуляции на GitHub вместе с таблицами зависимостей и этими графиками.

И он же ниже:

Используется синтаксис Python

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
rng = np.random.default_rng()
pd.options.display.float_format = '{:6.4f}'.format

sample_size = 1000000

sample_sizes = list()
inlegrals = list()
μ = 0.0
n_dimentions_list = range(2,15)
for n_dimentions in n_dimentions_list:
    sample = rng.normal(loc=μ, size=(sample_size, n_dimentions))
    norm_L1 = sample.cumsum(axis=1)
    norm_L1[:,-1] = sample[:,-1]
    filter_L1 = np.all(norm_L1 > 0.0, axis=1)
    sample_size_temp = filter_L1.sum()
    integral_temp = sample[filter_L1,-1].mean()
    sample_sizes.append(sample_size_temp)
    inlegrals.append(integral_temp)
relationship_to_dimentions = pd.DataFrame({'Sample size': sample_sizes, 'Integral': inlegrals}, index = n_dimentions_list)
relationship_to_dimentions.index.name = '# dimentions'

print(relationship_to_dimentions)

the_range = relationship_to_dimentions.Integral.max() - relationship_to_dimentions.Integral.min()
relationship_to_dimentions.Integral.plot(lw=2, 
                                         colormap='jet', 
                                         marker='.', 
                                         markersize=7, 
                                         title='Relationship semi-convolution integral to number of dimentions', 
                                         xlabel = 'Number of dimentions', ylabel = 'Positive integral', 
                                         ylim=(relationship_to_dimentions.Integral.mean() - 3.0 * the_range, 
                                               relationship_to_dimentions.Integral.mean() + 3.0 * the_range))
plt.show()

sample_sizes = list()
inlegrals = list()
n_dimentions = 2
μs = np.arange(-2.0,0.00001,0.02)
for μ in μs:
    sample = rng.normal(loc=μ, size=(sample_size, n_dimentions))
    norm_L1 = sample.cumsum(axis=1)
    norm_L1[:,-1] = sample[:,-1]
    filter_L1 = np.all(norm_L1 > 0.0, axis=1)
    sample_size_temp = filter_L1.sum()
    integral_temp = sample[filter_L1,-1].mean()
    sample_sizes.append(sample_size_temp)
    inlegrals.append(integral_temp)
relationship_to_μ = pd.DataFrame({'Sample size': sample_sizes, 'Integral': inlegrals}, index = μs)
relationship_to_μ.index.name = 'μ'

relationship_to_μ.index = relationship_to_μ.index.format(formatter='{:5.2f}'.format)
print(relationship_to_μ)

relationship_to_μ.Integral.plot(lw=2, 
                                colormap='jet', 
                                marker='.', 
                                markersize=7, 
                                title='Relationship semi-convolution integral to μ', 
                                xlabel = 'μ', 
                                ylabel = 'Positive integral')
plt.show()
 

А Вы програмируете на C, "C++" или R?

[ipgmvq]

Увы, маштабирование в расчете сделал неверное.
Переделал: код на GitHub
Конечно от n финальный интеграл зависит.

Используется синтаксис Python

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
rng = np.random.default_rng()
pd.options.display.float_format = '{:6.4f}'.format

sample_size = 1000000

sample_sizes = list()
inlegrals = list()
μ_list = list()
dimentions_list = list()
n_dimentions_list = np.arange(2,15)
μs = np.arange(-2.0,0.00001,0.05)

for n_dimentions in n_dimentions_list:
    for μ in μs:
        sample = rng.normal(loc=μ, size=(sample_size, n_dimentions))
        norm_L1 = sample.cumsum(axis=1)
        norm_L1[:,-1] = sample[:,-1]
        filter_L1 = np.all(norm_L1 > 0.0, axis=1)
        sample_size_temp = filter_L1.sum()
        sample[np.logical_not(filter_L1),-1] = 0.0
        integral_temp = sample[:,-1].mean()
        μ_list.append(μ)
        dimentions_list.append(n_dimentions)
        sample_sizes.append(sample_size_temp)
        inlegrals.append(integral_temp)
relationship = pd.DataFrame({'Sample size': sample_sizes, 'Integral': inlegrals, '# dimentions': dimentions_list, 'μ': μ_list})

print(relationship)

fig = plt.figure(figsize=(14,14))
ax = fig.add_subplot(111, projection='3d')
ax.plot_trisurf(
                relationship['# dimentions'], 
                relationship['μ'], 
                relationship['Integral'],
                cmap=cm.gist_ncar)

ax.view_init(20, -115)
ax.set_xlabel('Number of dimentions')
ax.set_ylabel('μ')
ax.set_zlabel('Integral')
plt.show()
 

[ipgmvq]

Увы, опять посчитал не то, что было Вам нужно.
Переделал опять: код на GitHub

Используется синтаксис Python

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
import networkx as nx
from matplotlib.animation import FuncAnimation, PillowWriter
rng = np.random.default_rng()
pd.options.display.float_format = '{:6.4f}'.format

sample_size = 1000000

integrals = list()
μ_list = list()
dimentions_list = list()
n_dimentions_list = np.arange(2,15)
μs = np.arange(-2.0,0.00001,0.05)
for n_dimentions in n_dimentions_list:
    for μ in μs:
        sample = rng.normal(loc=μ, size=(sample_size, n_dimentions))
        norm_L1 = sample.cumsum(axis=1)
        filter_L1 = np.all(norm_L1 > 0.0, axis=1)
        μ_list.append(μ)
        dimentions_list.append(n_dimentions)
        integrals.append(filter_L1.mean())
relationship = pd.DataFrame({'Integral': integrals, '# dimentions': dimentions_list, 'μ': μ_list})

print(relationship)

fig = plt.figure(figsize=(14,14))
ax = fig.add_subplot(111, projection='3d')
ax.plot_trisurf(
                relationship['# dimentions'], 
                relationship['μ'], 
                relationship['Integral'],
                cmap=cm.gist_ncar)

ax.view_init(20, -115)

ax.set_xlabel('Number of dimentions')
ax.set_ylabel('μ')
ax.set_zlabel('Integral')
plt.show()
 

[ipgmvq]

ipgmvq в сообщении #1546246 писал(а):

Задача на глаз похожа на то, как будто Вы хотите что-то оптимизировать на плотности распределения лог-доходностей портфеля out-of-the-money колл-опционов (без учета их стоимости и по одному страйку на акцию) для очень высокодивидентных акций.

И я понял, что это никакого отношения к портфелю не имеет, но всего ли поиск зависимости распределения supremum'а на дискретном во времени временном ряду от такого его параметра как длина временного ряда.
Эта публикация, мне кажется, как раз по теме (в т.ч. с негативным дрифтом, как Вы писали).
Choe J, Shroff NB. On the supremum distribution of integrated stationary Gaussian processes with negative linear drift. Advances in Applied Probability. 1999;31(1):135-57.
И глава 10-я (особенно 10.2) этого учебника тоже в тему.
Falk M, Hüsler J, Reiss RD. Extremes of gaussian processes. In: Laws of Small Numbers: Extremes and Rare Events. Springer, Basel; 2011:381-418. DOI: 10.1007/978-3-0348-0009-9_10

То, что интересует именно дискретное время, странно, потому что для вероятности stop loss или take profit по акции подошло бы именно непрерывное. А это больше подходит к структурным нотам по единственной акции, у которых купонные события идут раз в квартал...

[DeBill]

ipgmvq
А не могли бы Вы проверить численно полученную формулу для $p_n$$
$p_n\sim \frac{e^{-\frac{nu^2}{2}}}{\sqrt{2\pi n}\left\lvert u \right\rvert n$}$
(она должна работать при больших $n$, но также и при малых $n$ и больших отрицательных $u$).
Здесь $p_n$ - вероятность того, что за $n$ шагов случайного блуждания точка не покинет область $[0,\infty)$ (шаги - независимы, нормально распределены с отрицательным матожиданием $u$ и единичной дисперсией). Это можно считать без кратных интегралов, простой монтекарлой: интересует совпадение с указанной асимптотикой при больших по модулю $u$.

[ipgmvq]

DeBill
Т.к. тут в симуляции ищется пропорция (а не средняя), на очень низких значениях вероятности для сохранения разрешения требуются гиганские выборки, иначе в симуляции Монте Карло при нарастании $n$ и падении $\mu$ мы быстро начинаем получать пропорцию ноль из десяти миллионов, ноль из ста миллионов. При нарастании выборки время симуляции превращается из минут в дни. Это задача не для Питона, но для Си и распараллеливания симуляции по десяткам ядер ЦПУ.

На относительно небольших $n$ и не очень низких $\mu$ мы имеет такое соотношение симуляции (зеленый) и приведенной формулы (фиолетовый).

Код на GitHub.

[DeBill]

ipgmvq
Спасибо! Красивые картинки...
И, похоже (судя по первой), есть стабильное отставание - в лог-масштабе. Что, возможно, говорит о неправильной константе в формуле. А может, и нет - формула таки лишь асимптотическая.